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

using Markdown

# ╔═╡ 0d3bf5f6-1171-11ec-0fee-c73bb459dc3d
begin
	using HierarchicalTemporalMemory
	using Random, Chain, Setfield, Statistics, Test
end

# ╔═╡ 1bb0fcfc-2d7a-4634-9c93-263050c56a55
md"""
# Projection

This notebook implements the simplest operation of assembly calculus, *projection*, with HTM.

Assembly calculus and HTM share the concept of "Region" as a modular circuit of neurons.
HTM imposes specific structure in the Region.
The question is: _is the HTM structure behaving according to assembly calculus_?

Projection of assembly `x` from region A -> B is the operation where a stable assembly `y` is generated in B that fires consistently when `x` fires in A.
In assembly calculus this is shown to happen simply by repeated firing of `x` for a number of time steps:

```algo
repeat T:    # T: constant ~10
  fire x
```

Alternatively (to not require a constant T):

```algo
repeat until yₜ converges:
  fire x
```

Note that we don't necessarily expect point stability, but a limit circle of low enough radius. More on this with the convergence measure.

## Experiment

### The Plan

We will implement this simple algorithm with HTM regions and check the convergence of `y`.
If `yₜ` converges to a limit circle as $t→∞$ (practically 10-100 steps) of a small enough radius, we will consider the projection implemented.

But how small is small enough?
The algorithm parameters might have a big impact on the convergence radius, or might even prohibit convergence.
Intuitively in order to set an expectation, let the limit radius $|yₜ₁ - yₜ₂| < 5\%|y|$.

Let's define the relative distance behind this limit radius.
Since $yᵢ$ are sparse binary vectors:
- the core "similarity" kernel will be bitwise AND: $yᵢ .\& yⱼ$
- the norm $|yᵢ|$ will be the count of `true`

The relative distance $\frac{|yᵢ - yⱼ|}{|y|}$ can be:
"""

# ╔═╡ 70471bdc-660e-442e-b92a-f486abd120a7
reldistance(yᵢ,yⱼ)= 1 - count(yᵢ .& yⱼ) / min(count(yᵢ), count(yⱼ))

# ╔═╡ 7620c202-c6cf-44db-9cda-13642e28b45a
md"""
which is bounded by $0 ≤ \texttt{reldistance}(yᵢ, yⱼ) ≤ 1$
"""

# ╔═╡ 043f8da8-f17a-4486-9c49-7d5ecc75457f
md"""
### Collect the ingredients

Let's define the HTM regions A, B with the same properties (adapting the input size to match).
Minicolumns will be explained later.
"""

# ╔═╡ bbaf0a64-2471-4106-a27c-9dd5a4f58c7c
begin
	Nin= 1e3|> Int   # input size
	Nn= 4e5          # number of neurons in each area
	k= 10            # neurons per minicolumn
	_Nc()= floor(Int,Nn/k) # number of minicolumns
	sparsity()= 1/sqrt(_Nc());
end

# ╔═╡ d4683829-b214-4d79-890d-bcc902bcef2b
params_A= (
	sp= SPParams(szᵢₙ=Nin, szₛₚ=_Nc(), s= sparsity(), prob_synapse=1e-3, enable_local_inhibit=false),
	tm= TMParams(Nc=_Nc(), k=k, p⁺_01= .11, θ_stimulus_learn=15, θ_stimulus_activate= 22)
)

# ╔═╡ 5b0c7174-dc2b-4b28-8cda-fc2e8a37151c
params_B= @set params_A.sp.szᵢₙ= params_A.tm.Nₙ

# ╔═╡ 334b271d-247c-465c-ae82-f91007a6d9d0
md"The regions:"

# ╔═╡ 4833f12f-3eac-407c-a9ce-0d1aea216077
A= Region(params_A.sp, params_A.tm);

# ╔═╡ 45d8e50b-8baa-4da2-8cd0-5c81d634c450
B= Region(params_B.sp, params_B.tm);

# ╔═╡ 91fa8319-1702-4cad-8ae8-0c1ed51a7bb8
md"""
This is the assembly `x`, which first needs a random input.
"""

# ╔═╡ 3ee589ed-0fe0-4328-88e8-d2fb6a98bf28
x= @chain bitrand(Nin) A(_).active;

# ╔═╡ 6790a86a-3990-4c6a-878f-c18779f8d48d
md"""
Now all we need to do is to repeatedly stimulate `B` with `x` and let it learn.
Each time `B` is stimulated, it will produce a `yᵢ`:
"""

# ╔═╡ 5458c6ae-ace2-4928-8a9b-ef919a1e97cb
y₁= B(x).active;

# ╔═╡ 6b56422c-5f9b-41d8-a080-ef5ca3d7db7b
md"""
But in order to let region `B` learn from `x`, we need to call the `step!` function, which will return the same output as `yᵢ` the first time it's called, while at the same time adapting synapses:
"""

# ╔═╡ 3229da2a-92f3-4064-bb80-cbd9f5523d7d
begin
	reset!(B)
	md"`step!(B,x).active == y₁`: $((step!(B,x).active .== y₁) |> all)"
end

# ╔═╡ 5e0ac810-1689-4bfb-88a8-85504cd821b6
md"""
### Repeated stimulation

To plot the convergence of $yᵢ$ we will now produce them repeatedly and store them.
Each experiment will run for `T` steps and we will run 15 experiments with different `x`.
"""

# ╔═╡ 1015f629-9818-4dbb-b574-97c552e96164
T= 30; experiments=10;

# ╔═╡ 62d41be1-2970-48e1-b689-c2ecf1ab10ce
md"`y[T,experiments]` is the entire history of `yᵢ`"

# ╔═╡ 7670ffbc-3621-4acb-b1e3-b56668470c24
y= Matrix(undef,T,experiments);

# ╔═╡ d08cce9d-dd9c-4530-82ae-bf1db8be3281
# Send SDRs from A -> B
for e= 1:experiments
	@info "experiment $e"
	x= @chain bitrand(Nin) A(_).active;
	reset!(B)
	for t= 1:T
	t%10==0 && @info "t=$t"
	y[t,e]= step!(B,x).active
	end
end

# ╔═╡ 6a20715a-9dc9-461a-b6a2-f02522e277f0
md"We calculate the $Δyₜ=\texttt{reldistance}(yₜ,yₜ₊₁) ∀ t∈\{1..T-1\}$ and the experiment median. This metric will show whether `yₜ` converges."

# ╔═╡ e63a5860-587c-4aff-91be-a894b3443943
Δyₜₑ()= [reldistance(y[i,e], y[i+1,e]) for i=1:T-1, e=1:experiments];

# ╔═╡ c866c71e-5dce-4af1-84bb-045815d1acb9
expmedian(y)= median(y, dims=2);

# ╔═╡ 75560584-1631-4c93-8808-e269d345e1c8
Δȳₜ= Δyₜₑ()|> expmedian;

# ╔═╡ 456021aa-3dbf-4102-932f-43109905f9eb
begin
	# plot(Δyₜₑ(), minorgrid=true, ylim= [0,.3], opacity=0.3, labels=nothing,
	# 	title="Convergence of yₜ", ylabel="Relative distance per step Δyₜ",
	# 	xlabel="t")
	# plot!(Δȳₜ, linewidth=2, label="median Δȳ")
	# hline!(zeros(T) .+ .05, linecolor=:black, linestyle=:dash, label="5% limit")
	@test all(Δȳₜ[25:end] .< .05)
end

# ╔═╡ 56e0da0b-8858-459a-9aae-0eba4797cc06
md"""
## Understanding the convergence plot

As we see from the median (thick line), the experiments **converge after step $t=20$**.
This is double what assembly theory expects (10 steps), but it might be sensitive to the model's parameters.

So far, it looks like the "projection" operation has been implemented with HTM successfully.

#### Why is the distance exactly 0 in the beginning?

_Distal and proximal synapses_

This has to do with the biggest deviation of HTM from assembly calculus.
In assembly calculus, all synapses can trigger neurons to fire.
In HTM there are 2 kinds of synapses: _proximal_ and _distal_.
Only the proximal ones cause neurons to fire; distal synapses only cause neurons to win in competition among many that are proximally excited (keyword: *context*).

Neurons are grouped in "minicolumns", wherein they share the same proximal synapses, differing only in distal synapses.
Feedforward connectivity from A → B is through proximal synapses, while recurrent connectivity within B is through distal.
Every time `x` fires the same large set of neurons (minicolumns) is primed to fire in B.
In the beginning, before recurrent connections form through learning, every neuron in the primed minicolumns fires.
Gradually, recurrent connections form, which cause only a few neurons (often 1) in each minicolumn to be "excited" through distal connections.
But this distal excitement doesn't cause neurons to fire;
only, if the minicolumn happens to be activated in the following step, the distally excited neuron will win an in-column competition and be the only one to fire.

If this hypothesis is correct, then we expect to see a significant drop in the number of active neurons in B after the first few steps. Let's test this.
"""

# ╔═╡ cf83549c-03b1-4f7d-be29-d3a29da3e8f7
activity_n(y)= count.(y)|> expmedian

# ╔═╡ 85dfe3ed-cf70-4fde-a94a-f70c3fe4abf6
begin
	# using Plots.PlotMeasures
	# plot(activity_n(y)/activity_n(y)[1], linecolor=:blue, foreground_color_subplot=:blue, ylabel="Relative activity in B", rightmargin=14mm, legend=:none)
	# plot!(twinx(), Δȳₜ, linecolor=:red, foreground_color_subplot=:red, ylabel="median Δȳ", legend=:none)
end

# ╔═╡ 2b6d959f-63bb-4d42-a34d-70d93f6a1636
md"""
The difference in activity is much less than `k` (neurons per minicolumn), but it is still evident.
This means that after convergence there are still many neurons firing in each minicolumn;
in fact, the average number of neurons firing in each active minicolumn should be close to the region's activity in the previous plot.
"""

# ╔═╡ 3a5e3aa0-e34e-457c-a405-9bc95cd88910
activity_per_active_minicolumn(y)= @chain y reshape(_,k,:) count(_, dims=1) filter(x->x>0, _) mean

# ╔═╡ 5ec611f4-d3ea-4c5d-9221-4dcf2da3e18c
@test sum(abs.(expmedian(activity_per_active_minicolumn.(y)) .- activity_n(y)/activity_n(y)[1]*k)) < 1e-10

# ╔═╡ 3df40867-674b-405a-a05c-0733a4caaa67
md"The 2 are almost identical, therefore this explanation of the first few steps in the graph seems to hold up."

# ╔═╡ Cell order:
# ╟─0d3bf5f6-1171-11ec-0fee-c73bb459dc3d
# ╟─1bb0fcfc-2d7a-4634-9c93-263050c56a55
# ╠═70471bdc-660e-442e-b92a-f486abd120a7
# ╟─7620c202-c6cf-44db-9cda-13642e28b45a
# ╠═043f8da8-f17a-4486-9c49-7d5ecc75457f
# ╠═bbaf0a64-2471-4106-a27c-9dd5a4f58c7c
# ╠═d4683829-b214-4d79-890d-bcc902bcef2b
# ╠═5b0c7174-dc2b-4b28-8cda-fc2e8a37151c
# ╟─334b271d-247c-465c-ae82-f91007a6d9d0
# ╠═4833f12f-3eac-407c-a9ce-0d1aea216077
# ╠═45d8e50b-8baa-4da2-8cd0-5c81d634c450
# ╟─91fa8319-1702-4cad-8ae8-0c1ed51a7bb8
# ╠═3ee589ed-0fe0-4328-88e8-d2fb6a98bf28
# ╟─6790a86a-3990-4c6a-878f-c18779f8d48d
# ╠═5458c6ae-ace2-4928-8a9b-ef919a1e97cb
# ╟─6b56422c-5f9b-41d8-a080-ef5ca3d7db7b
# ╠═3229da2a-92f3-4064-bb80-cbd9f5523d7d
# ╠═5e0ac810-1689-4bfb-88a8-85504cd821b6
# ╠═1015f629-9818-4dbb-b574-97c552e96164
# ╠═62d41be1-2970-48e1-b689-c2ecf1ab10ce
# ╠═7670ffbc-3621-4acb-b1e3-b56668470c24
# ╠═d08cce9d-dd9c-4530-82ae-bf1db8be3281
# ╟─6a20715a-9dc9-461a-b6a2-f02522e277f0
# ╠═e63a5860-587c-4aff-91be-a894b3443943
# ╠═c866c71e-5dce-4af1-84bb-045815d1acb9
# ╠═75560584-1631-4c93-8808-e269d345e1c8
# ╠═456021aa-3dbf-4102-932f-43109905f9eb
# ╠═56e0da0b-8858-459a-9aae-0eba4797cc06
# ╠═cf83549c-03b1-4f7d-be29-d3a29da3e8f7
# ╠═85dfe3ed-cf70-4fde-a94a-f70c3fe4abf6
# ╠═2b6d959f-63bb-4d42-a34d-70d93f6a1636
# ╠═3a5e3aa0-e34e-457c-a405-9bc95cd88910
# ╠═5ec611f4-d3ea-4c5d-9221-4dcf2da3e18c
# ╠═3df40867-674b-405a-a05c-0733a4caaa67


text:
function IndicesDiagonal(v::AbstractVector)
    return IndicesArray(
        Diagonal(v),
        (axes(v, 1), axes(v, 1)),
        AllUnique,
        LengthChecked
       )
end

function LinearAlgebra.Diagonal(m::IndicesMatrix)
    return IndicesArray(
        diag(parent(m)),
        diagonal_indices(m),
        AllUnique,
        LengthChecked
       )
end

diagonal_indices(a::IndicesMatrix) = _diagonal_indices(a, findmin(size(a)))

function _diagional_indices(a::IndicesMatrix, imin::NTuple{2,Int})
    if last(imin) === 1
        return (axes(a, 1), unsafe_reindex(axes(a, 2), 1:first(imin)))
    else
        return (unsafe_reindex(axes(a, 1), 1:first(imin)), axes(a, 2))
    end
end

LinearAlgebra.diagm(v::IndicesVector) = IndicesArray(diagm(parent(v)), (axes(v, 1), axes(v, 1)))

ext: jl


text:
module DataTypes

using StaticArrays
using JuLIP: vecs, mat, JVec, energy, forces, virial
import IPFitting: vec_obs, devec_obs, eval_obs, weighthook, Dat, err_weighthook

export ENERGY, FORCES, VIRIAL

# some convenience functions to dispatch string arguments
vec_obs(s::AbstractString, args...) = vec_obs(Val(Symbol(s)), args...)
devec_obs(s::AbstractString, args...) = devec_obs(Val(Symbol(s)), args...)
eval_obs(s::AbstractString, args...) = eval_obs(Val(Symbol(s)), args...)
weighthook(s::AbstractString, args...) = weighthook(Val(Symbol(s)), args...)

# ------------------- ENERGY ------------------

const ENERGY = "E"
const ValE = Val{:E}
vec_obs(::ValE, E::Real) = [E]
devec_obs(::ValE, x::AbstractVector) = ((@assert length(x) == 1); x[1])
eval_obs(::ValE, B, at) = energy(B, at)
weighthook(::ValE, d::Dat) = 1.0 / length(d.at)
err_weighthook(::ValE, d::Dat) = 1.0 / length(d.at)

# ------------------- FORCES ------------------

const FORCES = "F"
const ValF = Val{:F}
vec_obs(v::ValF, F::AbstractVector{<:JVec}) = vec_obs(v, mat(F))
vec_obs(::ValF, F::AbstractMatrix) = vec(F)
devec_obs(::ValF, x::AbstractVector) = vecs(reshape(x, 3, :))
eval_obs(::ValF, B, at) = forces(B, at)

# ------------------- VIRIAL ------------------
# using Voigt convention for vectorising  symmetric 3 x 3 matrix
#  1  6  5
#  6  2  4
#  5  4  3
const _IV = [1,5,9,6,3,2]
const _IVst = SVector(1,5,9,6,3,2)
const VIRIAL = "V"
const ValV = Val{:V}
vec_obs(::ValV, v::AbstractVector) = (@assert length(v) == 6; collect(v))
vec_obs(::ValV, V::AbstractMatrix) = (@assert size(V) == (3,3); V[_IV])
devec_obs(::ValV, x::AbstractVector) =
   SMatrix{3,3}(x[1], x[6], x[5], x[6], x[2], x[4], x[5], x[4], x[3])
eval_obs(::ValV, B, at) = virial(B, at)
weighthook(::ValV, d::Dat) = 1.0 / length(d.at)
err_weighthook(::ValV, d::Dat) = 1.0 / length(d.at)

end


text:
## @lc code=start
using LeetCode
function minimum_time_2167(s::String)::Int
    right_cost = n = length(s)
    cost = s[1] == '0' ? 0 : 1
    for i in 2:n
        if s[i] == '1'
            right_cost = min(right_cost, cost + n - i + 1)
            cost = min(right_cost, i, cost + 2)
        end
    end
    cost
end
## @lc code=end

## @lc test=start
@testset "2167.minimum-time-to-remove-all-cars-containing-illegal-goods.jl" begin
    @test minimum_time_2167("1100101") == 5
    @test minimum_time_2167("0010") == 2
    @test minimum_time_2167("010110") == 5
    @test minimum_time_2167("0111001110") == 8
end
## @lc test=end

## @lc info=start
# link: solution to
## @lc info=end


text:
import Base: RefValue, @pure, show, clamp, findall
import Base: cumsum, cumprod, abs2
import LinearAlgebra: cholesky

export constant, get_last_error, err_to_string, sort_index, fir, iir
export mean_weighted, var_weighted, set_array_indexer, set_seq_param_indexer
export afeval, iota, sortbykey, select, device_mem_info, setafgcthreshold

const af_threshold = Ref(4*1024*1024*1024)

function device_gc()
    _error(ccall((:af_device_gc,af_lib),af_err,()), false)
end

function device_mem_info()
    alloc_bytes = RefValue{Csize_t}(0)
    alloc_buffers = RefValue{Csize_t}(0)
    lock_bytes = RefValue{Csize_t}(0)
    lock_buffers = RefValue{Csize_t}(0)
    _error(ccall((:af_device_mem_info,af_lib),af_err,
                 (Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t}),
                 alloc_bytes,alloc_buffers,lock_bytes,lock_buffers),
           false)
    (alloc_bytes[],alloc_buffers[],lock_bytes[],lock_buffers[])
end

function setafgcthreshold(threshold)
    af_threshold[] = threshold
end

function _afgc()
    for full in (false, true)
        alloc_bytes, alloc_buffers, lock_bytes, lock_buffers =  device_mem_info()
        if max(lock_bytes, alloc_bytes - lock_bytes) > af_threshold[]
            GC.gc(full)
            if full
                device_gc()
            end
        end
    end
    nothing
end

function release_array(arr::AFArray)
    _error(ccall((:af_release_array,af_lib),af_err,(af_array,),arr.arr), false)
end

toa(a) = issparse(a) ?  SparseMatrixCSC(a) : Array(a)
show(c::IOContext, a::AFArray) = (print(c, "AFArray: "); show(c, toa(a)))
show(io::IO, a::AFArray) = print(io, "AFArray: ", toa(a))
function show(io::IO, m::MIME"text/plain", a::AFArray)
    print(io, "AFArray: "); show(io, m, toa(a))
end

global const af_lib = Sys.isunix() ? "libaf" : "af"
global const bcast = Ref{Bool}(false)

function __init__()
    Libdl.dlopen(af_lib)

    backend_envvar_name = "JULIA_ARRAYFIRE_BACKEND"
    if haskey(ENV, backend_envvar_name)
        backend_str = lowercase(ENV[backend_envvar_name])
        backends = Dict(
            "cpu" => AF_BACKEND_CPU,
            "cuda" => AF_BACKEND_CUDA,
            "opencl" => AF_BACKEND_OPENCL
        )
        if haskey(backends, backend_str)
            set_backend(backends[backend_str])
        else
            error("Unknown arrayfire backend \"$backend_str\".")
        end
    end

    afinit()
    if !haskey(ENV, "JULIA_ARRAYFIRE_SILENT")
        afinfo()
    end
    set_seed(rand(RandomDevice(), UInt64))
    nothing
end

function _error(err::af_err, gc = true)
    if err == 0
        if gc
            _afgc()
        end
    elseif err == AF_ERR_NO_MEM
        error("GPU is out of memory, to avoid this in the future you can:
  setafgcthreshold(threshold) # lower garbage collect threshold (default 4Gb)
  finalize(array)             # manually free GPU memory")
    else
        str = err_to_string(err)
        str2 = get_last_error()
        error("ArrayFire Error ($err) : $str\n$str2")
    end
end

@pure batched(n1, n2) = max(n1, n2)

function typed(::Type{T1},::Type{T2}) where {T1,T2}
    if T1 == T2
        return T1
    elseif T1 == Complex{Float64} || T2 == Complex{Float64}
        return Complex{Float64}
    elseif T1 == Complex{Float32} || T2 == Complex{Float32}
        (T1 == Float64 || T2 == Float64) && return Complex{Float64}
        return Complex{Float32}
    elseif T1 == Float64 || T2 == Float64
        return Float64
    elseif T1 == Float32 || T2 == Float32
        return Float32
    elseif T1 == Float16 || T2 == Float16
        return Float16
    elseif T1 == UInt64 || T2 == UInt64
        return UInt64
    elseif T1 == Int64 || T2 == Int64
        return Int64
    elseif T1 == UInt32 || T2 == UInt32
        return UInt32
    elseif T1 == Int32 || T2 == Int32
        return Int32
    elseif T1 == UInt16 || T2 == UInt16
        return UInt16
    elseif T1 == Int16 || T2 == Int16
        return Int16
    elseif T1 == UInt8 || T2 == UInt8
        return UInt8
    elseif T1 == Bool || T2 == Bool
        return Bool
    else
        return Float32
    end
end

af_type(::Type{Float32})          = f32
af_type(::Type{Complex{Float32}}) = c32
af_type(::Type{Float64})          = f64
af_type(::Type{Complex{Float64}}) = c64
af_type(::Type{Bool})             = b8
af_type(::Type{Int32})            = s32
af_type(::Type{Int16})            = s16
af_type(::Type{UInt16})           = u16
af_type(::Type{UInt32})           = u32
af_type(::Type{UInt8})            = u8
af_type(::Type{Int64})            = s64
af_type(::Type{UInt64})           = u64
af_type(::Type{Float16})          = f16

function af_jltype(i::af_dtype)::Type
    if i == f32
        return Float32
    elseif i == c32
        return Complex{Float32}
    elseif i == f64
        return Float64
    elseif i == c64
        return Complex{Float64}
    elseif i == b8
        return Bool
    elseif i == s32
        return Int32
    elseif i == u32
        return UInt32
    elseif i == u8
        return UInt8
    elseif i == s64
        return Int64
    elseif i == u64
        return UInt64
    elseif i == s16
        return Int16
    elseif i == u16
        return UInt16
    elseif i == f16
        return Float16
    else
        error("Unknown type: $i")
    end
end

function get_numdims(arr::af_array)
    result = RefValue{UInt32}(0)
    _error(ccall((:af_get_numdims,af_lib),af_err,
                 (Ptr{UInt32},af_array),result,arr))
    Int(result[])
end

function get_type(arr::af_array)
    _type = RefValue{af_dtype}(0)
    _error(ccall((:af_get_type,af_lib),af_err,
                 (Ptr{af_dtype},af_array),_type,arr), false)
    af_jltype(_type[])
end

function check_type_numdims(arr::AFArray{T,N}) where {T,N}
    @assert get_type(arr) == T "type mismatch: $(get_type(arr)) != $T"
    @assert get_numdims(arr) == N "dims mismatch: $(get_numdims(arr)) != $N"
end

function convert_array(data::AbstractArray{T,N}) where {T,N}
    arr = RefValue{af_array}(0)
    sz = size(data)
    _error(ccall((:af_create_array,af_lib),af_err,
                 (Ptr{af_array},Ptr{Cvoid},UInt32,Ptr{dim_t},af_dtype),
                 arr,data,UInt32(length(sz)),[sz...],af_type(T)))
    AFArray{T,N}(arr[])
end

function convert_array(a::AFArray{T,N}) where {T,N}
    if issparse(a)
        a = full(a)
    end
    ret = Array{T,N}(undef, size(a))
    get_data_ptr(ret, a)
    ret
end

function convert_array_to_sparse(a::SparseMatrixCSC)
    sz = size(a)
    at = sparse(transpose(a))
    colptr = AFArray(Vector{Int32}(at.colptr.-1))
    rowval = AFArray(Vector{Int32}(at.rowval.-1))
    create_sparse_array(sz[1], sz[2], AFArray(at.nzval), colptr, rowval, AF_STORAGE_CSR)
end

function convert_array_to_sparse(a::AFArray)
    @assert issparse(a) "AFArray is not sparse"
    sz = size(a)
    @assert length(sz) == 2 "AFArray is not a matrix"
    nzval, colptr, rowval, d = sparse_get_info(a)
    if d == AF_STORAGE_CSR
        transpose(SparseMatrixCSC(sz[2], sz[1],
                                  Array(colptr.+1), Array(rowval.+1), Array(nzval)))
    else
        convert_array_to_sparse(sparse_convert_to(a, AF_STORAGE_CSR))
    end
end

function recast_array(::Type{AFArray{T1}},_in::AFArray{T2,N}) where {T1,N,T2}
    out = RefValue{af_array}(0)
    _error(ccall((:af_cast,af_lib),af_err,
                 (Ptr{af_array},af_array,af_dtype),out,_in.arr,af_type(T1)))
    AFArray{T1,N}(out[])
end

AFArray!(arr::af_array) = AFArray{get_type(arr), get_numdims(arr)}(arr)

function constant(val::T,sz::NTuple{N,Int}) where {T<:Real,N}
    arr = RefValue{af_array}(0)
    _error(ccall((:af_constant,af_lib),af_err,
                 (Ptr{af_array},Cdouble,UInt32,Ptr{dim_t},af_dtype),
                 arr,Cdouble(val),UInt32(N),[sz...],af_type(T)))
    AFArray{T,N}(arr[])
end

function constant(val::Complex{Bool},sz::NTuple{N,Int}) where {N}
    arr = RefValue{af_array}(0)
    _error(ccall((:af_constant_complex,af_lib),af_err,
                 (Ptr{af_array},Cdouble,Cdouble,UInt32,Ptr{dim_t},af_dtype),
                 arr,Cdouble(real(val)),Cdouble(imag(val)),UInt32(N),[sz...],c32))
    AFArray{Complex{Float32},N}(arr[])
end

function constant(val::T,sz::NTuple{N,Int}) where {T<:Complex,N}
    arr = RefValue{af_array}(0)
    _error(ccall((:af_constant_complex,af_lib),af_err,
                 (Ptr{af_array},Cdouble,Cdouble,UInt32,Ptr{dim_t},af_dtype),
                 arr,Cdouble(real(val)),Cdouble(imag(val)),UInt32(N),[sz...],af_type(T)))
    AFArray{T,N}(arr[])
end

function constant(val::Int,sz::NTuple{N,Int}) where {N}
    arr = RefValue{af_array}(0)
    _error(ccall((:af_constant_long,af_lib),af_err,
                 (Ptr{af_array},intl,UInt32,Ptr{dim_t}),
                 arr,val,UInt32(N),[sz...]))
    AFArray{Int,N}(arr[])
end

function constant(val::UInt,sz::NTuple{N,Int}) where {N}
    arr = RefValue{af_array}(0)
    _error(ccall((:af_constant_ulong,af_lib),af_err,
                 (Ptr{af_array},uintl,UInt32,Ptr{dim_t}),
                 arr,val,UInt32(N),[sz...]))
    AFArray{UInt,N}(arr[])
end

function select(cond::AFArray{Bool},a::AFArray{T1,N1},b::AFArray{T2,N2}) where {T1,N1,T2,N2}
    out = RefValue{af_array}(0)
    _error(ccall((:af_select,af_lib),af_err,
                 (Ptr{af_array},af_array,af_array,af_array),
                 out,cond.arr,a.arr,b.arr))
    AFArray{typed(T1,T2),batched(N1,N2)}(out[])
end

function select(cond::AFArray{Bool},a::AFArray{T1,N1},b::T2) where {T1,N1,T2<:Real}
    out = RefValue{af_array}(0)
    _error(ccall((:af_select_scalar_r,af_lib),af_err,
                 (Ptr{af_array},af_array,af_array,Cdouble),
                 out,cond.arr,a.arr,Cdouble(b)))
    AFArray{typed(T1,T2),N1}(out[])
end

function select(cond::AFArray{Bool},a::T1,b::AFArray{T2,N2}) where {T1,T2,N2}
    out = RefValue{af_array}(0)
    _error(ccall((:af_select_scalar_l,af_lib),af_err,
                 (Ptr{af_array},af_array,Cdouble,af_array),
                 out,cond.arr,Cdouble(a),b.arr))
    AFArray{typed(T1,T2),N2}(out[])
end

function err_to_string(err::af_err)
    unsafe_string(ccall((:af_err_to_string,af_lib),Cstring,(af_err,),err))
end

function get_last_error()
    msg = RefValue{Cstring}()
    len = RefValue{dim_t}(0)
    ccall((:af_get_last_error,af_lib),Cvoid,(Ptr{Cstring},Ptr{dim_t}),msg,len)
    unsafe_string(msg[])
end

function cat(dim::Integer,first::AFArray{T,N1},second::AFArray{T,N2}) where {T,N1,N2}
    out = RefValue{af_array}(0)
    _error(ccall((:af_join,af_lib),af_err,
                 (Ptr{af_array},Cint,af_array,af_array),
                 out,Cint(dim - 1),first.arr,second.arr))
    AFArray!(out[])
end

function cat(A::AFArray...; dims::Integer)
    out = RefValue{af_array}(0)
    _error(ccall((:af_join_many,af_lib),af_err,
                 (Ptr{af_array},Cint, Cuint, Ptr{af_array}),
                 out,Cint(dims - 1), Cuint(length(A)), [map(x->x.arr, A)...]))
    AFArray!(out[])
end

hcat(first::AFArray, second::AFArray) = cat(2, first, second)
vcat(first::AFArray, second::AFArray) = cat(1, first, second)

import DSP: conv

function conv(signal::AFArray{T,N}, filter::AFArray) where {T,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_convolve1,af_lib),af_err,
                 (Ptr{af_array},af_array,af_array,af_conv_mode,af_conv_domain),
                 out,signal.arr,filter.arr,AF_CONV_EXPAND,AF_CONV_AUTO))
    AFArray{T,N}(out[])
end

function conv_fft(signal::AFArray{T,N}, filter::AFArray) where {T,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft_convolve1,af_lib),af_err,(Ptr{af_array},af_array,af_array,af_conv_mode),
                 out,signal.arr,filter.arr,AF_CONV_EXPAND))
    AFArray{T,N}(out[])
end

norm(a::AFArray{T}) where T = T(norm(a, AF_NORM_EUCLID, 0, 0))
vecnorm(a::AFArray{T}) where T = T(norm(a, AF_NORM_EUCLID, 0, 0))

function svd(_in::AFArray{T,2}) where {T}
    u = RefValue{af_array}(0)
    s = RefValue{af_array}(0)
    vt = RefValue{af_array}(0)
    _error(ccall((:af_svd,af_lib),af_err,
                 (Ptr{af_array},Ptr{af_array},Ptr{af_array},af_array),
                 u,s,vt,_in.arr))
    (AFArray{T,2}(u[]),AFArray{T,1}(s[]),AFArray{T,2}(vt[]))
end

function sort(_in::AFArray{T,N},dim::Integer,isAscending::Bool=true) where {T,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_sort,af_lib),af_err,
                 (Ptr{af_array},af_array,UInt32,Bool),
                 out,_in.arr,UInt32(dim - 1),isAscending))
    AFArray{T,N}(out[])
end

function sort_index(_in::AFArray{T,N},dim::Integer=1,isAscending::Bool=true) where {T,N}
    out = RefValue{af_array}(0)
    indices = RefValue{af_array}(0)
    _error(ccall((:af_sort_index,af_lib),af_err,
                 (Ptr{af_array},Ptr{af_array},af_array,UInt32,Bool),
                 out,indices,_in.arr,UInt32(dim - 1),isAscending))
    (AFArray{T,N}(out[]),AFArray{UInt32,N}(indices[])+UInt32(1))
end

function sortperm(a::AFArray{T,N}, dim::Integer=1,isAscending::Bool=true) where {T,N}
    sort_index(a,dim,isAscending)[2]
end

function mean(_in::AFArray{T,N},dim::dim_t) where {T,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_mean,af_lib),af_err,(Ptr{af_array},af_array,dim_t),
                 out,_in.arr,dim-1))
    AFArray{T,N}(out[])
end

function mean_weighted(_in::AFArray{T,N},weights::AFArray,dim::dim_t) where {T,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_mean_weighted,af_lib),af_err,
                 (Ptr{af_array},af_array,af_array,dim_t),
                 out,_in.arr,weights.arr,dim-1))
    AFArray{T,N}(out[])
end

function var(_in::AFArray{T,N},isbiased::Bool,dim::dim_t) where {T,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_var,af_lib),af_err,
                 (Ptr{af_array},af_array,Bool,dim_t),
                 out,_in.arr,isbiased,dim-1))
    AFArray{T,N}(out[])
end

function var_weighted(_in::AFArray{T,N},weights::AFArray,dim::dim_t) where {T,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_var_weighted,af_lib),af_err,
                 (Ptr{af_array},af_array,af_array,dim_t),
                 out,_in.arr,weights.arr,dim-1))
    AFArray{T,N}(out[])
end

function std(_in::AFArray{T,N},dim::dim_t) where {T,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_stdev,af_lib),af_err,
                 (Ptr{af_array},af_array,dim_t),out,_in.arr,dim-1))
    AFArray{T,N}(out[])
end

function median(_in::AFArray{T,N},dim::dim_t) where {T,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_median,af_lib),af_err,
                 (Ptr{af_array},af_array,dim_t),
                 out,_in.arr,dim-1))
    AFArray{T,N}(out[])
end

function set_seq_param_indexer(_begin::Real,_end::Real,step::Real,dim::dim_t,is_batch::Bool)
    indexer = RefValue{af_index_t}(0)
    _error(ccall((:af_set_seq_param_indexer,af_lib),af_err,
                 (Ptr{af_index_t},Cdouble,Cdouble,Cdouble,dim_t,Bool),
                 indexer,Cdouble(_begin),Cdouble(_end),Cdouble(step),dim-1,is_batch))
    indexer[]
end

function clamp(_in::AFArray{T1,N1},lo::AFArray{T2,N2},hi::AFArray{T2,N2}) where {T1,N1,T2,N2}
    out = RefValue{af_array}(0)
    _error(ccall((:af_clamp,af_lib),af_err,
                 (Ptr{af_array},af_array,af_array,af_array,Bool),
                 out,_in.arr,lo.arr,hi.arr,bcast[]))
    AFArray{typed(T1,T2),batched(N1,N2)}(out[])
end

function af_where(_in::AFArray{T,N}) where {T,N}
    idx = RefValue{af_array}(0)
    _error(ccall((:af_where,af_lib),af_err,(Ptr{af_array},af_array),idx,_in.arr))
    return AFArray{UInt32,1}(idx[])
end

function findall(_in::AFArray{T,N}) where {T,N}
    out = af_where(_in)
    if length(out) > 0
        out = out + UInt32(1)
    end
    return out
end

cumsum(a::AFArray, dim::Int) = scan(a, dim, AF_BINARY_ADD, true)
cumprod(a::AFArray, dim::Int) = scan(a, dim, AF_BINARY_MUL, true)
cummin(a::AFArray, dim::Int) = scan(a, dim, AF_BINARY_MIN, true)
cummax(a::AFArray, dim::Int) = scan(a, dim, AF_BINARY_MAX, true)

function sync(a::AFArray)
    afeval(a)
    sync(get_device_id(a))
    a
end

function afeval(a::AFArray)
    _error(ccall((:af_eval,af_lib),af_err,(af_array,),a.arr))
    a
end

function afeval(A::AFArray...)
    _error(ccall((:af_eval_multiple,af_lib),af_err,(Cuint, Ptr{af_array}),
                 Cuint(length(A)), [map(x->x.arr, A)...]))
    A
end
afeval(A::Tuple) = afeval(A...)

function cholesky(_in::AFArray{T,N},is_upper::Bool=false) where {T,N}
    out = RefValue{af_array}(0)
    info = RefValue{Cint}(0)
    _error(ccall((:af_cholesky,af_lib),af_err,(Ptr{af_array},Ptr{Cint},af_array,Bool),out,info,_in.arr,is_upper))
    (AFArray{T,N}(out[]),info[])
end

abs2(a::AFArray{T}) where {T<:Real} = a.*a
abs2(a::AFArray{T}) where {T<:Complex} = (r = real(a); i = imag(a); r.*r+i.*i)

function complex(lhs::AFArray{T1,N1},rhs::AFArray{T2,N2}) where {T1,N1,T2,N2}
    batch = bcast[]
    out = RefValue{af_array}(0)
    _error(ccall((:af_cplx2,af_lib),af_err,(Ptr{af_array},af_array,af_array,Bool),out,lhs.arr,rhs.arr,batch))
    AFArray{Complex{typed(T1,T2)},batched(N1,N2)}(out[])
end

function fir(b::AFArray,x::AFArray{T,N}) where {T,N}
    y = RefValue{af_array}(0)
    _error(ccall((:af_fir,af_lib),af_err,(Ptr{af_array},af_array,af_array),y,b.arr,x.arr))
    AFArray{T,N}(y[])
end

function iir(b::AFArray{T},a::AFArray{T},x::AFArray{T,N}) where {T,N}
    y = RefValue{af_array}(0)
    _error(ccall((:af_iir,af_lib),af_err,(Ptr{af_array},af_array,af_array,af_array),y,b.arr,a.arr,x.arr))
    AFArray{T,N}(y[])
end

function iota(dims::NTuple{N,Int}, typ::Type{T} = Int32) where {T,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_iota,af_lib), af_err,
                 (Ptr{af_array},UInt32,Ptr{dim_t},UInt32,Ptr{dim_t},af_dtype),
                 out,UInt32(N),[dims...],UInt32(1),[1],af_type(T)))
    AFArray{T,N}(out[])+T(1)
end

function sortbykey(keys::AFArray{T1,N},values::AFArray{T2,N},dim::Integer,isAscending::Bool=true) where {T1,T2,N}
    out_keys = RefValue{af_array}(0)
    out_values = RefValue{af_array}(0)
    _error(ccall((:af_sort_by_key,af_lib),af_err,(Ptr{af_array},Ptr{af_array},af_array,af_array,UInt32,Bool),out_keys,out_values,keys.arr,values.arr,UInt32(dim - 1),isAscending))
    (AFArray{T1,N}(out_keys[]),AFArray{T2,N}(out_values[]))
end

function sum(_in::AFArray{UInt8,N},dim::Integer) where N
    out = RefValue{af_array}(0)
    _error(ccall((:af_sum,af_lib),af_err,(Ptr{af_array},af_array,Cint),out,_in.arr,Cint(dim - 1)))
    AFArray{UInt32,N}(out[])
end

function sum(_in::AFArray{Bool,N},dim::Integer) where N
    out = RefValue{af_array}(0)
    _error(ccall((:af_sum,af_lib),af_err,(Ptr{af_array},af_array,Cint),out,_in.arr,Cint(dim - 1)))
    AFArray{UInt32,N}(out[])
end


text:
@testset "System Test Core Creation" begin
    brk = h.helicsCreateBroker("inproc", "gbrokerc", "--root")

    argv = ["", "--name=gcore", "--timeout=2000", "--broker=gbrokerc"]

    cr = h.helicsCreateCoreFromArgs("inproc", "", argv)
    @test h.helicsCoreGetIdentifier(cr) == "gcore"

    argv = ["--name=gcore2", "--log-level=what_logs?"]

    @test_throws h.HELICSErrorInvalidArgument cr2 = h.helicsCreateCoreFromArgs("inproc", "", argv)

    h.helicsBrokerDisconnect(brk)
    h.helicsCoreDisconnect(cr)

    @test h.helicsBrokerIsConnected(brk) == false
end

@testset "System Test Broker Creation" begin

    argv = ["", "--name=gbrokerc", "--timeout=2000", "--root"]

    brk = h.helicsCreateBrokerFromArgs("inproc", "", argv)
    @test h.helicsBrokerGetIdentifier(brk) == "gbrokerc"

    argv[2] = "--name=gbrokerc2"
    argv[3] = "--log-level=what_logs?"

    @test_throws h.HELICSErrorInvalidArgument brk2 = h.helicsCreateBrokerFromArgs("inproc", "", argv)

    h.helicsBrokerDisconnect(brk)

end

@testset "System Test Broker Global Value" begin

    brk = h.helicsCreateBroker("inproc", "gbroker", "--root")
    globalVal = "this is a string constant that functions as a global"
    globalVal2 = "this is a second string constant that functions as a global"
    h.helicsBrokerSetGlobal(brk, "testglobal", globalVal)
    q = h.helicsCreateQuery("global", "testglobal")
    res = h.helicsQueryBrokerExecute(q, brk)
    @test res == globalVal

    h.helicsBrokerSetGlobal(brk, "testglobal2", globalVal2)
    h.helicsQueryFree(q)
    q = h.helicsCreateQuery("global", "testglobal2")
    res = h.helicsQueryBrokerExecute(q, brk)
    @test res == globalVal2

    h.helicsBrokerDisconnect(brk)
    h.helicsQueryFree(q)
    @test h.helicsBrokerIsConnected(brk) == false

end

@testset "System Test Core Global Value" begin

    brk = h.helicsCreateBroker("zmq", "gbrokerc", "--root")
    cr = h.helicsCreateCore("zmq", "gcore", "--broker=gbrokerc")

    globalVal = "this is a string constant that functions as a global"
    globalVal2 = "this is a second string constant that functions as a global"

    h.helicsCoreSetGlobal(cr, "testglobal", globalVal)

    # q = h.helicsCreateQuery("global", "testglobal")
    # TODO: This hangs on core execute
    # res = h.helicsQueryCoreExecute(q, cr)
    # @test res == globalVal
    # h.helicsQueryFree(q)
    @test_broken false

    h.helicsCoreDisconnect(cr)
    h.helicsBrokerDisconnect(brk)

    @test h.helicsBrokerIsConnected(brk) == false

end
@testset "System Test Federate Global Value" begin

    brk = h.helicsCreateBroker("inproc", "gbrokerc", "--root")
    cr = h.helicsCreateCore("inproc", "gcore", "--broker=gbrokerc")

    # test creation of federateInfo from command line arguments
    argv = [
        ""
        "--corename=gcore"
        "--type=inproc"
        "--period=1.0"
    ]

    fi = h.helicsCreateFederateInfo()
    h.helicsFederateInfoLoadFromArgs(fi, argv)

    fed = h.helicsCreateValueFederate("fed0", fi)

    argv[4] = "--period=frogs" #this is meant to generate an error in command line processing

    fi2 = h.helicsFederateInfoClone(fi)
    @test_throws h.HELICSErrorInvalidArgument h.helicsFederateInfoLoadFromArgs(fi2, argv)

    h.helicsFederateInfoFree(fi2)
    h.helicsFederateInfoFree(fi)

    globalVal = "this is a string constant that functions as a global"
    globalVal2 = "this is a second string constant that functions as a global"
    h.helicsFederateSetGlobal(fed, "testglobal", globalVal)
    q = h.helicsCreateQuery("global", "testglobal")
    res = h.helicsQueryExecute(q, fed)
    @test res == globalVal
    h.helicsFederateSetGlobal(fed, "testglobal2", globalVal2)
    h.helicsQueryFree(q)
    q = h.helicsCreateQuery("global", "testglobal2")
    h.helicsQueryExecuteAsync(q, fed)
    while (h.helicsQueryIsCompleted(q) == false)
        sleep(0.20)
    end
    res = h.helicsQueryExecuteComplete(q)
    @test res == globalVal2

    q2 = h.helicsCreateQuery("", "isinit")
    h.helicsQueryExecuteAsync(q2, fed)
    while (h.helicsQueryIsCompleted(q2) == false)
        sleep(0.20)
    end
    res = h.helicsQueryExecuteComplete(q2)
    @test res == "false"

    h.helicsFederateFinalize(fed)

    h.helicsCoreDisconnect(cr)
    h.helicsBrokerDisconnect(brk)

    h.helicsQueryFree(q)
    h.helicsQueryFree(q2)
    @test h.helicsBrokerIsConnected(brk) == false

    h.helicsBrokerDisconnect(brk)
    h.helicsCoreDisconnect(cr)

    @test h.helicsBrokerIsConnected(brk) == false

end

@testset "System Tests Core Logging" begin
    lfile = "log.txt";
    rm(lfile, force=true)
    core = h.helicsCreateCore("inproc", "clog", "--autobroker --log_level=trace")
    h.helicsCoreSetLogFile(core, lfile)
    h.helicsCoreDisconnect(core)
    h.helicsCloseLibrary()
    @test isfile(lfile)
    rm(lfile, force=true)
end

@testset "System Tests Broker Logging" begin
    lfile = "log.txt";
    rm(lfile, force=true)
    broker = h.helicsCreateBroker("inproc", "blog", "--log_level=trace")
    h.helicsBrokerSetLogFile(broker, lfile)
    h.helicsBrokerDisconnect(broker)
    h.helicsCloseLibrary()
    @test isfile(lfile)
    rm(lfile, force=true)
end

@testset "System Tests Federate Logging" begin

    lfile = "log.txt"
    rm(lfile, force=true)
    core = h.helicsCreateCore("inproc", "clogf", "--autobroker --log_level=trace")

    fi = h.helicsCreateFederateInfo()
    h.helicsFederateInfoSetBrokerKey(fi, "key")
    h.helicsFederateInfoSetCoreName(fi, "clogf")
    fed = h.helicsCreateValueFederate("f1", fi)
    h.helicsFederateSetLogFile(fed, lfile)
    h.helicsFederateLogLevelMessage(fed, 7, "hello")
    h.helicsFederateLogErrorMessage(fed, "hello")
    h.helicsFederateLogDebugMessage(fed, "hello")
    h.helicsFederateLogWarningMessage(fed, "hello")
    h.helicsFederateClearMessages(fed)
    h.helicsCoreSetLogFile(core, lfile)
    h.helicsCoreDisconnect(core)
    h.helicsFederateFinalize(fed)
    h.helicsFederateInfoFree(fi)
    h.helicsCloseLibrary()

    @test isfile(lfile)
    rm(lfile, force=true)

end


text:
module ParticleDetection
    using KernelOps

    include("radius_find.jl")
    export find_best_radius

    @inline function op_mean(A,Is,I)
        m=zero(eltype(A))
        @simd for i in Is
            @inbounds m+=A[i]
        end
        m/length(Is)
    end
    @inline function rep_tuple(val,A::AbstractArray{T,N}) where {T,N}
        ntuple(N) do x
            val
        end
    end
    struct op_gaussian{T,N} <: Function
        sn::Int
        w::Array{T,N}
        function op_gaussian(A,sn)
            w=similar(eltype(A)[],rep_tuple(sn*2+1,A))
            for I in CartesianIndices(w)
                deltaX=Tuple(I).-sn.-1
                delta=sum(pow2,deltaX)/(sn*sn)
                w[I]=exp(-2*delta)
            end
            new{eltype(A),ndims(A)}(sn,w)
        end
    end

    @inline function (p::op_gaussian)(A::AbstractArray{T}, Is, I) where T
        ws=zero(T)
        s=ws
        rc=CartesianIndex(rep_tuple(p.sn+1,A))
        @inbounds for i in Is
            ni=i-I+rc
            s+=A[i]*p.w[ni]
            ws+=p.w[ni]
        end
        s/ws
    end

    @inline function op_max(A,Is,I)
        @inbounds m=A[I]
        @simd for i in Is
            @inbounds m=max(m,A[i])
        end
        m
    end

    @inline function pow2(x)
        x*x
    end

    @inline function bp_filter(image,sn,so)
        T=eltype(image)
        G=op_gaussian(image,sn)
        Z=zero(T)
        fobject=KernelOp(op_mean,image,rep_tuple(so,image))
        fnoise=KernelOp(G,image,rep_tuple(sn,image))
        pre=similar(T[],size(image))
        Threads.@threads for i in CartesianIndices(pre)
            @inbounds pre[i]=max(fnoise[i]-fobject[i],Z)
        end
        pre
    end

    @inline function local_maxima(image,abs_th=zero(eltype(image)))
        abs_th=convert(eltype(image),abs_th)
        max_im=KernelOp(op_max,image,rep_tuple(1,image))
        mask=falses(size(image)...)
        Threads.@threads for i in CartesianIndices(mask)
            @inbounds a=max_im[i]
            @inbounds b=image[i]
            @fastmath c=((a == b) && (b > abs_th))
            mask[i]=c
        end
        mask
    end
    export detect_particles
"""
detect_particles(image, noise_scale, object_scale, abs_threshold)
"""
    function detect_particles(image::AbstractArray{T},noise_scale::Int,object_scale::Int,abs_threshold::T) where T
        lm=local_maxima(bp_filter(image,noise_scale,object_scale),abs_threshold)
        [Tuple(p) for p in CartesianIndices(lm) if lm[p]]
    end
end # module


text:
using Test
using SparseArrays
using LinearAlgebra
import sbp
using sbp.CollocatedAcoustic: init_operators, pressure_norm, velocity_norm,
                             energy_norm, contravariant_metric_tensor,
                             divergence, gradient,
                             spatial_discretization, pressure_sat, grid,
                             pressure_source
using sbp.Sparse: diag_block_matrix_2x2
using sbp.Grid: grid_xp, grid_xm, grid_2d_x, grid_2d_y

st=sbp.Strand1994

nx=10
ny=20

function mapping(nx, ny)
        x1, h = grid_xp(nx)
        y1, h = grid_xp(ny)

        p = grid_2d_x(x1, ny)
        q = grid_2d_y(y1, nx)
        x = p
        y = (1 .+ 0.1 * sin.(p)).*q
        return x, y
end

fx, fy = mapping(nx, ny)
desc_D, desc_H = st.load_description(order=2)
ops = init_operators(nx, ny, st, desc_D, desc_H, fx, fy)

HJp = pressure_norm(ops)
HJ = velocity_norm(ops)
Gp = contravariant_metric_tensor(ops)
Ek = HJ*Gp
Ap = divergence(ops)
Av = gradient(ops, Gp)
S = pressure_sat(ops, Gp)
H = energy_norm(HJp, HJ, Gp)
A = spatial_discretization(Ap, Av - S)
d = pressure_source(ops, 0.5, 0.5, 4, 4)

@testset "Grids" begin
        x, y = grid(nx, ny)
        @test size(x,1) == nx * ny
        @test size(y,1) == nx * ny
end

@testset "Kinetic energy tensor" begin
        # Check symmetry
        @test isapprox(Ek, (Ek)')
        # Check positivity
        lam = real(eigvals(Matrix(Ek)))
        @test minimum(lam) > 0

end

@testset "Stability" begin
        # Check that the discretization is energy conservative to rounding error
        lam = real(eigvals(Matrix(A)))
        @test minimum(lam) > -1e13
        @test maximum(lam) < 1e13

        E = H*A + (H*A)'
        lam = real(eigvals(Matrix(E)))
        @test minimum(lam) > -1e13
        @test maximum(lam) < 1e13

end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
#
# This example shows the use of processes and channels with Simulate.jl
#
# pb, 2019-11-09
#

using Simulate, Printf, Random

function simple(c::Clock, input::Channel, output::Channel, name, id, op)
    token = take!(input)         # take something from the input
    # @printf("%5.2f: %s %d took token %d\n", tau(), name, id, token)
    now!(c, SF(println, @sprintf("%5.2f: %s %d took token %d", tau(c), name, id, token)))
    d = delay!(c, rand())           # after a delay
    put!(output, op(token, id))  # put it out with some op applied
end

clk = Clock()
Random.seed!(123)

ch1 = Channel(32)  # create two channels
ch2 = Channel(32)

for i in 1:2:8    # create and register 8 SimProcesses SP
    process!(clk, SP(i, simple, ch1, ch2, "foo", i, +))
    process!(clk, SP(i+1, simple, ch2, ch1, "bar", i+1, *))
end

put!(ch1, 1)  # put first token into channel 1
yield()
run!(clk, 10)   # and run for 10 time units


text:
#TODO use Base.@kwdef ?
mutable struct SaveLocation
    prefix::String
    suffix::Maybe{String}
    SaveLocation(prefix, suffix) = new(prefix, suffix)
end

SaveLocation(;prefix=path_prefix(), suffix=missing) = SaveLocation(prefix, suffix)
SaveLocation(prefix::AbstractString) = SaveLocation(prefix, "")
SaveLocation(sl::SaveLocation) = sl

"""
    abstract AbstractDataDescriptor

Supertype for all `DataDescriptor`s.
A `DataDescriptor` holds information about a set parameters that can be used to obtain the corresponding `Data`. Either by `load`ing or by `create`ing it.

## Necessary methods to implement
 - `_filename(descriptor)` - return something like "mydata/p1_p2_p3". The prefix/suffix/file ending will all be automatically added.
 - `create(descriptor)` - Create the corresponding `Data` type and compute it's data. If additional data is needed use `load_or_create` to get it.

## Optional methods
 - `load_mydata(param1, param2,...)` - construct a Descriptor with the necessary fields and `load` it
"""
abstract type AbstractDataDescriptor end

"""
    _filename(descriptor)
Generate the middle part of the save location. A prefix will by prepended by `joinpath` and a suffix with filetype will be appended.
"""
function _filename end

"""
    abstract AbstractData

Supertype for simulation data. Subtypes consist of an appropriate `DataDescriptor` and the actual data.
The `DataDescriptor` object should be stored in a field called `descriptor` or else they need to overwrite the function `descriptor(::MyDataType)` to return the appropriate `DataDescriptor`.
"""
abstract type AbstractData end

descriptor(data::AbstractData) = data.descriptor

#forward properties to descriptor
Base.getproperty(data::AbstractData, s::Symbol) = hasfield(typeof(data), s) ? getfield(data, s) : getproperty(data.descriptor, s)

"""
    AbstractSimpleData <: AbstractData

Supertype `Data` types that only consist of a single array and the descriptor.
"""
abstract type AbstractSimpleData <: AbstractData end

data(data::AbstractSimpleData) = data.data
Base.getindex(sdata::AbstractSimpleData, inds...) = getindex(data(sdata), inds...)
Base.setindex!(sdata::AbstractSimpleData, args...) = setindex!(data(sdata), args...)
Base.size(sdata::AbstractSimpleData) = size(data(sdata))
Base.size(sdata::AbstractSimpleData, dim) = size(data(sdata), dim)
Base.axes(sdata::AbstractSimpleData) = axes(data(sdata))
Base.axes(sdata::AbstractSimpleData, dim) = axes(data(sdata), dim)

"""
    datapath(dataOrDesc)
    datapath(dataOrDesc, path)
    datapath(dataOrDesc, prefix, suffix)
    datapath(dataOrDesc[, savelocation]; prefix, suffix)

Construct the path to where the data for this object should be saved to/loaded from.
There are a few different ways the path can be constructed.
 1. If a path is specified directly, take that
 2. Take prefix/suffix kwargs and fill in missing information from either the directly specified `SaveLocation` or take the one from the descriptor
"""
datapath(data::AbstractData, args...; kwargs...) = datapath(descriptor(data), args...; kwargs...) # unwrap Data -> DataDescriptor
# directly specified path takes precedence
datapath(::AbstractDataDescriptor, path::AbstractString) = path
# else merge SaveLocation and keyword args (latter have precedence)
datapath(desc::AbstractDataDescriptor, pathdata::SaveLocation=desc.pathdata; prefix=pathdata.prefix, suffix=pathdata.suffix) = datapath(desc, prefix, suffix)
# construct path
function datapath(desc::AbstractDataDescriptor, prefix::AbstractString, suffix::Maybe{<:AbstractString})
    if !ismissing(suffix) && length(suffix) > 0
        joinpath(prefix, "$(_filename(desc))-$(suffix).jld2")
    else
        joinpath(prefix, "$(_filename(desc)).jld2")
    end
end

"""
    create(descriptor)

Perform the calculation described by the `DataDescriptor`.
"""
function create end

"""
    save(data)
    save(data, path)
    save(data[, savelocation]]; prefix, suffix)

Save the `Data` to the specified location.
See [`datapath`](@ref) for documention on the different variants to specifiy a path.
"""
save(data::AbstractData, args...; kwargs...) = save(data, datapath(data, args...; kwargs...))

function save(data::AbstractData, path::AbstractString)
    dname = dirname(path)
    if !isdir(dname)
        logmsg("Save directory: $dname does not exist. Creating!")
        mkpath(dname)
    end
    logmsg("Saving file: $path")
    JLD2.jldsave(path; data, version="0.3")
end

"""
    load(descriptor)
    load(descriptor, path)
    load(descriptor[, savelocation]]; prefix, suffix)

Load the data described by the `DataDescriptor`.
See [`datapath`](@ref) for documention on the different variants to specifiy a path.
"""
function load(desc::AbstractDataDescriptor, args...; kwargs...)
    data = load(datapath(desc, args...; kwargs...))
    # copy savelocation data over in case it's a legacy file
    data.pathdata.prefix = desc.pathdata.prefix
    data.pathdata.suffix = desc.pathdata.suffix
    data
end

function load(path::AbstractString; throwerror=true)
    if !isfile(path)
        if !throwerror
            return nothing
        else
            error("$(abspath(path)) does not exist!")
        end
    end
    data = JLD2.load(path)
    if !haskey(data, "version")
        error("File too old. Try version 0.2 to load")
    elseif data["version"] != "0.3"
        error("File too old. Try version $(data["version"]) to load")
    end
    return data["data"]
end

"""
    load_or_create(descriptor; dosave)
    load_or_create(descriptor, path; dosave)
    load_or_create(descriptor[, savelocation]; dosave, prefix, suffix)

First try to `load` the data. If the file is not present or its contents do not match the requirements of the `descriptor`
then `create` the data.
If `save` is true, also `save` the data (possibly overwriting what was there). Default: true
"""
function load_or_create(desc::AbstractDataDescriptor, pathargs...; dosave=true, kwargs...)
    logmsg("load_or_create with spec $(desc)")
    p = datapath(desc, pathargs...; kwargs...)
    if !isfile(p)
        logmsg("No data found at $(p).")
    else
        logmsg("Found existing data for $(desc).")
        data = load(p)
        try
            compat = descriptor(data) == desc
            (compat || ismissing(compat)) && return data
        catch e;
            logmsg("Some error occured during loading: $(string(e))")
            display(stacktrace(catch_backtrace()))
        end
        # only if compat == false continue
        logmsg("Loaded data does not fit requirements.")
    end
    data = create(desc)
    dosave && save(data, p)
    data
end


text:
export SplitBregman

mutable struct SplitBregman{matT,vecT,opT,ropT,rvecT,preconT} <: AbstractLinearSolver
  # oerators and regularization
  A::matT
  reg::Vector{Regularization}
  regTrafo::Vector{ropT}
  y::vecT
  # fields and operators for x update
  op::opT
  β::vecT
  β_yj::vecT
  y_j::vecT
  # fields for primal & dual variables
  u::vecT
  v::Vector{vecT}
  vᵒˡᵈ::Vector{vecT}
  b::Vector{vecT}
  # other parameters
  precon::preconT
  ρ::rvecT
  iterations::Int64
  iterationsInner::Int64
  iterationsCG::Int64
  # state variables for CG
  cgStateVars::CGStateVariables
  # convergence parameters
  rk::rvecT
  sk::vecT
  eps_pri::rvecT
  eps_dt::vecT
  # eps_dual::Float64
  σᵃᵇˢ::Float64
  absTol::Float64
  relTol::Float64
  tolInner::Float64
  #counter for internal iterations
  iter_cnt::Int64
  normalizeReg::Bool
  regFac::Float64
end

"""
    SplitBregman(A::matT, x::vecT=zeros(eltype(A),size(A,2))
                ; reg=nothing, regName=["L1","TV"], λ=[0.0,0.0], kargs...) where {matT,vecT}

creates a `SplitBregman` object for the system matrix `A`.

# Arguments
* `A::matT`                     - system matrix
* `x::vecT`                     - Array with the same type and size as the solution
* (`reg=nothing`)               - Regularization object
* (`regName=["L1"]`)            - name of the regularizations to use (if reg==nothing)
* (`λ=[0.0]`)                   - Regularization paramters
* (`regTrafo=nothing`)          - transformations applied inside each regularizer
* (`precon=Identity()`)         - preconditionner for the internal CG algorithm
* (`ρ=[1.e2]`)                  - weights for condition on regularized variables
* (`iterations::Int64=10`)      - number of outer iterations
* (`iterationsInner::Int64=50`) - maximum number of inner iterations
* (`iterationsCG::Int64=10`)    - maximum number of CG iterations
* (`absTol::Float64=eps()`)     - abs tolerance for stopping criterion
* (`relTol::Float64=eps()`)     - rel tolerance for stopping criterion
* (`tolInner::Float64=1.e-5`)   - tolerance for CG stopping criterion
"""
function SplitBregman(A::matT, x::vecT=zeros(eltype(A),size(A,2)), b=nothing; reg=nothing, regName=["L1"]
                    , λ=[0.0]
                    , regTrafo=nothing
                    , precon=Identity()
                    , ρ=[1.e2]
                    , iterations::Int64=10
                    , iterationsInner::Int64=50
                    , iterationsCG::Int64=10
                    , absTol::Float64=eps()
                    , relTol::Float64=eps()
                    , tolInner::Float64=1.e-6
                    , normalizeReg::Bool=false
                    , kargs...) where {matT<:Trafo, vecT<:AbstractVector}

  if reg == nothing
    reg = Regularization(regName, λ, kargs...)
  end

  if regTrafo == nothing
    regTrafo = [opEye(eltype(x),size(A,2)) for i=1:length(vec(reg))]
  end

  if b==nothing
    y = similar(x,size(A,1))
  else
    y = b
  end

  # operator and fields for the update of x
  op = A'*A
  for i=1:length(vec(reg))
    op += ρ[i]*adjoint(regTrafo[i])*regTrafo[i]
  end
  β = similar(x)
  β_yj = similar(x)
  y_j = similar(x, size(A,1))

  # fields for primal & dual variables
  u = similar(x)
  v = [similar(x, size(regTrafo[i],1)) for i=1:length(vec(reg))]
  vᵒˡᵈ = [similar(v[i]) for i=1:length(vec(reg))]
  b = [similar(v[i]) for i=1:length(vec(reg))]

  # statevariables for CG
  # we store them here to prevent CG from allocating new fields at each call
  statevars = CGStateVariables(zero(u),similar(u),similar(u))

  # convergence parameters
  rk = similar( real.(x), length(vec(reg)) ) #[0.0 for i=1:length(vec(reg))]
  sk = similar(x)
  eps_pri = similar( real.(x), length(vec(reg)) ) # [0.0 for i=1:length(vec(reg))]
  eps_dt = similar(x)

  iter_cnt = 1

  # make sure that ρ is a vector
  if typeof(ρ) <: Real
    ρ_vec = similar(x, real(eltype(x)), 1)
    ρ_vec .= ρ
  else
    ρ_vec = typeof(real.(x))(ρ)
  end

  return SplitBregman(A,vec(reg),regTrafo,y,op,β,β_yj,y_j,u,v,vᵒˡᵈ,b,precon,ρ_vec
              ,iterations,iterationsInner,iterationsCG,statevars, rk,sk
              ,eps_pri,eps_dt,0.0,absTol,relTol,tolInner,iter_cnt,normalizeReg,1.0)
end

"""
  init!(solver::SplitBregman{matT,vecT,opT,rvecT,preconT}, b::vecT
              ; A::matT=solver.A
              , u::vecT=similar(b,0)
              , kargs...) where {matT,vecT,opT,rvecT,preconT}

(re-) initializes the SplitBregman iterator
"""
function init!(solver::SplitBregman{matT,vecT,opT,ropT,rvecT,preconT}, b::vecT
              ; A::matT=solver.A
              , u::vecT=similar(b,0)
              , kargs...) where {matT,vecT,opT,ropT,rvecT,preconT}

  # operators
  if A != solver.A
    solver.A = A
    solver.op = A'*A
    for i=1:length(vec(reg))
      solver.op += ρ[i]*adjoint(regTrafo[i])*regTrafo[i]
    end
  end
  solver.y = b

  # start vector
  if isempty(u)
    if !isempty(b)
      solver.u[:] .= adjoint(A) * b
    else
      solver.u[:] .= 0.0
    end
  else
    solver.u[:] .= u
  end

  # primal and dual variables
  for i=1:length(solver.reg)
    solver.v[i][:] .= solver.regTrafo[i]*solver.u
    solver.vᵒˡᵈ[i][:] .= 0
    solver.b[i][:] .= 0
  end

  # right hand side for the x-update
  solver.y_j[:] .= b
  solver.β_yj .= adjoint(A) * b

  # convergence parameter
  solver.σᵃᵇˢ = sqrt(length(b))*solver.absTol

  # normalization of regularization parameters
  if solver.normalizeReg
    solver.regFac = norm(b,1)/length(b)
  else
    solver.regFac = 1.0
  end

  # reset interation counter
  solver.iter_cnt = 1
end

"""
    solve(solver::SplitBregman, b::vecT; kargs...) where vecT

solves an inverse problem using the Split Bregman method.

# Arguments
* `solver::SplitBregman`              - the solver containing both system matrix and regularizer
* `b::vecT`                           - data vector
* (`A::matT=solver.A`)                - operator for the data-term of the problem
* (`startVector::vecT=similar(b,0)`)  - initial guess for the solution
* (`solverInfo=nothing`)              - solverInfo for logging

when a `SolverInfo` objects is passed, the primal residuals `solver.rk`
and the dual residual `norm(solver.sk)` are stored in `solverInfo.convMeas`.
"""
function solve(solver::SplitBregman, b::vecT; A::matT=solver.A, startVector::vecT=similar(b,0), solverInfo=nothing, kargs...) where {vecT,matT}
  # initialize solver parameters
  init!(solver, b; A=A, u=startVector)

  # log solver information
  solverInfo != nothing && storeInfo(solverInfo,solver.v,solver.rk...,norm(solver.sk))

  # perform SplitBregman iterations
  for (iteration, item) = enumerate(solver)
    solverInfo != nothing && storeInfo(solverInfo,solver.v,solver.rk...,norm(solver.sk))
  end

  return solver.u
end

"""
    splitBregman(A, y::vecT, reg::Vector{Regularization}; kargs...) where vecT

Split Bregman method

Solve the problem: min_x r1(x) + λ*r2(x) such that Ax=b
Here:
  x: variable (vector)
  b: measured data
  A: a general linear operator
  g(X): a convex but not necessarily a smooth function

For details see:
T. Goldstein, S. Osher ,
The Split Bregman Method for l1 Regularized Problems

  # Arguments
  * `A`                           - system matrix
  * `y::vecT`                     - data vector (right hand size)
  * `reg::Vector{Regularization}` - Regularization objects
  * `regTrafo::Vector{Trafo}`     - transformations applied inside each regularizer
  * (`ρ1::Float64=1.0`)           - weighting factor for constraint u=v1
  * (`ρ2::Float64=1.0`)           - weighting factor for constraint u=v2
  * (`precon=Identity()`)         - precondionner to use with CG
  * (`startVector=nothing`)       - start vector
  * (`iterations::Int64=50`)      - maximum number of iterations
  * (`iterationsInner::Int64=50`) - maximum number of inner iterations
  * (`iterationsCG::Int64=10`)    - maximum number of CG iterations
  * (`absTol::Float64=eps()`)     - absolute tolerance for stopping criterion
  * (`relTol::Float64=eps()`)     - relative tolerance for stopping criterion
  * (`tolInner::Float64=1.e-3`)   - relative tolerance for CG
  * (`solverInfo = nothing`)      - `solverInfo` object used to store convergence metrics
"""
function iterate(solver::SplitBregman{matT,vecT,opT,rvecT,preconT}, iteration::Int=1) where {matT,vecT,opT,rvecT,preconT}
  if done(solver, iteration) return nothing end

  # update u
  solver.β[:] .= solver.β_yj
  for i=1:length(solver.reg)
    solver.β[:] .+= solver.ρ[i]*adjoint(solver.regTrafo[i])*(solver.v[i].-solver.b[i])
  end
  cg!(solver.u,solver.op,solver.β,Pl=solver.precon,maxiter=solver.iterationsCG,reltol=solver.tolInner)

  #  proximal map for regularization terms
  for i=1:length(solver.reg)
    copyto!(solver.vᵒˡᵈ[i], solver.v[i])
    solver.v[i][:] .= solver.regTrafo[i]*solver.u .+ solver.b[i]
    if solver.ρ[i] != 0
      solver.reg[i].prox!(solver.v[i],solver.regFac*solver.reg[i].λ/solver.ρ[i]; solver.reg[i].params...)
    end
  end

  # update b
  for i=1:length(solver.reg)
    solver.b[i] .+= solver.regTrafo[i]*solver.u .- solver.v[i]
  end

  # update convergence criteria
  # primal residuals norms (one for each constraint)
  for i=1:length(solver.reg)
    solver.rk[i] = norm(solver.regTrafo[i]*solver.u-solver.v[i])
    solver.eps_pri[i] = solver.σᵃᵇˢ + solver.relTol*max( norm(solver.regTrafo[i]*solver.u), norm(solver.v[i]) )
  end
  # accumulated dual residual
  # effectively this corresponds to combining all constraints into one larger constraint.
  solver.sk[:] .= 0.0
  solver.eps_dt[:] .= 0.0
  for i=1:length(solver.reg)
    solver.sk[:] .+= solver.ρ[i]*adjoint(solver.regTrafo[i])*(solver.v[i].-solver.vᵒˡᵈ[i])
    solver.eps_dt[:] .+= solver.ρ[i]*adjoint(solver.regTrafo[i])*solver.b[i]
  end

  if update_y(solver,iteration)
    solver.y_j[:] .+= solver.y .- solver.A*solver.u
    solver.β_yj[:] .= adjoint(solver.A) * solver.y_j
    # reset v and b
    for i=1:length(solver.reg)
      solver.v[i][:] .= solver.regTrafo[i]*solver.u
      solver.b[i] .= 0
    end
    solver.iter_cnt += 1
    iteration = 0
  end

  return solver.rk[1], iteration+1

end

function converged(solver::SplitBregman)
  if norm(solver.sk) >= solver.σᵃᵇˢ+solver.relTol*norm(solver.eps_dt)
    return false
  else
    for i=1:length(solver.reg)
      (solver.rk[i] >= solver.eps_pri[i]) && return false
    end
  end

  return true
end

@inline done(solver::SplitBregman,iteration::Int) = (iteration==1 && solver.iter_cnt>solver.iterations)

function update_y(solver::SplitBregman,iteration::Int)
  conv = converged(solver)
  return conv || iteration >= solver.iterationsInner
end


text:
using Distributions

function forward(V::Array{Float64,1},A::Array{Float64,2}, mu::Array{Float64,1}, sigma::Real)
    T = length(V)
    nstates = size(A,1)
    a = zeros(nstates, T)
    a[1,1] = 1
    LP = Array(Normal, nstates)
    for i=1:nstates
        LP[i] = Normal(mu[i],sigma)
    end
    aa = 0.0
    Z = 0.0
    for i=2:T
        Z = 0.0
        for j=1:nstates
            b = pdf(LP[j],V[i])
            aa = 0.0
            for k=1:nstates
                aa += a[k,i-1]*A[k,j]
            end
            a[j,i] = b*aa
            Z += a[j,i]
        end
        for j=1:nstates
            a[j,i] /= Z
        end
    end
    return a
end

function backward(V::Array{Float64,1},A::Array{Float64,2}, mu::Array{Float64,1}, sigma::Real)
    T = length(V)
    nstates = size(A,1)
    a = zeros(nstates, T)
    a[1,T] = 1
    LP = Array(Normal, nstates)
    for i=1:nstates
        LP[i] = Normal(mu[i],sigma)
    end
    aa = 0.0
    Z = 0.0
    for i=T-1:-1:1
        Z = 0.0
        for j=1:nstates
            aa = 0.0
            for k=1:nstates
                b = pdf(LP[k],V[i+1])
                aa += a[k,i+1]*A[j,k]
            end
            a[j,i] = aa
            Z += a[j,i]
        end
        for j=1:nstates
            a[j,i] /= Z
        end
    end
    return a
end

function gamma(a,b,A,mu,sigma,x)
    nstates = size(A,1)
    n = size(a,2)
    G = zeros(nstates,nstates,n)
    LP = Array(Normal, nstates)
    for i=1:nstates
        LP[i] = Normal(mu[i],sigma)
    end
    for t=2:n
        for j=1:nstates
            for i=1:nstates
                G[i,j,t] = a[i,t-1]*A[i,j]*b[j,t]*pdf(LP[j],x[t])
            end
        end
    end
    return G
end

function prepA(p,n)
    A = log(zeros(n,n))
    A[1,1] = log(1.-p)
    A[1,2] = log(p)
    for i=2:n-1
        A[i,i+1] = 0
    end
    A[n,1] = 0
    return A
end

function test()
    A = prepA(1e-9)
end


text:
using Gadfly, Colors, BlackBoxOptimizationBenchmarking
using Distributed, Statistics

using BlackBoxOptimizationBenchmarking
const BBOB = BlackBoxOptimizationBenchmarking

bbob_dir = dirname(pathof(BBOB))
include(joinpath(bbob_dir,"../scripts/optimizers_interface.jl"))

## implement interface for IntervalOptimisation
using IntervalArithmetic, IntervalOptimisation

struct IntervalOptimisationMethod
end
string(opt::IntervalOptimisationMethod) = "IntervalOptimisation"

function optimize(NLmeth::IntervalOptimisationMethod,f,D,run_length)
    minf, minx = minimise(f, IntervalBox(-5.5..5.5, D), f_calls_limit=run_length)
    return NLmeth, Vector(mean(mid.(minx))), mid(minf)
end
minimum(  mfit::Tuple{IntervalOptimisationMethod, Vector{Float64}, Float64}) = mfit[3]
minimizer(mfit::Tuple{IntervalOptimisationMethod, Vector{Float64}, Float64}) = mfit[2]

## Setup

dimensions = 3
Ntrials = 30
Δf = 1e-6

c = opt -> Chain(opt, NelderMead(), 0.9)

funcs = 1:14#length(enumerate(BBOBFunction))

run_lengths = round.(Int, range(20, stop = 10_000, length=20))

optimizers = [
    c(BlackBoxOptimMethod(:adaptive_de_rand_1_bin)),
    NelderMead(),
    c(IntervalOptimisationMethod()),
]
opt_strings = map(string, optimizers)

mean_succ, mean_dist, mean_fmin, runtime = BBOB.benchmark(
    optimizers, funcs, run_lengths, Ntrials, dimensions, Δf,
)

cols = Colors.distinguishable_colors(size(mean_succ,1)+1,colorant"white")[2:end]
line_style = i -> mod(i,2)==0 ? [:solid] : [:dash]

## Plots

idx = sortperm([mean(mean_succ[i,:,end,:]) for i=1:length(optimizers)],rev=true)

p = plot(
    [layer(
        x=collect(run_lengths),y=mean(mean_succ[i,:,:,:],dims=(1,3)),Geom.line,
        Theme(default_color=cols[i], line_style=line_style(i), line_width=2pt)
    ) for i in 1:size(mean_succ,1)]...,
    Coord.cartesian(ymax=1),
    Guide.title("D: $(dimensions), $funcs"), Guide.xlabel("Run Length"), Guide.ylabel("Success rate"),
    Guide.manual_color_key("", opt_strings[idx], cols[idx]),
    Theme(background_color="white"),
)

draw(PNG("mean_succ_$(dimensions).png",16cm,10cm),p)

## runtime

idx = sortperm(mean(runtime,dims=2)[:],rev=true)

p = plot(
    layer(
        y=opt_strings[idx], x=(mean(runtime,dims=2)/Base.minimum(mean(runtime,dims=2)))[idx],
        color = cols[idx],
        Geom.bar(orientation=:horizontal)
    ),
    Guide.title("All functions"), Guide.ylabel(""),
    Guide.xlabel("Relative Run Time",orientation=:horizontal),
    Scale.x_log10,
#    Scale.y_log10,
)


text:
# ---------------- logsumexp ------------------

@testset "testing logsumexp3" begin

    i = ShaleDrillingModel._nS(wp)-3
    dmaxp1 = ShaleDrillingModel._dmax(wp,i)+1

    tmpv = dcdp_tmpvars(prim)
    ubV  = @view(tmpv.ubVfull[ :,:,1:dmaxp1])
    dubV = @view(tmpv.dubVfull[:,:,:,1:dmaxp1])

    let geoid = 2,
        roy = 0.25,
        itype = (geoid, roy,)

        @views fillflows!(      ubV,   flow, prim, θt, σv, i, itype...)
        @views fillflows_grad!(dubV, flowdθ, prim, θt, σv, i, itype...)
    end

    tst = zeros(size(tmpv.lse))
    lsetest = zeros(size(tmpv.lse))
    @views qvw =  tmpv.q[:,:,1:dmaxp1]
    @views EV0  = evs.EV[:,:,i]
    @views dEV0 = evs.dEV[:,:,:,i]
    lse = tmpv.lse
    tmp = tmpv.tmp

    lse .= 0.0
    tmp .= 0.0
    @test all(lse.==0.0)
    @test all(tmp.==0.0)

    ShaleDrillingModel.logsumexp3!(lse,tmp,ubV)
    for CI in CartesianIndices(size(lsetest))
        @views lsetest[CI] = logsumexp(ubV[CI,:])
    end
    @show findmax(abs.(lsetest .- lse))
    @test lsetest ≈ lse

    # logsumexp_and_softmax3
    logsumexp_and_softmax3!(lse, qvw, tmp, ubV)
    @test all(sum(qvw, dims=3) .≈ 1.0)
    @show findmax(abs.(lsetest .- lse))
    @test lsetest ≈ lse

    # logsumexp_and_softmax3 - q0
    logsumexp_and_softmax3!(lse, tst, tmp, ubV)
    @test all(tst .== qvw[:,:,1])
    @test all(lsetest .== lse)
end


text:
struct NNPDEHan{C1,C2,O} <: NeuralPDEAlgorithm
    u0::C1
    σᵀ∇u::C2
    opt::O
end
NNPDEHan(u0,σᵀ∇u;opt=Flux.ADAM(0.1)) = NNPDEHan(u0,σᵀ∇u,opt)

function DiffEqBase.solve(
    prob::TerminalPDEProblem,
    alg::NNPDEHan;
    abstol = 1f-6,
    verbose = false,
    maxiters = 300,
    save_everystep = false,
    dt,
    give_limit = false,
    trajectories,
    sdealg = EM(),
    ensemblealg = EnsembleThreads(),
    trajectories_upper = 1000,
    trajectories_lower = 1000,
    maxiters_upper = 10,
    )

    X0 = prob.X0
    ts = prob.tspan[1]:dt:prob.tspan[2]
    d  = length(X0)
    g,f,μ,σ,p = prob.g,prob.f,prob.μ,prob.σ,prob.p

    data = Iterators.repeated((), maxiters)

    #hidden layer
    opt = alg.opt
    u0 = alg.u0
    σᵀ∇u = alg.σᵀ∇u
    ps = Flux.params(u0, σᵀ∇u...)

    function sol()
        map(1:trajectories) do j
            u = u0(X0)[1]
            X = X0
            for i in 1:length(ts)-1
                t = ts[i]
                _σᵀ∇u = σᵀ∇u[i](X)
                dW = sqrt(dt)*randn(d)
                u = u - f(X, u, _σᵀ∇u, p, t)*dt + _σᵀ∇u'*dW
                X  = X .+ μ(X,p,t)*dt .+ σ(X,p,t)*dW
            end
            X,u
        end
    end

    function loss()
        mean(sum(abs2,g(X) - u) for (X,u) in sol())
    end

    iters = eltype(X0)[]

    cb = function ()
        save_everystep && push!(iters, u0(X0)[1])
        l = loss()
        verbose && println("Current loss is: $l")
        l < abstol && Flux.stop()
    end

    Flux.train!(loss, ps, data, opt; cb = cb)

    if give_limit == false
        save_everystep ? iters : u0(X0)[1]
    else
        A = prob.A
        u_domain = prob.u_domain

        ## UPPER LIMIT
        sdeProb = SDEProblem(μ , σ , X0 , prob.tspan)
        ensembleprob = EnsembleProblem(sdeProb)
        sim = solve(ensembleprob, EM(), ensemblealg, dt=dt, trajectories=trajectories_upper,prob.kwargs...)
        function sol_high()
            map(sim.u) do u
                xsde = u.u
                U = g(xsde[end])
                u = u0(X0)[1]
                for i in length(ts):-1:3
                    t = ts[i]
                    _σᵀ∇u = σᵀ∇u[i-1](xsde[i-1])
                    dW = sqrt(dt)*randn(d)
                    U = U .+ f(xsde[i-1], U, _σᵀ∇u, p, t)*dt .- _σᵀ∇u'*dW
                end
                U
            end
        end

        loss_() = sum(sol_high())/trajectories_upper

        ps = Flux.params(u0, σᵀ∇u...)
        cb = function ()
            l = loss_()
            verbose && println("Current loss is: $l")
            l < abstol && Flux.stop()
        end
        dataS = Iterators.repeated((), maxiters_upper)
        Flux.train!(loss_, ps, dataS, ADAM(0.01); cb = cb)
        u_high = loss_()
        ##Lower Limit

        # Function to precalculate the f values over the domain
        function give_f_matrix(X,urange,σᵀ∇u,p,t)
          map(urange) do u
            f(X,u,σᵀ∇u,p,t)
          end
        end

        #The legendre transform that uses the precalculated f values.
        function legendre_transform(f_matrix , a , urange)
            le = a.*(collect(urange)) .- f_matrix
            return maximum(le)
        end

        function sol_low()
            map(1:trajectories_lower) do j
                u = u0(X0)[1]
                X = X0
                I = zero(eltype(u))
                Q = zero(eltype(u))
                for i in 1:length(ts)-1
                    t = ts[i]
                    _σᵀ∇u = σᵀ∇u[i](X)
                    dW = sqrt(dt)*randn(d)
                    u = u - f(X, u, _σᵀ∇u, p, t)*dt + _σᵀ∇u'*dW
                    X  = X .+ μ(X,p,t)*dt .+ σ(X,p,t)*dW
                    f_matrix = give_f_matrix(X , u_domain, _σᵀ∇u, p, ts[i])
                    a_ = A[findmax(collect(A).*u .- collect(legendre_transform(f_matrix, a, u_domain) for a in A))[2]]
                    I = I + a_*dt
                    Q = Q + exp(I)*legendre_transform(f_matrix, a_, u_domain)
                end
                I , Q , X
            end
        end
        u_low = sum(exp(I)*g(X) - Q for (I ,Q ,X) in sol_low())/(trajectories_lower)
        save_everystep ? iters : u0(X0)[1] , u_low , u_high
    end

end #pde_solve


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations2D(1.4)

"""
    initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations2D)

The classical isentropic vortex test case of
- Chi-Wang Shu (1997)
  Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory
  Schemes for Hyperbolic Conservation Laws
  [NASA/CR-97-206253](https://ntrs.nasa.gov/citations/19980007543)
"""
function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations2D)
  # needs appropriate mesh size, e.g. [-10,-10]x[10,10]
  # make sure that the inicenter does not exit the domain, e.g. T=10.0
  # initial center of the vortex
  inicenter = SVector(0.0, 0.0)
  # size and strength of the vortex
  iniamplitude = 0.2
  # base flow
  rho = 1.0
  v1 = 1.0
  v2 = 1.0
  vel = SVector(v1, v2)
  p = 10.0
  rt = p / rho                  # ideal gas equation
  cent = inicenter + vel*t      # advection of center
  cent = x - cent               # distance to centerpoint
  #cent=cross(iniaxis,cent)     # distance to axis, tangent vector, length r
  # cross product with iniaxis = [0,0,1]
  cent = SVector(-cent[2], cent[1])
  r2 = cent[1]^2 + cent[2]^2
  du = iniamplitude/(2*π)*exp(0.5*(1-r2)) # vel. perturbation
  dtemp = -(equations.gamma-1)/(2*equations.gamma*rt)*du^2            # isentrop
  rho = rho * (1+dtemp)^(1\(equations.gamma-1))
  vel = vel + du*cent
  v1, v2 = vel
  p = p * (1+dtemp)^(equations.gamma/(equations.gamma-1))
  prim = SVector(rho, v1, v2, p)
  return prim2cons(prim, equations)
end
initial_condition = initial_condition_isentropic_vortex

volume_flux = flux_shima_etal
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs,
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

coordinates_min = (-10.0, -10.0)
coordinates_max = ( 10.0,  10.0)
refinement_patches = (
  (type="box", coordinates_min=(0.0, -10.0), coordinates_max=(10.0, 10.0)),
)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=4,
                refinement_patches=refinement_patches,
                n_cells_max=10_000,)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true,
                                     extra_analysis_errors=(:conservation_error,),
                                     extra_analysis_integrals=(entropy, energy_total,
                                                               energy_kinetic, energy_internal))

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=1.4)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
@testset "Sorts" begin
    sorts = [
        bubble_sort!
        bucket_sort!
        counting_sort!
        exchange_sort!
        heap_sort!
        insertion_sort!
        merge_sort!
        quick_sort!
        selection_sort!
    ]

    for f in sorts
        x = [3, 5, 1, 4, 2]
        f(x)
        @test x == [1, 2, 3, 4, 5]
    end
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
export malloc

# Stub implementation
malloc(::Csize_t) = C_NULL


text:
"""
Main module for `KirchMig.jl` -- a Julia package for Kirchhoff migration.
"""
module KirchMig

include("migration.jl")
include("map.jl")
include("eikonal.jl")
include("optimization.jl")
include("regularization.jl")

end


text:
using LinearAlgebra
using Statistics

include("algo_util.jl")

abstract type FirstOrder <: DescentDirectionMethod end
function solve(M::FirstOrder, f, ∇f, x, max_iters; num_eval_termination=true)
    init!(M, f, ∇f, x)
    x_hist = [x]
    while true
        x = step!(M, f, ∇f, x)
        push!(x_hist, x)
        if num_eval_termination && (count(∇f) >= max_iters/2) # 2 calls per iteration
            break
        end
    end
    return x, x_hist
end

# Vanilla Gradient Descent
Base.@kwdef mutable struct GradientDescent <: FirstOrder
    α = 1e-3
end
init!(M::GradientDescent, f, ∇, x) = M
function step!(M::GradientDescent, f, ∇f, x)
    α = M.α
    return x - α * ∇f(x)
end

# GD + Momentum
Base.@kwdef mutable struct GDMomentum <: FirstOrder
    α = 1e-3
    β = 0.9
    v = nothing
end
function init!(M::GDMomentum, f, ∇f, x)
    M.v = zeros(length(x))
    return M
end
function step!(M::GDMomentum, f, ∇f, x)
    α, β, = M.α, M.β
    M.v[:] = β * M.v - α * ∇f(x)
    return x + M.v
end

# GD + Nesterov Momentum
Base.@kwdef mutable struct GDMomentumNesterov <: FirstOrder
    α = 1e-3
    β = 0.9
    v = nothing
end
function init!(M::GDMomentumNesterov, f, ∇f, x)
    M.v = zeros(length(x))
    return M
end
function step!(M::GDMomentumNesterov, f, ∇f, x)
    α, β = M.α, M.β
    M.v[:] = β * M.v - α * ∇f(x + β * M.v)
    return x + M.v
end

# Adam
Base.@kwdef mutable struct Adam <: FirstOrder
    α = 1e-3
    β1 = 0.9
    β2 = 0.999
    ϵ = 1e-8
    k = 0
    v = nothing
    s = nothing
end
function init!(M::Adam, f, ∇f, x)
    M.v = zeros(length(x))
    M.s = zeros(length(x))
    return M
end
function step!(M::Adam, f, ∇f, x)
    α, β1, β2, ϵ = M.α, M.β1, M.β2, M.ϵ
    g = ∇f(x)
    M.v[:] = β1 * M.v + (1 - β1) * g
    M.s[:] = β2 * M.s + (1 - β2) * g .* g
    M.k += 1
    v̂ = M.v ./ (1 - β1 ^ M.k)
    ŝ = M.s ./ (1 - β2 ^ M.k)
    return x - (α * v̂) ./ (sqrt.(ŝ) .+ ϵ)
end

# Gradient Descent + Backtracking Line Search
Base.@kwdef mutable struct GDApproxLineSearch <: FirstOrder
    α = 3e-3
    approx_line_search = backtracking_line_search
end
init!(M::GDApproxLineSearch, f, ∇, x) = M
function step!(M::GDApproxLineSearch, f, ∇f, x)
    g = ∇f(x)
    α = M.approx_line_search(f, ∇f, x, -g, M.α)
    return x - α * g
end


text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.

# OptVariable is inspired by the AxisArrays package.
# OptVariable can be replaced with AxisArray once integer indices are no
# longer a special case. See discussions at:
# https://github.com/JuliaArrays/AxisArrays.jl/issues/117
# https://github.com/JuliaArrays/AxisArrays.jl/issues/84
function build_lookup(ax)
    d = Dict{eltype(ax),Int}()
    cnt = 1
    for el in ax
        if haskey(d, el)
            error("Repeated index $el. Index sets must have unique elements.")
        end
        d[el] = cnt
        cnt += 1
    end
    d
end

"""
    OptVariable(cep::OptModelCEP, variable::Symbol, type::String)
Constructor for OptVariable taking JuMP Array and type (ov-operational variable or dv-decision variable)
"""
function OptVariable(cep::OptModelCEP,
                     variable::Symbol,
                     type::String,
                     scale::Dict{Symbol,Int};
                     round_sigdigits::Int=8)
  jumparray=value.(cep.model[variable])
  #Get the scaled optvar from model and turn DenseAxisArray and SparseAxisArray into Dense OptVariable
  scaled_optvar=OptVariable(jumparray,type,cep.set)
  #Unscale the data based on the scaling parameters in Dictionary scaleend
  unscaled_data=scaled_optvar.data*scale[variable]
  #Return Optvariable with the unscaled data
  return OptVariable(round.(unscaled_data;sigdigits=round_sigdigits), scaled_optvar.axes... ;axes_names=scaled_optvar.axes_names ,type=scaled_optvar.type)
end

function OptVariable(jumparray::JuMP.Containers.SparseAxisArray,
                    type::String,
                    set::Dict{String,Dict{String,Array}})
    axes_number=length(first(keys(jumparray.data)))
    axes_names=Array{String,1}()
    for axe_number in 1:axes_number
        #Get the unique values in each dimension of the indexing tuple
        v=unique(getfield.(keys(jumparray.data),axe_number))
        push!(axes_names, get_axes_name(set,v))
    end
    #Get the axes for the optVariable
    axs=get_axes(set, axes_names)
    #Create an OptVariable with zeros
    optvar=OptVariable(zeros(length.(axs)...), axs...; axes_names=axes_names, type=type)
    #Fill data from spare jumparray into dense array
    for idx in eachindex(jumparray)
        #Find corresponding index of dense array
        dense_idx=Array{Any,1}()
        for i in 1:length(idx)
            push!(dense_idx,findfirst(optvar.axes[i].==idx[i]))
        end
        #Overwrite the
        optvar.data[dense_idx...] = jumparray[idx]
    end
    return optvar
end

function OptVariable(jumparray::JuMP.Containers.DenseAxisArray,
                    type::String,
                    set::Dict{String,Dict{String,Array}})
    axes_names=Array{String,1}()
    for axe in jumparray.axes
      for (name, val) in set
          for (n,v) in val
            if axe==v
              push!(axes_names, name)
              break
            end
        end
      end
    end
    return OptVariable(jumparray.data, jumparray.axes...; axes_names=axes_names, type=type)
end

"""
    OptVariable(data::Array{T, N}, axes...) where {T, N}

Construct a OptVariable array with the underlying data specified by the `data` array
and the given axes. Exactly `N` axes must be provided, and their lengths must
match `size(data)` in the corresponding dimensions.

# Example
```jldoctest
julia> cap = OptVariable([1 2; 3 4],["pv", "wind"], 1:2; axes_names=["tech","time"],type="dv")
2-dimensional OptVariable{Int64,2,...} of type dv with index sets:
    Dimension 1 - tech, ["pv", "wind"]
    Dimension 2 - time, 1:2
And data, a 2×2 Array{Int64,2}:
 1  2
 3  4

julia> cap["pv", 2]
2
```
"""
function OptVariable(data::Array{T,N}, axs...;axes_names::Array=repeat([""],N), type="") where {T,N}
    @assert length(axs) == N
    return OptVariable(data, axs, build_lookup.(axs), axes_names, type)
end

"""
    OptVariable{T}(undef, axes...) where T

Construct an uninitialized OptVariable with element-type `T` indexed over the
given axes.

# Example
```jldoctest
julia> cap = OptVariable{Float64}(undef,["pv", "wind"], 1:2; axes_names=["tech","time"],type="dv");

julia> fill!(cap, 1.0)
2-dimensional OptVariable{Float64,2,...} of type dv with index sets:
    Dimension 1 - tech, ["pv", "wind"]
    Dimension 2 - time, 1:2
And data, a 2×2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0

julia> cap["wind", 2] = 5.0
5.0

julia> cap["wind", 2]
5.0

julia> cap
2-dimensional OptVariable{Float64,2,...} of type dv with index sets:
    Dimension 1 - tech, ["pv", "wind"]
    Dimension 2 - time, 1:2
And data, a 2×2 Array{Float64,2}:
 1.0  1.0
 1.0  5.0
```
"""
function OptVariable{T}(::UndefInitializer, axs...; axes_names::Array{String,1}=repeat([""],length(axs)), type="") where T
    return construct_undef_array(T, axs; axes_names=axes_names, type=type)
end

function construct_undef_array(::Type{T}, axs::Tuple{Vararg{Any, N}}; axes_names::Array{String,1}, type::String
                               ) where {T, N}
    return OptVariable(Array{T, N}(undef, length.(axs)...), axs...; axes_names=axes_names, type=type)
end

Base.isempty(A::OptVariable) = isempty(A.data)

# TODO: similar

# AbstractArray interface

Base.size(A::OptVariable) = size(A.data)
Base.LinearIndices(A::OptVariable) = error("OptVariable does not support this operation.")
Base.axes(A::OptVariable) = A.axes
Base.axes(A::OptVariable, dims::String) = A.axes[findfirst(A.axes_names.==dims)]
Base.CartesianIndices(a::OptVariable) = CartesianIndices(a.data)

############
# Indexing #
############

Base.isassigned(A::OptVariable{T,N}, idx...) where {T,N} = length(idx) == N && all(t -> haskey(A.lookup[t[1]], t[2]), enumerate(idx))
# For ambiguity
Base.isassigned(A::OptVariable{T,N}, idx::Int...) where {T,N} = length(idx) == N && all(t -> haskey(A.lookup[t[1]], t[2]), enumerate(idx))

Base.eachindex(A::OptVariable) = CartesianIndices(size(A.data))

lookup_index(i, lookup::Dict) = isa(i, Colon) ? Colon() : lookup[i]

# Lisp-y tuple recursion trick to handle indexing in a nice type-
# stable way. The idea here is that `_to_index_tuple(idx, lookup)`
# performs a lookup on the first element of `idx` and `lookup`,
# then recurses using the remaining elements of both tuples.
# The compiler knows the lengths and types of each tuple, so
# all of the types are inferable.
function _to_index_tuple(idx::Tuple, lookup::Tuple)
    tuple(lookup_index(first(idx), first(lookup)),
          _to_index_tuple(Base.tail(idx), Base.tail(lookup))...)
end

# Handle the base case when we have more indices than lookups:
function _to_index_tuple(idx::NTuple{N}, ::NTuple{0}) where {N}
    ntuple(k -> begin
        i = idx[k]
        (i == 1) ? 1 : error("invalid index $i")
    end, Val(N))
end

# Handle the base case when we have fewer indices than lookups:
_to_index_tuple(idx::NTuple{0}, lookup::Tuple) = ()

# Resolve ambiguity with the above two base cases
_to_index_tuple(idx::NTuple{0}, lookup::NTuple{0}) = ()

to_index(A::OptVariable, idx...) = _to_index_tuple(idx, A.lookup)

# Doing `Colon() in idx` is relatively slow because it involves
# a non-unrolled loop through the `idx` tuple which may be of
# varying element type. Another lisp-y recursion trick fixes that
has_colon(idx::Tuple{}) = false
has_colon(idx::Tuple) = isa(first(idx), Colon) || has_colon(Base.tail(idx))

# TODO: better error (or just handle correctly) when user tries to index with a range like a:b
# The only kind of slicing we support is dropping a dimension with colons
function Base.getindex(A::OptVariable, idx...)
    #if has_colon(idx)
    #    OptVariable(A.data[to_index(A,idx...)...], (ax for (i,ax) in enumerate(A.axes) if idx[i] == Colon())...)
    #else
        return A.data[to_index(A,idx...)...]
    #end
end
Base.getindex(A::OptVariable, idx::CartesianIndex) = A.data[idx]

Base.setindex!(A::OptVariable, v, idx...) = A.data[to_index(A,idx...)...] = v
Base.setindex!(A::OptVariable, v, idx::CartesianIndex) = A.data[idx] = v

Base.IndexStyle(::Type{OptVariable{T,N,Ax}}) where {T,N,Ax} = IndexAnyCartesian()

########
# Keys #
########

"""
    OptVariableKey

Structure to hold a OptVariable key when it is viewed as key-value collection.
"""
struct OptVariableKey{T<:Tuple}
    I::T
end
Base.getindex(k::OptVariableKey, args...) = getindex(k.I, args...)

struct OptVariableKeys{T<:Tuple}
    product_iter::Base.Iterators.ProductIterator{T}
end
Base.length(iter::OptVariableKeys) = length(iter.product_iter)
function Base.eltype(iter::OptVariableKeys)
    return OptVariableKey{eltype(iter.product_iter)}
end
function Base.iterate(iter::OptVariableKeys)
    next = iterate(iter.product_iter)
    return next == nothing ? nothing : (OptVariableKey(next[1]), next[2])
end
function Base.iterate(iter::OptVariableKeys, state)
    next = iterate(iter.product_iter, state)
    return next == nothing ? nothing : (OptVariableKey(next[1]), next[2])
end
function Base.keys(a::OptVariable)
    return OptVariableKeys(Base.Iterators.product(a.axes...))
end
Base.getindex(a::OptVariable, k::OptVariableKey) = a[k.I...]

########
# Show #
########

# Adapted printing from Julia's show.jl

# Copyright (c) 2009-2016: Jeff Bezanson, Stefan Karpinski, Viral B. Shah,
# and other contributors:
#
# https://github.com/JuliaLang/julia/contributors
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.

function Base.summary(io::IO, A::OptVariable)
    _summary(io, A)
    for (k,ax) in enumerate(A.axes)
        print(io, "    Dimension $k - $(A.axes_names[k]), ")
        show(IOContext(io, :limit=>true), ax)
        println(io)
    end
    print(io, "And data, a ", summary(A.data))
end
_summary(io, A::OptVariable{T,N}) where {T,N} = println(io, "$N-dimensional OptVariable{$T,$N,...} of type $(A.type) with index sets:")

function Base.summary(A::OptVariable)
    io = IOBuffer()
    Base.summary(io, A)
    String(io)
end

if isdefined(Base, :print_array) # 0.7 and later
    Base.print_array(io::IO, X::OptVariable{T,1}) where {T} = Base.print_matrix(io, X.data)
    Base.print_array(io::IO, X::OptVariable{T,2}) where {T} = Base.print_matrix(io, X.data)
end

# n-dimensional arrays
function Base.show_nd(io::IO, a::OptVariable, print_matrix::Function, label_slices::Bool)
    limit::Bool = get(io, :limit, false)
    if isempty(a)
        return
    end
    tailinds = Base.tail(Base.tail(axes(a.data)))
    nd = ndims(a)-2
    for I in CartesianIndices(tailinds)
        idxs = I.I
        if limit
            for i = 1:nd
                ii = idxs[i]
                ind = tailinds[i]
                if length(ind) > 10
                    if ii == ind[4] && all(d->idxs[d]==first(tailinds[d]),1:i-1)
                        for j=i+1:nd
                            szj = size(a.data,j+2)
                            indj = tailinds[j]
                            if szj>10 && first(indj)+2 < idxs[j] <= last(indj)-3
                                @goto skip
                            end
                        end
                        #println(io, idxs)
                        print(io, "...\n\n")
                        @goto skip
                    end
                    if ind[3] < ii <= ind[end-3]
                        @goto skip
                    end
                end
            end
        end
        if label_slices
            print(io, "[:, :, ")
            for i = 1:(nd-1); show(io, a.axes[i+2][idxs[i]]); print(io,", "); end
            show(io, a.axes[end][idxs[end]])
            println(io, "] =")
        end
        slice = view(a.data, axes(a.data,1), axes(a.data,2),
                     idxs...)
        Base.print_matrix(io, slice)
        print(io, idxs == map(last,tailinds) ? "" : "\n\n")
        @label skip
    end
end

function Base.show(io::IO, array::OptVariable)
    summary(io, array)
    isempty(array) && return
    println(io, ":")
    Base.print_array(io, array)
end

###################
# Sparse to Dense #
###################

"""
    get_axes(set::Dict{String,Dict{String,Array}}, axes_names::Array{String,1})

get the axes defined by the `axes_names` from the `set`
"""
function get_axes(set::Dict{String,Dict{String,Array}}, axes_names::Array{String,1})
    axes=Array{Array,1}()
    for axes_name in axes_names
        push!(axes,set[axes_name]["all"])
    end
    return axes
end

"""
    get_axes_name(set::Dict{String,Any},values::Array)

Figure out the set within the dictionary `set`, which has equivalent elements to the provided `values`.
The `set` has to be organized as follows: Each entry `set[set_name]` can either be:
- a set-element itself, which is an Array or UnitRange
- or a dictionary with set-subgroups for this set. The set-subgroup has to have a set element called `set[set_name]["all"]`, which contains an Array or UnitRange containing all values for the set_name
"""
function get_axes_name(set::Dict{String,Dict{String,Array}},values::Array)
    for (k,v) in set
        #Check the group `all` first
        if get_axes_name(v["all"],values)==true
            return k
        end
        #Check the subsets second
        for (kk,vv) in v
            if get_axes_name(vv,values)==true
                return k
            end
        end
    end
    return error("The values $values were not found in set $set")
end

function get_axes_name(set_element::Array,values::Array)
    if sort(set_element)==sort(values)
        return true
    end
end

function get_axes_name(set_element::UnitRange,values::Array)
    if collect(set_element)==sort(values)
        return true
    end
end


text:
struct EternityFeaturization <: SeaPearl.AbstractFeaturization end

#this is the best featurization for eternity2 so far.
#it uses FullFeaturedCPNN
#it uses the adhocinfo from the instance generator to instantiate the values' features with the edges' colors.

#when we place a piece on the board, update_features will add the edges' colors to the features.

function SeaPearl.featurize(sr::SeaPearl.DefaultStateRepresentation{EternityFeaturization,TS}) where TS
    n = eternity2_generator.n
    m = eternity2_generator.m
    k = eternity2_generator.k
    g = sr.cplayergraph
    pieces  = g.cpmodel.adhocInfo
    features = zeros(Float32, 6, nv(g))
    for i in 1:nv(g)
        cp_vertex = SeaPearl.cpVertexFromIndex(g, i)
        if isa(cp_vertex, SeaPearl.VariableVertex)
            id = cp_vertex.variable.id
            #println(id)
            prefix, suffix = split(id,'_')
            if prefix=="id"
                a,b = parse(Int, suffix[1]), parse(Int, suffix[2])
                features[1,i] = a
                features[2,i] = b
                features[3,i] = -1 #for now no neighboors
                features[4,i] = -1
                features[5,i] = -1
                features[6,i] = -1
                #features[7,i] = length(cp_vertex.variable.domain)
                #features[8,i] = length(cp_vertex.variable.onDomainChange)
            end
        elseif isa(cp_vertex, SeaPearl.ValueVertex)
            v = cp_vertex.value
            piece = pieces[v,:]
            features[1,i] = piece[1]
            features[2,i] = piece[2]
            features[3,i] = piece[3]
            features[4,i] = piece[4]
        end
    end
    features
end

function SeaPearl.update_features!(sr::SeaPearl.DefaultStateRepresentation{EternityFeaturization,TS}, model::SeaPearl.CPModel) where TS
    g = sr.cplayergraph
    features = sr.nodeFeatures
    for i in 1:nv(g)
        cp_vertex = SeaPearl.cpVertexFromIndex(g, i)
        if isa(cp_vertex, SeaPearl.VariableVertex)
            id = cp_vertex.variable.id
            prefix, suffix = split(id,'_')
            if prefix=="id"
                a,b = parse(Int, suffix[1]), parse(Int, suffix[2])
                var = sr.cplayergraph.cpmodel.variables["src_v"*string(a)*string(b)]
                features[3,i] = SeaPearl.isbound(var) ? SeaPearl.assignedValue(var) : -1
                var = sr.cplayergraph.cpmodel.variables["src_v"*string(a)*string(b+1)]
                features[4,i] = SeaPearl.isbound(var) ? SeaPearl.assignedValue(var) : -1
                var = sr.cplayergraph.cpmodel.variables["src_h"*string(a)*string(b)]
                features[5,i] = SeaPearl.isbound(var) ? SeaPearl.assignedValue(var) : -1
                var = sr.cplayergraph.cpmodel.variables["src_h"*string(a+1)*string(b)]
                features[6,i] = SeaPearl.isbound(var) ? SeaPearl.assignedValue(var) : -1
            end
        end
    end
end

function SeaPearl.feature_length(::Type{SeaPearl.DefaultStateRepresentation{EternityFeaturization, TS}}) where TS
    return 6
end

function SeaPearl.global_feature_length(::Type{SeaPearl.DefaultStateRepresentation{EternityFeaturization, TS}}) where TS
    return 0
end

max_stars_repo_name: corail-research/SeaPearlZoo

max_issues_repo_name: corail-research/SeaPearlZoo

max_forks_repo_name: corail-research/SeaPearlZoo


text:
export Times
"""
    Times(param; ns, Ns)

Constructs a mesh of times, to be used in initial-value problems (see `Problem`).

# Arguments
`param` is either
- a `NamedTuple` containing `dt` the timestep and `T` the final time of comuptation; or
- a vector of computed times.

## Optional keyword arguments
- `ns`  : data are stored every `ns` computations (optional, default = 1).
- `Ns`  : `Ns` data (in addition to the initial datum) are stored (optional, by default `floor( tfin/dt)).
If both `Ns` and `ns` are given, `Ns` overrules `ns`.

# Return values
`t=Times(args)` is of parametric type and offers
- `t.Nc`: number of computed times (including initial datum);
- `t.Ns`: number of stored times (including initial datum);
- `t.ns`: number of computed times between two stored times;
- `t.tfin`: the final time;
- `t.dt`: the timestep;
- `t.tc` : the vector of computed times;
- `t.ts`: the vector of stored times.

"""
struct Times

    Nc   :: Int
    Ns   :: Int
    ns   :: Vector{Int}
    tfin :: Float64
    dt   :: Float64
    tc    :: Vector{Float64}
    ts   :: Vector{Float64}

    function Times( param::NamedTuple; ns = 1, Ns=nothing)
        tc = range(0, stop=param.T, step = param.dt)
        Nc = length(tc)
        if !isnothing(Ns)
            ind_stored = round.(Int,range(1, stop=Nc, length = minimum([Ns+1 Nc])))
        else
            ind_stored = round.(Int,range(1, stop=Nc, step = ns))
        end
        ts = tc[ind_stored]
        Ns = length(ts)
        ns = ind_stored[2:end]-ind_stored[1:end-1]
        new( Nc, Ns, ns, param.T, param.dt, tc, ts)
    end

    function Times( t ; ns = 1, Ns=nothing)
        if !(t[2:end].-t[2]≈t[1:end-1].-t[1]) || t[1]!=0
            @error("Computed times must be equally spaced, and start at the origin.")
        else
            Times((dt=t[2]-t[1],T=t[end]);ns=ns,Ns=Ns)
        end
    end

end

show(io::IO, t::Times) =
    print(io,"Mesh of times on [0, $(t.tfin)], with timestep dt=$(t.dt).\n\
    ├─Number of computed times: $(t.Nc),\n\
    └───Number of stored times: $(t.Ns).")

Base.:(==)(a::Times, b::Times) =
     a.dt == b.dt &&
     a.tfin == a.tfin &&
     a.ns == b.ns &&
     a.Ns == b.Ns &&
     a.Nc == b.Nc &&
     a.tc == b.tc &&
     a.ts == b.ts


text:
using Test, BasicInterpolators
using Random: MersenneTwister

@testset "Accuracy" begin include("test_accuracy.jl") end
@testset "Boundary" begin include("test_boundary.jl") end
@testset "Nevilles" begin include("test_nevilles.jl") end
@testset "Vandermonde" begin include("test_vandermonde.jl") end
@testset "Scattered" begin include("test_scattered.jl") end
@testset "Other" begin include("test_other.jl") end


text:
using Hyperscript
using JSServe, Observables
using JSServe: Application, Session, evaljs, linkjs, update_dom!, div, active_sessions
using JSServe: @js_str, font, onjs, Button, TextField, Slider, JSString, Dependency, with_session

# Javascript & CSS dependencies can be declared locally and
# freely interpolated in the DOM / js string, and will make sure it loads
const THREE = JSServe.Dependency(
    :THREE,
    ["https://cdn.jsdelivr.net/gh/mrdoob/three.js/build/three.min.js"]
)

d = with_session() do session
    width = 500; height = 500
    dom = div(width = width, height = height)
    JSServe.onload(session, dom, js"""
        function (container){
            var renderer = new $(THREE).WebGLRenderer({antialias: true});
            renderer.setSize($width, $height);
            renderer.setClearColor("#ffffff");
            container.appendChild(renderer.domElement);
            var scene = new $THREE.Scene();
            var camera = new THREE.PerspectiveCamera(75, $width / $height, 0.1, 1000);
            camera.position.z = 4;
            var ambientLight = new THREE.AmbientLight( 0xcccccc, 0.4 );
            scene.add( ambientLight );
            var pointLight = new THREE.PointLight( 0xffffff, 0.8 );
            camera.add( pointLight );
            scene.add( camera );
            var geometry = new THREE.SphereGeometry(1.0, 32, 32);
            var material = new THREE.MeshPhongMaterial({color: 0xffff00});
            var sphere = new THREE.Mesh(geometry, material);
            scene.add(sphere);
            renderer.render(scene, camera);
        }
    """)
    return dom
end

# Writing it out to inspect the generated html
open("test.html", "w") do io
    dom = JSServe.jsrender(d.session, d.dom)
    JSServe.dom2html(io, d.session, d.sessionid, dom)
end

# Example of dynamic dome:

color = Observable("red")
with_session() do session
    return font("hi", color = color)
end
# Now you can just update the color wherever you have the observable
color[] = "blue"