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

using Markdown
using InteractiveUtils

# ╔═╡ 943cd3a0-b9a7-11eb-0852-fb73c86a7ebe
using CUDA, Zygote, Optim, Napari, TestImages, LinearAlgebra, FFTW, Noise, LineSearches, Tullio, Statistics, DeconvOptim, FourierTools, Colors, ImageShow, Images

# ╔═╡ 11c1d382-a7b0-4c84-b2fa-ee21ea2da201
md"## Import Packages"

# ╔═╡ 26439a7f-d22c-4b1e-a29c-5124c0a68094
FFTW.set_num_threads(6)

# ╔═╡ 0554827f-6a54-42d6-8a03-fd3a2d1112f5
md"## How does CUDA work?"

# ╔═╡ a14bacc8-a080-430b-a8fc-8084ce7a0e4d
begin
	x = rand(10, 10, 10)
	x_c = CuArray(x)
end;

# ╔═╡ e3d6933a-3ef4-4067-bc6a-570b5f5b1bb0
typeof(x)

# ╔═╡ b5e85f9d-e389-456a-b415-6303cffe7bdb
typeof(x_c)

# ╔═╡ 178b4a21-dd8b-4858-9337-5cf062fe4b09
fft(x);

# ╔═╡ d7347964-eda2-4ae2-80f9-ea2eeca640d8
fft(x_c);

# ╔═╡ fd0311fb-ea86-4e07-b3f9-0038fb539cc8
begin
	# prevent slow scalar indexing on GPU
	CUDA.allowscalar(false);

	# we need to fix some operations so that they are fast o GPUs
	# Reference: https://discourse.julialang.org/t/cuarray-and-optim/14053
	LinearAlgebra.norm1(x::CUDA.CuArray{T,N}) where {T,N} = sum(abs, x); # specializes the one-norm
	LinearAlgebra.normInf(x::CUDA.CuArray{T,N}) where {T,N} = maximum(abs, x); # specializes the one-norm
	Optim.maxdiff(x::CUDA.CuArray{T,N},y::CUDA.CuArray{T,N}) where {T,N} = maximum(abs.(x-y));
end

# ╔═╡ 7fcfc7b7-74b9-436b-9de9-e4b13da0087b
testimage("simple_3d_ball")

# ╔═╡ 8e481367-6fd5-4f63-b4e8-e5f84a70acd5
md"## Load example 3D image"

# ╔═╡ 5a242207-d570-43c7-bfef-20dc02b5a9ea
begin
	img = Float32.(testimage("simple_3d_ball"))
	img = abs.(resample(img, (128, 128, 128)))
	psf = permutedims(ifftshift(Float32.(testimage("simple_3d_psf").parent)), (3, 1, 2))
	psf = abs.(resample(psf, (128, 128, 128)))
	psf ./= sum(psf)
end;

# ╔═╡ b20ca70f-a979-4a53-ba9c-bf3b91dc5261
@view_image Float32.(testimage("simple_3d_psf").parent)

# ╔═╡ c979f867-0830-4e0b-aaca-3ccdc3ac2bc9
md"## Define a Loss Function"

# ╔═╡ eeb36b76-e9bc-4734-95d7-b1296f4050b1
function loss(rec, meas)
	return sum(abs2, rec .- meas)
end

# ╔═╡ 0f29fcc1-90f5-4446-b563-7765de99c177
md"## Create a forward deconvolution model

$$Y = (S * \text{PSF}) + N$$

measurement $Y$, sample $S$,  point spread function $\text{PSF}$, Noise $N$
"

# ╔═╡ d2dbe922-9dc1-4996-b753-62c8f7ff5a6b
begin
	otf = rfft(psf)

	function forward(img)
		return irfft(rfft(img) .* otf, size(img, 1))
	end

	otf_cuda = CuArray(otf)
	p = plan_rfft(CuArray(psf))
	p_inv = inv(p)

	function forward(img::CuArray)
		return p_inv.scale .* (p_inv.p * ((p * img) .* otf_cuda))
	end

	#return conv

end

# ╔═╡ cfc95c98-3f72-4b3a-86d6-6c73806330c6
forward(img);

# ╔═╡ 67e0992e-c504-4631-a510-3adb564a434b
img_c = CuArray(img);

# ╔═╡ 877a5075-0dac-456b-9539-3712abd884e0
forward(img_c);

# ╔═╡ ecf713de-7c69-4fc2-8bf2-00bdff7ddf04
md"## Create a noisy measurement"

# ╔═╡ a4a8e9e9-d68e-43f2-b575-f6b30861753b
begin
	measurement = poisson(forward(img), 50)
	measurement_cuda = CuArray(measurement)
end;

# ╔═╡ 09b5ae2e-0031-44fe-a29d-2c0b0105cb65
@view_image measurement

# ╔═╡ 88498069-35df-489d-9973-57695163fcdb
md"## Regularizer"

# ╔═╡ 7f9b85e3-a3e6-4120-8fd5-ce30debfbf82
reg(rec) = @tullio r = sqrt(1f-8 + abs2(rec[i, j, k] - rec[i + 1, j, k]) + abs2(rec[i, j, k] - rec[i, j + 1, k]) + abs2(rec[i, j, k] - rec[i, j, k + 1]))

# ╔═╡ fdb9ea82-2965-4993-9872-919bc2ed53e8
reg(rec::CuArray) = DeconvOptim.TV_cuda(num_dims=3)(rec)

# ╔═╡ dcc35fc0-b894-42aa-b67e-61bd66f1d7a2
begin
	f(rec) = loss(forward(rec.^2), measurement) + 0.001 * reg(rec.^2)
	f(rec::CuArray) = loss(forward(rec.^2), measurement_cuda) + 0.001 * reg(rec.^2)
	g!(G, rec) = G .= gradient(f, rec)[1]
end

# ╔═╡ 043efd53-e032-4f01-ab13-37a17b9d944c
begin
	rec0 = mean(measurement) .* ones(Float32, size(measurement));
	rec0_cuda = CuArray(rec0);
end;

# ╔═╡ f198cd1d-1fb1-4e00-b135-a69b9c54872b
begin
	f(rec0)
	f(rec0_cuda)
	g!(copy(rec0), rec0)
	g!(copy(rec0_cuda), rec0_cuda)
end;

# ╔═╡ eb9aaea3-b521-4b40-a3c8-04c893c1bf9c
md"## Optimizer"

# ╔═╡ beef10c8-2c24-4be2-8396-3024ff899d09
@time res = Optim.optimize(f, g!, rec0, LBFGS(linesearch=LineSearches.BackTracking()), Optim.Options(iterations=40))

# ╔═╡ f55b0aa2-4de9-43ab-9a30-3d1a7bc85984
CUDA.@time res_cuda = Optim.optimize(f, g!, rec0_cuda, LBFGS(linesearch=LineSearches.BackTracking()), Optim.Options(iterations=40))

# ╔═╡ 214fa627-9dd0-45c1-b766-77cbc2494a7d
begin
	rec = Optim.minimizer(res);
	rec_cuda = Optim.minimizer(res_cuda);
end

# ╔═╡ 4fc1d818-f2e8-428d-bd96-0fbfa9af3954
md"## Viewer"

# ╔═╡ de6df8ce-a6de-4910-a423-ec4e058bd79f
begin
	vr= napari.view_image(rec.^2)
	@add_image vr Array(rec_cuda.^2)
	@add_image vr measurement
	@add_image vr img
	vr.grid.enabled = true
end

# ╔═╡ Cell order:
# ╟─11c1d382-a7b0-4c84-b2fa-ee21ea2da201
# ╠═943cd3a0-b9a7-11eb-0852-fb73c86a7ebe
# ╠═26439a7f-d22c-4b1e-a29c-5124c0a68094
# ╠═0554827f-6a54-42d6-8a03-fd3a2d1112f5
# ╠═a14bacc8-a080-430b-a8fc-8084ce7a0e4d
# ╠═e3d6933a-3ef4-4067-bc6a-570b5f5b1bb0
# ╠═b5e85f9d-e389-456a-b415-6303cffe7bdb
# ╠═178b4a21-dd8b-4858-9337-5cf062fe4b09
# ╠═d7347964-eda2-4ae2-80f9-ea2eeca640d8
# ╟─fd0311fb-ea86-4e07-b3f9-0038fb539cc8
# ╠═7fcfc7b7-74b9-436b-9de9-e4b13da0087b
# ╟─8e481367-6fd5-4f63-b4e8-e5f84a70acd5
# ╠═5a242207-d570-43c7-bfef-20dc02b5a9ea
# ╠═b20ca70f-a979-4a53-ba9c-bf3b91dc5261
# ╟─c979f867-0830-4e0b-aaca-3ccdc3ac2bc9
# ╠═eeb36b76-e9bc-4734-95d7-b1296f4050b1
# ╟─0f29fcc1-90f5-4446-b563-7765de99c177
# ╠═d2dbe922-9dc1-4996-b753-62c8f7ff5a6b
# ╠═cfc95c98-3f72-4b3a-86d6-6c73806330c6
# ╠═877a5075-0dac-456b-9539-3712abd884e0
# ╠═67e0992e-c504-4631-a510-3adb564a434b
# ╠═ecf713de-7c69-4fc2-8bf2-00bdff7ddf04
# ╠═a4a8e9e9-d68e-43f2-b575-f6b30861753b
# ╠═09b5ae2e-0031-44fe-a29d-2c0b0105cb65
# ╠═88498069-35df-489d-9973-57695163fcdb
# ╠═7f9b85e3-a3e6-4120-8fd5-ce30debfbf82
# ╠═fdb9ea82-2965-4993-9872-919bc2ed53e8
# ╠═dcc35fc0-b894-42aa-b67e-61bd66f1d7a2
# ╠═f198cd1d-1fb1-4e00-b135-a69b9c54872b
# ╠═043efd53-e032-4f01-ab13-37a17b9d944c
# ╠═eb9aaea3-b521-4b40-a3c8-04c893c1bf9c
# ╠═beef10c8-2c24-4be2-8396-3024ff899d09
# ╠═f55b0aa2-4de9-43ab-9a30-3d1a7bc85984
# ╠═214fa627-9dd0-45c1-b766-77cbc2494a7d
# ╠═4fc1d818-f2e8-428d-bd96-0fbfa9af3954
# ╠═de6df8ce-a6de-4910-a423-ec4e058bd79f


text:
module ParallelMapAccumulate

using Distributed

export pmapaccumulate

using Base: AsyncGenerator
using Distributed: wrap_on_error, wrap_retry

pmapaccumulate(f, op, p::AbstractWorkerPool, c1, c...; kwargs...) = pmapaccumulate(
    a->f(a...), op, p, zip(c1, c...); kwargs...
)
pmapaccumulate(f, op, c; kwargs...) = pmapaccumulate(
    f, op, default_worker_pool(), c; kwargs...
)
pmapaccumulate(f, op, c1, c...; kwargs...) = pmapaccumulate(
    a->f(a...), op, zip(c1, c...); kwargs...
)

function pmapaccumulate(f, op, p::AbstractWorkerPool, c;
                        distributed=true, on_error=nothing,
                        retry_delays=[], retry_check=nothing)

    # Don't do remote calls if there are no workers.
    if (length(p) == 0) || (length(p) == 1 && fetch(p.channel) == myid())
        distributed = false
    end

    if on_error !== nothing
        f = wrap_on_error(f, on_error)
    end

    if distributed
        f = remote(p, f)
    end

    if length(retry_delays) > 0
        f = wrap_retry(f, retry_delays, retry_check)
    end

    iter = AsyncGenerator(f, c; ntasks = ()->nworkers(p))

    return Iterators.accumulate(op, iter)
end

end # module

ext: jl


text:
@testset "testing simple matrix model" begin
    layerbuilder(k) = Flux.Dense(k, 2, NNlib.relu)
    x = ArrayNode(randn(Float32, 4, 5))
    m = reflectinmodel(x, layerbuilder)
    @test size(m(x).data) == (2, 5)
    @test m isa ArrayModel
    @test eltype(m(x).data) == Float32
end

@testset "testing simple aggregation model" begin
    layerbuilder(k) = Flux.Dense(k, 2, NNlib.relu)
    x = BagNode(ArrayNode(randn(Float32, 4, 4)), [1:2, 3:4])
    m = reflectinmodel(x, layerbuilder)
    @test size(m(x).data) == (2, 2)
    @test m isa BagModel
    @test eltype(m(x).data) == Float32
end

@testset "testing simple tuple models" begin
    layerbuilder(k) = Flux.Dense(k, 2, NNlib.relu)
    x = ProductNode((ArrayNode(randn(Float32, 3, 4)), ArrayNode(randn(Float32, 4, 4))))
    m = reflectinmodel(x, layerbuilder)
    @test eltype(m(x).data) == Float32
    @test size(m(x).data) == (2, 4)
    @test m isa ProductModel
    @test m.ms[1] isa ArrayModel
    @test m.ms[2] isa ArrayModel

    x = ProductNode((BagNode(ArrayNode(randn(Float32, 3, 4)), [1:2, 3:4]),
                  BagNode(ArrayNode(randn(Float32, 4, 4)), [1:1, 2:4])))
    m = reflectinmodel(x, layerbuilder)
    @test size(m(x).data) == (2, 2)
    @test m isa ProductModel
    @test m.ms[1] isa BagModel
    @test m.ms[1].im isa ArrayModel
    @test m.ms[1].bm isa ArrayModel
    @test m.ms[2] isa BagModel
    @test m.ms[2].im isa ArrayModel
    @test m.ms[2].bm isa ArrayModel
end

@testset "testing nested bag model" begin
    bn = BagNode(ArrayNode(randn(Float32, 2, 8)), [1:1, 2:2, 3:6, 7:8])
    x = BagNode(bn, [1:2, 3:4])
    m = reflectinmodel(x, d -> Flux.Dense(d, 2))
    @test size(m(x).data) == (2, 2)
    @test m isa BagModel
    @test m.im isa BagModel
    @test m.im.im isa ArrayModel
    @test m.im.bm isa ArrayModel
    @test m.bm isa ArrayModel
    @test eltype(m(x).data) == Float32

    a = BagNode(BagNode(ArrayNode(randn(2,2)),[1:2]),[1:1])
    b = BagNode(missing,[0:-1])
    c = BagNode(a.data[1:0], [0:-1])
    m = reflectinmodel(a, d -> Dense(d,2), d -> SegmentedMeanMax(d))
    abc = catobs(a, b, c)
    bca = catobs(b, c, a)
    ma = m(a).data
    mb = m(b).data
    mc = m(c).data
    mabc = m(abc).data
    mbca = m(bca).data
    @test mb ≈ mc
    @test mabc[:,1] ≈ ma
    @test mabc[:,2] ≈ mb
    @test mabc[:,3] ≈ mc
    @test mbca[:,1] ≈ mb
    @test mbca[:,2] ≈ mc
    @test mbca[:,3] ≈ ma
end


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export rm_state!
function rm_state!(fst::WFST, idx)
  # map to new indices with (-1) where rm
  new_id = zeros(Int,length(fst))
  s = get_states(fst)

  for i in idx
    new_id[i] = -1
    rm_initial!(fst,i)
    rm_final!(fst,i)
  end

  nstates = 0
  # this pushes all states to be removed at the end of s
  for i in eachindex(s)
    if new_id[i] != -1
      nstates += 1
      new_id[i] = nstates
      if isfinal(fst,i)
        w = get_final_weight(fst,i)
        rm_final!(fst,i)
        final!(fst,nstates,w)
      end
      if isinitial(fst,i)
        w = get_initial_weight(fst,i)
        rm_initial!(fst,i)
        initial!(fst,nstates,w)
      end
      if i != nstates
        s[i], s[nstates] = s[nstates], s[i]
      end
    end
  end
  # remove states
  deleteat!(s,(length(s)-length(idx)+1):length(s))

  # remove arcs and update idx
  for si in s
    to_rm = Int[]
    for (ii, arc) in enumerate(si)
      t = new_id[get_nextstate(arc)]
      if t != -1
        new_arc = change_next(arc,t)
        si[ii] = new_arc
      else
        push!(to_rm,ii)
      end
    end
    deleteat!(si,to_rm)
  end
end


text: Entity(dio, Gl.Text(), Gl.Spatial(), Gl.UniformColor(RGBAf0(1.0, 0.0, 0.0, 1.0)))


text:
"""
$(SIGNATURES)

Sort the vector `v` in-place using `Base.sort!`, then combine all entries which compare equal
(as determined using the `lt` and `by` functions), using the function `combined = combine(x, y)`.

See documentation for [`Base.sort!`](@ref) regarding the `alg`, `by` and `lt` keyword arguments.
"""
function sort_and_combine!(v::AbstractVector;
        alg::Base.Sort.Algorithm=Base.Sort.defalg(v), combine, by=identity, lt=isless)
    # For context, see: https://github.com/JuliaOpt/MathOptInterface.jl/issues/429#issuecomment-406232629.
    isempty(v) && return v
    sort!(v, alg=alg, lt=lt, by=by)
    j = 1
    @inbounds for i in eachindex(v)[2 : end]
        x = v[i]
        if lt(by(v[j]), by(x))
            j += 1
            v[j] = x
        else
            v[j] = combine(v[j], x)
        end
    end
    resize!(v, j)
    v
end


text:
using Test, DiscreteExteriorCalculus, DiscretePDEs
using UniqueVectors: UniqueVector

file_name = joinpath(@__DIR__, "superconductor.geo")
geo_write!(file_name, characteristic_length_factor=.7,
    footer="""
    Box(1) = {-5, -5, -1, 10, 10, 2};
    Box(2) = {-10, -10, -5, 20, 20, 10};
    BooleanDifference{Volume{2}; Delete;}{Volume{1};}
    Point(17) = {10, 0, 0, 1.0};
    Point(18) = {5, 0, 0, 1.0};
    Point(19) = {-5, 0, 0, 1.0};
    Point(20) = {-10, 0, 0, 1.0};
    Line(25) = {17, 18};
    Line(26) = {19, 20};
    Point{17} In Surface{12};
    Point{18} In Surface{2};
    Point{19} In Surface{1};
    Point{20} In Surface{7};
    Curve{25} In Volume{2};
    Curve{26} In Volume{2};
    Physical Curve("1") = {25};
    Physical Curve("2") = {26};
    Physical Volume("Superconductor") = {1};
    """)
initialize!()
gmsh_open!(file_name)
N, K = 3, 4
mesh!(K)

node_tags, points, tcomp = get_triangulated_complex(N, K)
group_dict = get_physical_groups(node_tags, points)
@test typeof(tcomp) <: TriangulatedComplex{N, K}
comp = tcomp.complex
orient!(comp)
sources = [get_current_source(comp, group_dict[k])
    for k in ["1", "2"]]
bound = boundary(comp)
m = Metric(N)
mesh = Mesh(tcomp, circumcenter(m))

μ⁻, Λ⁻ = 2, 0
μ⁻_form = get_material(comp, μ⁻, 3)
Λ⁻_form = get_material(comp, Λ⁻, 2)

bbox, null_basis = magnetostatics_blackbox(m, mesh, sources,
    append!(subcomplex(comp, group_dict["Superconductor"]), bound), μ⁻_form, Λ⁻_form)

if false # set to true to display plots with Gmsh
    A, I = solve_statics(bbox, null_basis, [1.0, 1.0])
    vec_A = sharp(m, comp, A)
    comp_points = UniqueVector([c.points[1] for c in comp.cells[1]])
    ordering = [findfirst(isequal(p), comp_points) for p in points]
    add_field!("Vector potential", node_tags, vec_A[ordering])
    J = source_density(m, mesh, I, 2)
    vec_J = sharp(m, comp, J)
    # current on superconductor
    sub_points = UniqueVector([c.points[1] for c in subcomplex(comp, group_dict["Superconductor"]).cells[1]])
    sub_inds = findall([p in sub_points for p in points])
    add_field!("Current density superconductor", node_tags[sub_inds], vec_J[ordering[sub_inds]])
    # current on boundary
    sub_points = UniqueVector([c.points[1] for c in bound.cells[1]])
    sub_inds = findall([p in sub_points for p in points])
    add_field!("Current density boundary", node_tags[sub_inds], vec_J[ordering[sub_inds]])
    gui!()
end
rm(file_name)


text:
using Documenter, Arborist

makedocs(;
    modules=[Arborist],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/oxinabiox/Arborist.jl/blob/{commit}{path}#L{line}",
    sitename="Arborist.jl",
    authors="Lyndon White",
    assets=[],
)

deploydocs(;
    repo="github.com/oxinabiox/Arborist.jl",
)


text:
using DiscreteExteriorCalculus: Metric, Simplex, Cell, CellComplex, Mesh, sharp,
    subcomplex, orient!, differential_operator_sequence, differential_operator,
    zero_constraint, constant_constraint
using SparseArrays: sparse, spdiagm
using AdmittanceModels: Blackbox, PSOModel, apply_transform, sparse_nullbasis,
    impedance_matrices, admittance_matrices

export electrostatics_blackbox
"""
    electrostatics_blackbox(m::Metric{N}, mesh::Mesh{N, K},
        sources::AbstractVector{<:AbstractVector{Cell{N}}},
        boundary_points::AbstractVector{Cell{N}},
        ϵ::AbstractVector{<:Real}) where {N, K}
    electrostatics_blackbox(m::Metric{N}, mesh::Mesh{N, K},
        sources::AbstractVector{<:AbstractVector{Cell{N}}},
        boundary::CellComplex{N},
        ϵ::AbstractVector{<:Real}) where {N, K}

Each source is a vector of cells with `K == 1` representing a conductor. Discretize the
Poisson equation to find a Blackbox that maps charge to voltage as well as a sparse
constraint matrix that asserts that the voltage is constant on each source.
"""
function electrostatics_blackbox(m::Metric{N}, mesh::Mesh{N, K},
    sources::AbstractVector{<:AbstractVector{Cell{N}}},
    boundary_points::AbstractVector{Cell{N}},
    ϵ::AbstractVector{<:Real}) where {N, K}
    k = 1
    ★★, d₁, ★, d₀ = differential_operator_sequence(m, mesh, "★★d★d", k, true)
    Y = -★★ * d₁ * ★ * spdiagm(0 => ϵ) * d₀
    comp = mesh.primal.complex
    row_inds = [findfirst(isequal(s[1]), comp.cells[k]) for s in sources]
    n = length(sources)
    P = sparse(row_inds, 1:n, ones(n), length(comp.cells[k]), n)
    bbox = Blackbox([0.0], [Y], P, P, collect(1:n))
    conductors = filter(s -> length(s) > 1, sources)
    constraint = vcat(zero_constraint(comp, boundary_points, k),
        [constant_constraint(comp, s, k) for s in conductors]...)
    return bbox, sparse_nullbasis(constraint)
end

electrostatics_blackbox(m::Metric{N}, mesh::Mesh{N, K},
    sources::AbstractVector{<:AbstractVector{Cell{N}}},
    boundary::CellComplex{N},
    ϵ::AbstractVector{<:Real}) where {N, K} =
    electrostatics_blackbox(m, mesh, sources, boundary.cells[1], ϵ)

export get_charge_source
"""
    get_charge_source(comp::CellComplex{N},
        group::AbstractVector{Simplex{N, K}}) where {N, K}

Return a vector of cells representing the points in the given simplices.
"""
get_charge_source(comp::CellComplex{N}, group::AbstractVector{Simplex{N, K}}) where {N, K} =
    subcomplex(comp, group).cells[1]

export magnetostatics_blackbox
"""
    magnetostatics_blackbox(m::Metric{N}, mesh::Mesh{N, K},
        sources::AbstractVector{<:AbstractVector{Cell{N}}},
        boundary_points::AbstractVector{Cell{N}},
        boundary_edges::AbstractVector{Cell{N}},
        μ⁻::AbstractVector{<:Real}, Λ⁻::AbstractVector{<:Real}) where {N, K}
    magnetostatics_blackbox(m::Metric{N}, mesh::Mesh{N, K},
        sources::AbstractVector{<:AbstractVector{Cell{N}}},
        boundary::CellComplex{N}, μ⁻::AbstractVector{<:Real},
        Λ⁻::AbstractVector{<:Real}) where {N, K}

Each source is a vector of cells with `K == 2` representing a current line. In order to get
correct results, the current specified must have divergence 0. Discretize the Maxwell
equations in the Coulomb gauge to find a Blackbox mapping current to the line integral of
the vector potential (for a closed loop this is magnetic flux) as well as a sparse
constraint matrix that asserts that the vector potential is orthogonal to the boundary
and has divergence 0 in the interior.
"""
function magnetostatics_blackbox(m::Metric{N}, mesh::Mesh{N, K},
    sources::AbstractVector{<:AbstractVector{Cell{N}}},
    boundary_points::AbstractVector{Cell{N}},
    boundary_edges::AbstractVector{Cell{N}},
    μ⁻::AbstractVector{<:Real}, Λ⁻::AbstractVector{<:Real}) where {N, K}
    pso, null_basis = electrodynamics_pso(m, mesh, sources, boundary_points,
        boundary_edges, μ⁻, Λ⁻, zeros(size(Λ⁻)), ones(size(Λ⁻)))
    bbox = Blackbox([0.0], [pso.K], pso.P, pso.Q, pso.ports)
    return bbox, null_basis
end

magnetostatics_blackbox(m::Metric{N}, mesh::Mesh{N, K},
    sources::AbstractVector{<:AbstractVector{Cell{N}}},
    boundary::CellComplex{N}, μ⁻::AbstractVector{<:Real},
    Λ⁻::AbstractVector{<:Real}) where {N, K} = magnetostatics_blackbox(m, mesh,
    sources, boundary.cells[1], boundary.cells[2], μ⁻, Λ⁻)

export get_current_source
"""
    get_current_source(comp::CellComplex{N},
        group::AbstractVector{Simplex{N, K}}) where {N, K}

Return a vector of cells representing the lines in the given simplices.
"""
get_current_source(comp::CellComplex{N}, group::AbstractVector{Simplex{N, K}}) where {N, K} =
    orient!(subcomplex(comp, group).cells[2])

export electrodynamics_pso
"""
    electrodynamics_pso(m::Metric{N}, mesh::Mesh{N, K},
        sources::AbstractVector{<:AbstractVector{Cell{N}}},
        boundary_points::AbstractVector{Cell{N}},
        boundary_edges::AbstractVector{Cell{N}},
        μ⁻::AbstractVector{<:Real}, Λ⁻::AbstractVector{<:Real},
        σ::AbstractVector{<:Real}, ϵ::AbstractVector{<:Real}) where {N, K}
    electrodynamics_pso(m::Metric{N}, mesh::Mesh{N, K},
        sources::AbstractVector{<:AbstractVector{Cell{N}}},
        boundary::CellComplex{N},
        μ⁻::AbstractVector{<:Real}, Λ⁻::AbstractVector{<:Real},
        σ::AbstractVector{<:Real}, ϵ::AbstractVector{<:Real}) where {N, K}

Each source is a vector of cells with `K == 2` representing a current line. In order to get
correct results, the current specified must have divergence 0. Discretize the Maxwell
equations in the Coulomb gauge to find a Positive Second Order model mapping current to
voltage as well as a sparse constraint matrix that asserts that the vector potential is
orthogonal to the boundary and has divergence 0 in the interior.
"""
function electrodynamics_pso(m::Metric{N}, mesh::Mesh{N, K},
    sources::AbstractVector{<:AbstractVector{Cell{N}}},
    boundary_points::AbstractVector{Cell{N}},
    boundary_edges::AbstractVector{Cell{N}},
    μ⁻::AbstractVector{<:Real}, Λ⁻::AbstractVector{<:Real},
    σ::AbstractVector{<:Real}, ϵ::AbstractVector{<:Real}) where {N, K}
    k = 2
    ★★, d₁, ★, d₀ = differential_operator_sequence(m, mesh, "★★d★d", k, true)
    K₀ = ★★ * d₁ * ★ * spdiagm(0 => μ⁻) * d₀
    d, ★ = differential_operator_sequence(m, mesh, "d★", k, true)
    K₁ = ★ * spdiagm(0 => Λ⁻)
    G = ★ * spdiagm(0 => σ)
    C = ★ * spdiagm(0 => ϵ)
    comp = mesh.primal.complex
    row_inds = vcat([[findfirst(isequal(c), comp.cells[k]) for c in s] for s in sources]...)
    col_inds = vcat([[i for _ in s] for (i, s) in enumerate(sources)]...)
    P = sparse(row_inds, col_inds, ones(length(row_inds)), length(comp.cells[k]), length(sources))
    pso = PSOModel(K₀ + K₁, G, C, P, P, collect(1:length(sources)))
    interior_inds = findall(i -> !(comp.cells[1][i] in boundary_points), 1:length(comp.cells[1]))
    d★ϵ_constraint = (d * C/maximum(abs.(ϵ)))[interior_inds, :] # only enforce d★ϵA = 0 in the interior
    constraint = vcat(zero_constraint(comp, boundary_edges, k), d★ϵ_constraint)
    return pso, sparse_nullbasis(constraint)
end

electrodynamics_pso(m::Metric{N}, mesh::Mesh{N, K},
    sources::AbstractVector{<:AbstractVector{Cell{N}}},
    boundary::CellComplex{N},
    μ⁻::AbstractVector{<:Real}, Λ⁻::AbstractVector{<:Real},
    σ::AbstractVector{<:Real}, ϵ::AbstractVector{<:Real}) where {N, K} =
    electrodynamics_pso(m, mesh, sources, boundary.cells[1], boundary.cells[2], μ⁻, Λ⁻, σ, ϵ)

export solve_statics
"""
    solve_statics(bbox::Blackbox, null_basis::AbstractMatrix{<:Real},
        source_magnitudes::AbstractVector{<:Real})

Given a Blackbox and constraint matrix from `electrostatics_blackbox` or
`magnetostatics_blackbox`, return the scalar potential and charge or vector potential and
current, respectively.
"""
function solve_statics(bbox::Blackbox, null_basis::AbstractMatrix{<:Real},
    source_magnitudes::AbstractVector{<:Real})
    constrained_bbox = apply_transform(bbox, null_basis)
    field = null_basis * (collect(constrained_bbox.Y[1]) \ collect(constrained_bbox.P * source_magnitudes))
    source = bbox.Y[1] * field
    return field, source
end

export impedance_matrix
"""
    impedance_matrix(bbox::Blackbox, null_basis::AbstractMatrix{<:Real})

Given a Blackbox and constraint matrix from `electrostatics_blackbox` or
`magnetostatics_blackbox`, return the elastance matrix or inductance matrix, respectively.
"""
impedance_matrix(bbox::Blackbox, null_basis::AbstractMatrix{<:Real}) =
    impedance_matrices(apply_transform(bbox, null_basis))[1]

export admittance_matrix
"""
    admittance_matrix(bbox::Blackbox, null_basis::AbstractMatrix{<:Real})

Given a Blackbox and constraint matrix from `electrostatics_blackbox` or
`magnetostatics_blackbox`, return the capacitance matrix or inverse inductance matrix,
respectively.
"""
admittance_matrix(bbox::Blackbox, null_basis::AbstractMatrix{<:Real}) =
    admittance_matrices(apply_transform(bbox, null_basis))[1]

export electric_field
"""
    electric_field(m::Metric{N}, mesh::Mesh{N}, φ::AbstractVector{<:Real}) where N

Given the scalar potential from `solve_statics`, approximate the electric field at each
primal vertex.
"""
electric_field(m::Metric{N}, mesh::Mesh{N}, φ::AbstractVector{<:Real}) where N =
    sharp(m, mesh.primal.complex, -differential_operator(m, mesh, "d", 1, true, φ))

export magnetic_field
"""
    magnetic_field(m::Metric{N}, mesh::Mesh{N}, A::AbstractVector{<:Real}) where N

Given the scalar potential from `solve_statics`, approximate the magnetic field at each
dual vertex.
"""
magnetic_field(m::Metric{N}, mesh::Mesh{N}, A::AbstractVector{<:Real}) where N =
    sharp(m, mesh.dual.complex, differential_operator(m, mesh, "★d", 2, true, A))

export source_density
"""
    source_density(m::Metric{N}, mesh::Mesh{N, K},
        source::AbstractVector{<:Real}, k::Int) where {N, K}

Given the charge or current from `solve_statics`, find the charge density or current
density, respectively.
"""
function source_density(m::Metric{N}, mesh::Mesh{N, K},
    source::AbstractVector{<:Real}, k::Int) where {N, K}
    return differential_operator(m, mesh, "★★★", K-k+1, false, source)
end

export get_material
"""
    get_material(comp::CellComplex{N}, group::AbstractVector{Simplex{N, K}},
        value::Real, k::Int) where {N, K}

Return a vector representing the material distribution that is 0 everywhere except on the
simplices in `group` where it is `value`.
"""
function get_material(comp::CellComplex{N}, group::AbstractVector{Simplex{N, K}},
    value::Real, k::Int) where {N, K}
    material = zeros(length(comp.cells[k]))
    for c in subcomplex(comp, group).cells[k]
        material[findfirst(isequal(c), comp.cells[k])] = value
    end
    return material
end

"""
    get_material(comp::CellComplex, value::Real, k::Int)

Return a vector representing the material distribution that is `value` everywhere.
"""
get_material(comp::CellComplex, value::Real, k::Int) = ones(length(comp.cells[k])) * value

const speed_of_light = 299792458.0
const μ₀ = 4π/1e7
const ϵ₀ = 1/(μ₀ * speed_of_light^2)


text:
module Dataset
export TrainTest, get_MINST

using MLDatasets

struct TrainTest
  train_x
  train_y
  test_x
  test_y
end

function get_MINST()
  test_x, test_y = MNIST.testdata()
  train_x, train_y = MNIST.traindata()
  n_sample_train = 60000
  n_sample_test = 10000
  x_train_whcn = reshape(float.(train_x[:,:,1:n_sample_train]), 28, 28, 1, n_sample_train)
  x_test_whcn = reshape(float.(test_x[:,:,1:n_sample_test]), 28, 28, 1, n_sample_test)

  return TrainTest(
    float.(x_train_whcn .> 0.5),
    train_y[1:n_sample_train],
    float.(x_test_whcn .> 0.5),
    test_y[1:n_sample_test])
end

# take a smaller sample of
function get_MINST(n_sample :: Int)
  @assert n_sample < 10000
  test_x, test_y = MNIST.testdata()
  train_x, train_y = MNIST.traindata()

  x_train_whcn = reshape(float.(train_x[:,:,1:n_sample]), 28, 28, 1, n_sample)
  x_test_whcn = reshape(float.(test_x[:,:,1:n_sample]), 28, 28, 1, n_sample)

  return TrainTest(
    float.(x_train_whcn .> 0.5),
    train_y[1:n_sample],
    float.(x_test_whcn .> 0.5),
    test_y[1:n_sample])
end

function split_X_by_batches(
  X :: Array{T, 4},
  batch_size :: Int64
  ) where {T}
  """
  split_X_by_batches(X, batch_size)
  given X is an Array in (width, height, channel, samples)
  we will further slice X into batches and return a result
  directly passable to Flux.Train!
  """
  N_total = size(X, 4)
  batch_idx = collect(Base.Iterators.partition(1:N_total, batch_size))
  return zip([X[:,:,:,idx] for idx in batch_idx])
end

function repr_digit(
  X :: U,
  y :: Array{T, 1},
  digit :: T
) where {T, U <: AbstractArray}
  @assert digit in unique(y) #ugh, how can this be a total fn?

  digit_idx = findall(y .== digit)
  Xdigit = X[:,:,:, digit_idx]
  n_examples = size(Xdigit,4)
  mat = reshape(Xdigit, n_examples, 28*28)
  dist_mat = zeros(n_examples, n_examples)
  distance(x,y) = sum((x .- y) .^ 2)
  for i in range(1, stop = n_examples)
    for j in range(1, stop = n_examples)
      dist_mat[i,j] = distance(mat[i,:], mat[j,:])
    end
  end

  avg_dist = reshape(sum(dist_mat, dims = 1), n_examples)
  digit_idx[findall(avg_dist .== minimum(avg_dist))[1]]
end

function min_distance_digits()
  """
  min_distance_digits
  Returns the index of a sample set of digits 0-9 in the test set
  with the condition that the example for each digit is the digit
  with the minimum average distance (L2) from that digit and all
  other digits of the same numeral.
  """
  test_x, test_y = MNIST.testdata()
  x_test_whcn = reshape(float.(test_x[:,:,:]), 28, 28, 1, :)

  [repr_digit(x_test_whcn, test_y, i) for i in sort(unique(test_y))]
end

function test_viz_idx()
  """
  test_viz_idx()
  returns the results of min_distance_digits for 0-9
  when applied to MINST test set
  """
  return [
    9426, # 0
    2380,
    1304,
    9008,
    766,
    711,
    2472,
    7761,
    7756,
    7953 # 9
  ]
end

end

max_stars_repo_name: adamwespiser/variational-autoencoders

max_issues_repo_name: adamwespiser/variational-autoencoders

max_forks_repo_name: adamwespiser/variational-autoencoders


text:
##################################################
"""
    ref_map_dict(fmap::String)
---
This create a dictionary of SNP -> [chr, bp] and return it

## Notes
The 3 columns in `fmap`:

1. chromosome
2. SNP-name
3. base pair position
"""
function ref_map_dict(fmap)     # fmap: the physical map
    snpdic = Dict()
    for line in eachline(fmap)
        snp, chr, bp = split(line)
        snpdic[snp] = [chr, bp]
    end

    return snpdic
end

################################################################################
"""
    create_plink_map(sample::String, dic::Dict, to::String)
---
Given a `sample` in final report format and a SNP -> [chr, bp] dictionary, `dic`,
which is a super set of SNP in the final report, this subroutine write a plink
map to `to`.
"""
function create_plink_map(sample::String, dic::Dict, to::String)
    open(to, "w") do foo
        open(sample, "r") do sp
            for line in eachline(sp)
                occursin("SNP Name", line) && break
            end
            for line in eachline(sp)
                snp = split(line)[1]
                chr, bp = dic[snp]
                write(foo, chr, ' ', snp, " 0 ", bp, '\n')
            end
        end
    end
end

##################################################
"""
    clean_all(; force=false)
---
Remove all intermediate data results.  Take care to call this function
"""
function clean_all(; force::Bool = false)
    list = ["data/maps/origin",
            "data/maps/updated",
            "data/genotypes/step-1.plk"]
    if force
        for l in list
            item("Removing '$l'")
            rm(l, force=true, recursive=true)
            done()
        end
    else
        title("Remove all intermediate results")
        message("Are you sure to remove all intermediate results?\n" *
                "If so, call `clean_all(force = true)\n" *
                "And these directories will be removed:")
        for l in list
            println("- $l")
        end
    end
end

"""
    mkdir_tmp(dir::AbstractString = "tmp")
---
- If not exists, make it.
- If exists, clean all files and folders in it.
"""
function mkdir_tmp(dir::AbstractString = "tmp")
    isdir(dir) && rm(dir, force=true, recursive=true)
    mkdir(dir)
end

"""
Create SNP Set from a bim file
"""
function get_bim_snp_set(bim::AbstractString)
    ssnp = Set()
    vsnp = String[]
    for line in eachline(bim)
        push!(vsnp, split(line)[2])
    end
    ssnp = Set(vsnp)
    return ssnp
end


text:
using ADSeismic
using ADCME
using PyPlot
using DelimitedFiles
using JLD2
r = mpi_init()

M = N = Int(round(sqrt(mpi_size())))
n = 200÷M
@info M
param = MPIAcousticPropagatorParams(NX = n*M, NY = n*N, n = n,
     Rcoef=0.2, DELTAX=20, DELTAY=20, DELTAT=0.05)
compute_PML_Params!(param)

srci = [n*param.M÷5]
srcj = [n*param.N÷2]
srcv = hcat([Ricker(param, 100.0, 500.0)]...)
src = MPIAcousticSource(param, srci,srcj,srcv)

rcvj = Array(20:n*param.M-19)
rcvi = ones(Int64, length(rcvj))*((n*param.N)÷5)
receiver =  MPIAcousticReceiver(param, rcvi, rcvj)

c = 1000*ones(n, n)
center = n*M÷2
width = n*M ÷ 8
for i = 1:n
    for j = 1:n
        ii = i + (param.II-1) * n
        jj = j + (param.JJ-1) * n
        x = i*param.DELTAX + (param.II-1) * n * param.DELTAX
        y = j*param.DELTAY + (param.JJ-1) * n * param.DELTAY
        if center -width <= ii <= center + width && center -width <=jj<=center +width
            c[i, j] = 2000.
        end
    end
end

propagator = MPIAcousticPropagatorSolver(param, src, c)
MPISimulatedObservation!(propagator, receiver)
sess = Session(); init(sess)

U = run(sess, propagator.u)
rcvv = run(sess, receiver.rcvv)

if mpi_size()==1
    close("all")
    pcolormesh(rcvv, cmap="gray")
    colorbar()
    xlabel("Location")
    ylabel("Time")
    savefig("receiver.png")

    close("all")
    p = visualize_wavefield(U, param)
    saveanim(p, "wavefield.gif")

    close("all")
    visualize_model(c, param)
    savefig("model.png")
    xlabel("x")
    ylabel("y")
    savefig("model.pdf")
end

@info "Saving..."
@save "data/dat$(mpi_rank())-$(mpi_size()).jld2" U rcvv

mpi_finalize()


text:
using FactCheck, Images, Color, FixedPointNumbers

facts("Overlay") do
    gray = linspace(0.0, 1.0, 5)
    context("One") do
        ovr = Images.Overlay((2gray, 2gray), (RGB(1, 0, 0), RGB(0, 0, 1)), (Clamp{Float64}(), Clamp{Float64}()))
        @fact ovr[1] --> RGB(0, 0, 0)
        @fact ovr[2] --> RGB(0.5, 0, 0.5)
        @fact ovr[3] --> ovr[4]
        @fact ovr[4] --> ovr[5]
        @fact ovr[5] --> exactly(RGB(1, 0, 1))
        @fact eltype(ovr) --> RGB{Float64}
        @fact length(ovr) --> 5
        @fact size(ovr) --> (5,)
        @fact size(ovr, 1) --> 5
        @fact size(ovr, 2) --> 1
        @fact nchannels(ovr) --> 2
        iob = IOBuffer()
        show(iob, ovr)  # exercise only
    end

    context("Two") do
        local ovr = Images.Overlay((gray, 0*gray), (RGB{Ufixed8}(1, 0, 1), RGB{Ufixed8}(0, 1, 0)), ((0, 1), (0, 1)))
        @fact eltype(ovr) --> RGB{Ufixed8}
        s = similar(ovr)
        @fact isa(s, Vector{RGB{Ufixed8}}) --> true
        @fact length(s) --> 5
        s = similar(ovr, RGB{Float32})
        @fact isa(s, Vector{RGB{Float32}}) --> true
        @fact length(s) --> 5
        s = similar(ovr, RGB{Float32}, (3, 2))
        @fact isa(s, Matrix{RGB{Float32}}) --> true
        @fact size(s) --> (3, 2)
        buf = Images.uint32color(ovr)
        gray8 = round(UInt8, 255*gray)
        nogreen = [convert(UInt32, g)<<16 | convert(UInt32, g) for g in gray8]
        @fact buf --> nogreen
    end

    context("Three") do
        ovr = Images.Overlay((gray, 0*gray), [RGB(1, 0, 1), RGB(0, 1, 0)], ([0, 1], [0, 1]))
        @fact_throws ErrorException Images.Overlay((gray, 0*gray), (RGB(1, 0, 1), RGB(0, 1, 0)), ((0,), (0, 1)))
        @fact_throws ErrorException ovr[1] = RGB(0.2, 0.4, 0.6)
    end

    context("Four") do
        img1 = Images.Image(gray, Dict{ASCIIString, Any}([("colorspace", "Gray"), ("spatialorder",["x"])]))
        ovr = Images.OverlayImage((2gray, img1), (RGB{Float32}(1, 0, 1), RGB{Float32}(0, 1, 0)), ((0, 1),(0, 1)))
        @fact isa(ovr, Images.Image) --> true
        @fact haskey(ovr, "colorspace") --> false
        @fact Images.colorspace(ovr) --> "RGB"
        @fact ovr[2] --> RGB{Float32}(0.5, 0.25, 0.5)
        a = rand(Float32, 3, 2)
        b = rand(Float32, 3, 2)
        ovr = Images.OverlayImage((a, b), (RGB{Float32}(1, 0, 1), RGB{Float32}(0, 1, 0)), ((0, 1), (0, 1)))
        @fact isa(ovr, Images.Image) --> true
        @fact abs(ovr[1, 2] - RGB{Float32}(a[1, 2], b[1, 2], a[1, 2])) --> roughly(0, atol=1e-5)
    end
end


text:
using Distributed
addprocs(16)
@everywhere include("../src/get_sens_solenoid.jl")
si = LinRange(0.,0.,1)
s = zeros(3,size(si)[1])
s[1,:] .= 1.0
s[3,:] .= si
s[2,:] .= 4.0
get_sens(s)

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
abstract type AbstractFeedForwardAffect end

struct UpperBoundFF <: AbstractFeedForwardAffect
    vars_prefix::Vector{Symbol}
    cache::Union{Nothing, Type{<:AbstractCache}}
end

UpperBoundFF(var::Symbol,
             cache::Union{Nothing, Type{<:AbstractCache}}=nothing) = UpperBoundFF([var], cache)

get_vars_prefix(p::UpperBoundFF) = p.vars_prefix

struct RangeFF <: AbstractFeedForwardAffect
    lb_vars_prefix::Vector{Symbol}
    ub_vars_prefix::Vector{Symbol}
    cache::Union{Nothing, Type{<:AbstractCache}}
end

RangeFF(var_lb::Symbol,
        var_ub::Symbol,
        cache::Union{Nothing, Type{<:AbstractCache}}=nothing) = RangeFF([var_lb], [var_ub], cache)

get_vars_prefix(p::RangeFF) = (p.lb_var_prefix, p.lb_var_prefix)

struct SemiContinuousFF <: AbstractFeedForwardAffect
    vars_prefix::Vector{Symbol}
    bin_prefix::Symbol
    cache::Union{Nothing, Type{<:AbstractCache}}
end

SemiContinuousFF(var::Symbol,
    bin_var::Symbol,
    cache::Union{Nothing, Type{<:AbstractCache}}=nothing) = SemiContinuousFF([var], bin_var, cache)

get_bin_prefix(p::SemiContinuousFF) = p.bin_prefix
get_vars_prefix(p::SemiContinuousFF) = p.vars_prefix


text:
# 多次元対応？

# 確率密度関数を求める
function calc_densities(p, data, sigma)
    density = []
    for i = 1:length(p)
        push!(density, calc_density(p[i], data, sigma))
    end

    return density
end

# 各点の密度を求める
function calc_density(p, data, sigma)
    density = 0.0
    for i = 1:length(data)

        # 各次元のガウス値を掛け合わせる
        gaussian = 1.0
        for j = 1:length(data[i])
            gaussian *= get_gaussian(p[j], data[i][j], sigma)
        end

        density += gaussian
    end

    density /= length(data)
    return density
end

# ガウス関数の値を求める
function get_gaussian(dx, mean, sigma)
    return exp((-(mean - dx)^2)/(2*(sigma^2)))/(sqrt(2.0*pi)*sigma)

    # sigma_matrix = (sigma^2) * collect(Diagonal(ones(length(dx))))

    # gaussian = exp(-(1/2) * ((dx - mean) * inv(sigma_matrix) * (dx - mean)'))
    # return gaussian[1, 1]
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

@test big(2.0)^big(3) == 8

for T in [Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128, BigInt]
    @test T(2)^big(3.0) == 8
    @test big(2.0)^T(3) == 8
end

for x in (2f0, pi, 7.8, big(e))
    @test big(typeof(x)) == typeof(big(x))
    @test big(typeof(complex(x, x))) == typeof(big(complex(x, x)))
end

# issue 15659
@test (setprecision(53) do; big(1/3); end) < 1//3

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
const _allowslow = Ref(true)

allowslow(flag = true) = (_allowslow[] = flag)

function assertslow(op = "Operation")
    _allowslow[] || error("$op is disabled")
    return
end

Base.IndexStyle(::Type{<:GPUArray}) = IndexLinear()

function _getindex(xs::GPUArray{T}, i::Integer) where T
    x = Array{T}(1)
    copy!(x, 1, xs, i, 1)
    return x[1]
end

function Base.getindex(xs::GPUArray{T}, i::Integer) where T
    assertslow("getindex")
    _getindex(xs, i)
end

function _setindex!(xs::GPUArray{T}, v::T, i::Integer) where T
    x = T[v]
    copy!(xs, i, x, 1, 1)
    return v
end

function Base.setindex!(xs::GPUArray{T}, v::T, i::Integer) where T
    assertslow("setindex!")
    _setindex!(xs, v, i)
end

Base.setindex!(xs::GPUArray, v, i::Integer) = xs[i] = convert(eltype(xs), v)

# Vector indexing

using Base.Cartesian
to_index(a, x) = x
to_index(::A, x::Array{ET}) where {A, ET} = copy!(similar(A, ET, size(x)), x)
to_index(a, x::UnitRange{<: Integer}) = convert(UnitRange{UInt32}, x)
to_index(a, x::Base.LogicalIndex) = error("Logical indexing not implemented")

@generated function index_kernel(state, dest::AbstractArray, src::AbstractArray, idims, Is)
    N = length(Is.parameters)
    quote
        i = linear_index(state)
        i > length(dest) && return
        is = gpu_ind2sub(idims, i)
        @nexprs $N i -> @inbounds I_i = Is[i][Int(is[i])]
        @inbounds dest[i] = @ncall $N getindex src i -> I_i
        return
    end
end

function Base._unsafe_getindex!(dest::GPUArray, src::GPUArray, Is::Union{Real, AbstractArray}...)
    if length(Is) == 1 && isa(first(Is), Array) && isempty(first(Is)) # indexing with empty array
        return dest
    end
    idims = map(length, Is)
    gpu_call(index_kernel, dest, (dest, src, UInt32.(idims), map(x-> to_index(dest, x), Is)))
    return dest
end


text:
#!/bin/env julia

# The MIT License (MIT)
#
# Copyright (c) 2021 Alexander Samoilov
#
# 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.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
# FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
# COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
# IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

#using DocOpt
using Printf
using DataFrames

macro p_str(s) s end
const flt_regexp = p"\-?\d+(?:\.\d*)?(?:[eE][+\-]?\d+)?"

mutable struct TimingsSummary
    matrix_name::AbstractString
    rows::Int
    solver_type::AbstractString
    processor::AbstractString
    setup_time::Float64
    solve_time::Float64
    iterations::Int
    error::Float64

    # a constructor
    TimingsSummary(s::AbstractString = "", m::AbstractString = "") =
        new(m, 0, s, "CPU", 0, 0, 0, 0)
end

@enum State none test_name rows iterations error setup_time solve_time

function parse_case(case_params::Vector{Tuple{String, String}}, precision_data::AbstractString, case_data::AbstractString)
    #println("case_data: ", case_data)
    #println("precision_data: ", precision_data)
    precision = match(r"(float|double)", precision_data)
    d = split(case_data)
    matrix_name = basename(d[1])
    matrix_name_elems = split(first(split(matrix_name, ".")),"_")
    matrix_size = last(matrix_name_elems)
    polynomial_degree = matrix_name_elems[3]
    push!(case_params, ("precision", precision[1]))
    push!(case_params, ("matrix_name", matrix_name))
    push!(case_params, ("matrix_size", matrix_size))
    push!(case_params, ("polynomial_degree", polynomial_degree))
    params = [Tuple(String(x) for x in split(z, "=")) for z in d if occursin("=", z)]
    for z in params
        push!(case_params, z)
    end
end

function main()
    if length(ARGS) < 1
        println("gimme input")
        exit(0)
    end

    ifname = ARGS[1]
    state = none
    summary = TimingsSummary[]
    dfs = DataFrame()
    open(ifname, "r") do f
        timing_summary = TimingsSummary()
        case_params = Vector{Tuple{String, String}}()
        for line in eachline(f)
            ##println(line)
            m = match(r"(.+)(solver.*)\s+\-A\s+(.+)", line)
            if m != nothing
                parse_case(case_params, m[1], m[3])
                insert!(case_params, 3, ("solver_program", m[2]))
                #@printf("solver: %s matrix: %s\n", m[2], m[3])
                #println("case_params: ", case_params)
                timing_summary = TimingsSummary(m[2], m[3])
                state = test_name
                continue
            end
            m = match(r"rows:\s*(\d+)", line)
            if m != nothing
                timing_summary.rows = parse(Int, m[1])
                push!(case_params, ("rows", m[1]))
                #@printf("error: %d\n", timing_summary.rows)
                state = rows
                continue
            end
            m = match(r"Iterations:\s*(\d+)", line)
            if m != nothing
                timing_summary.iterations = parse(Int, m[1])
                push!(case_params, ("iterations", m[1]))
                #@printf("error: %d\n", timing_summary.iterations)
                state = iterations
                continue
            end
            m = match(Regex("Error:\\s*($flt_regexp)"), line)
            if m != nothing
                timing_summary.error = parse(Float64, m[1])
                push!(case_params, ("error", m[1]))
                #@printf("error: %g\n", timing_summary.error)
                state = error
                continue
            end
            m = match(Regex("\\[\\s*setup:\\s*($flt_regexp)\\s*s\\](.+)"), line)
            if m != nothing
                timing_summary.setup_time = parse(Float64, m[1])
                push!(case_params, ("c.setup_time", m[1]))
                #insert!(case_params, 4, ("setup_time", m[1]))
                #@printf("setup_time: %g\n", timing_summary.setup_time)
                state = setup_time
                continue
            end
            m = match(Regex("\\[\\s*solve:\\s*($flt_regexp)\\s*s\\](.+)"), line)
            if m != nothing
                timing_summary.solve_time = parse(Float64, m[1])
                push!(case_params, ("c.solve_time", m[1]))
                #insert!(case_params, 3, ("solve_time", m[1]))
                #@printf("solve_time: %g\n", timing_summary.solve_time)
                state = solve_time
                continue
            end

            if state == solve_time # finished parsing the case
                #println("entry: $timing_summary")
                push!(summary, timing_summary)
                #println("case_params: ", case_params)
                df = DataFrame(Dict(case_params))
                #println(df)
                if ncol(df) >= 15
                    dfs = vcat(dfs, df)
                end
                case_params = Vector{Tuple{String, String}}()
                state = none
            end
        end
    end

    #println(describe(dfs))
    #println(dfs)

    speedup_dict = Dict{AbstractString, Vector{DataFrameRow}}()

    dfs_good = filter(row -> parse(Float64, row.error) < 1e-3 && parse(Int, row.iterations) < 100, dfs)
    println(dfs_good)

    for r in eachrow(dfs_good)
        #println(r)
        local m_name = r.matrix_name
        #println(matrix_name)
        if !haskey(speedup_dict, m_name)
            speedup_dict[m_name] = [r]
        else
            push!(speedup_dict[r.matrix_name], r)
        end
    end

    for (k, v) in speedup_dict
        println(k, " -> ", length(v))
        sort!(v, lt = (x,y) -> x."c.solve_time" < y."c.solve_time")
        println(NamedTuple(v[1]))
        println(NamedTuple(v[2]))
        println(NamedTuple(v[3]))
        println("===================================================================================================================================================================================\n\n")
    end

    sorted = sort!(summary, lt = (a,b) -> a.rows < b.rows)

    #println("sorted: $sorted")

    indices = Dict{AbstractString, Vector{UInt}}()
    for (index, entry) in enumerate(sorted)
        #println("$index: $entry")
        local test_name = entry.matrix_name
        if !haskey(indices, test_name)
            # a new one
            indices[test_name] = [index]
        else
            push!(indices[test_name], index)
        end
    end

    for (k, v) in indices
        #println("$k : $v")
        local (cpu_result_idx, gpu_result_idx) = (0, 0)
        for index in v
            sum = summary[index]
            if sum.solver_type == "solver_cuda"
                gpu_result_idx = index
            else
                cpu_result_idx = index
            end
        end
        if gpu_result_idx > 0 && cpu_result_idx > 0
            local (cpu_info, gpu_info) = (summary[cpu_result_idx], summary[gpu_result_idx])
            #@printf("%s: CPU solve_time: %g GPU solve_time: %g speedup: %g CPU error: %g GPU error: %g\n",
            #        k, cpu_info.solve_time, gpu_info.solve_time,
            #        cpu_info.solve_time / gpu_info.solve_time, cpu_info.error, gpu_info.error)
        end
    end

end

@time main()
#@time main()

max_stars_repo_name: alsam/scriptics-samples

max_issues_repo_name: alsam/scriptics-samples

max_forks_repo_name: alsam/scriptics-samples


text:
# core functions for evaluating tests

"""
    bayestest(x, y)

Compares two variants `x` and `y` and prints summary results

- Assumes uniform prior
- auto-detects distribution type
- return a BayesTestEval summary object
"""
bayestest(x, y)

"""
    bayestest(x, y, p)

Compares two variants `x` and `y` using prior of `p` and prints summary results

- auto-detects distribution type
- return a BayesTestEval summary object
"""
bayestest(x, y, p)

"""
    bayestest(x, y, p, d)

Compares two variants `x` and `y` using prior of `p` and distribution of `d`
and prints summary results.

- return a BayesTestEval summary object
"""
bayestest(x, y, p, d)


text:
# ---
# title: 1296. Divide Array in Sets of K Consecutive Numbers
# id: problem1296
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Greedy
# link: <https://leetcode.com/problems/divide-array-in-sets-of-k-consecutive-numbers/description/>
# hidden: true
# ---
#
# Given an array of integers `nums` and a positive integer `k`, find whether
# it's possible to divide this array into sets of `k` consecutive numbers
# Return `True` if its possible ** ** otherwise return `False`.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,3,3,4,4,5,6], k = 4
#     Output: true
#     Explanation: Array can be divided into [1,2,3,4] and [3,4,5,6].
#
#
# **Example 2:**
#
#
#
#     Input: nums = [3,2,1,2,3,4,3,4,5,9,10,11], k = 3
#     Output: true
#     Explanation: Array can be divided into [1,2,3] , [2,3,4] , [3,4,5] and [9,10,11].
#
#
# **Example 3:**
#
#
#
#     Input: nums = [3,3,2,2,1,1], k = 3
#     Output: true
#
#
# **Example 4:**
#
#
#
#     Input: nums = [1,2,3,4], k = 3
#     Output: false
#     Explanation: Each array should be divided in subarrays of size 3.
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 10^5`
#   * `1 <= nums[i] <= 10^9`
#   * `1 <= k <= nums.length`
#
# **Note:** This question is the same as 846:
# <https://leetcode.comhttps://leetcode.com/problems/hand-of-straights/>
#
#
## @lc code=start
using LeetCode

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


text:
function sequentialSearch(A, t)

  for (index, value) in enumerate(A)
    if value == t
      return index
    end
  end
  return false

end

function main()

  # list to be sorted
  rng = MersenneTwister(1234);
  list_size = 1e7
  A = shuffle(rng, Vector(1:list_size))

  # Search
  search_num = list_size/2
  index_found = @time sequentialSearch(A, search_num)

  print("Search for : ")
  println(search_num)

  print("Found index: ")
  println(index_found)

  print("Found value is ... ")
  println(A[index_found])

end

if contains(@__FILE__, PROGRAM_FILE)
    main()
end

max_stars_repo_name: szmlb/algorithmsInJulia

max_issues_repo_name: szmlb/algorithmsInJulia

max_forks_repo_name: szmlb/algorithmsInJulia


text:
"""
    rand(::Type{MC}, m::Model)

Draw random configuration.
"""
Base.rand(::Type{MC}, m::Model) = error("Model $(typeof(m)) doesn't implement `rand(::Type{MC}, m::Model)`!")

"""
    propose_local(mc::MC, m::Model, i::Int, conf) -> delta_E, delta_i

Propose a local move for lattice site `i` of current configuration `conf`.
Returns local move information `delta_i` (e.g. `new[i] - conf[i]`, will be
forwarded to `accept_local!`) and energy difference `delta_E = E_new - E_old`.

See also [`accept_local!`](@ref).
"""
propose_local(mc::MC, m::Model, i::Int, conf) = error("Model has no implementation of `propose_local(mc::MC, m::Model, i::Int, conf)`!")

"""
    accept_local(mc::MC, m::Model, i::Int, conf, delta_i, delta_E::Float64)

Accept a local move for site `i` of current configuration `conf`.
Arguments `delta_i` and `delta_E` correspond to output of `propose_local()`
for that local move.

See also [`propose_local`](@ref).
"""
accept_local!(mc::MC, m::Model, i::Int, conf, delta_i, delta_E::Float64) = error("Model has no implementation of `accept_local!(m::Model, i::Int, conf, delta_i, delta_E::Float64)`!")


text:
# test initilization functions

for T in (Float64, Float32)

    X = rand(T, 8, 12)

    # randinit

    W, H = NMF.randinit(X, 5)
    @test size(W) == (8, 5)
    @test size(H) == (5, 12)
    @test all(W .>= zero(T))
    @test all(H .>= zero(T))

    W, H = NMF.randinit(X, 5; zeroh=true)
    @test size(W) == (8, 5)
    @test size(H) == (5, 12)
    @test all(W .>= zero(T))
    @test all(H .== zero(T))

    W, H = NMF.randinit(X, 5; normalize=true)
    @test size(W) == (8, 5)
    @test size(H) == (5, 12)
    @test all(W .>= zero(T))
    @test all(H .>= zero(T))
    @test vec(sum(W, dims=1)) ≈ ones(5)

    ## nndsvd

    Random.seed!(5678)
    W, H = NMF.nndsvd(X, 5)
    @test size(W) == (8, 5)
    @test size(H) == (5, 12)
    @test all(W .>= zero(T))
    @test all(H .>= zero(T))

    Random.seed!(5678)
    W2, H2 = NMF.nndsvd(X, 5; zeroh=true)
    @test size(W) == (8, 5)
    @test size(H) == (5, 12)
    @test all(W2 .== W)
    @test all(H2 .== zero(T))

    W, H = NMF.nndsvd(X, 5; variant=:ar)
    @test all(W .> zero(T))
    # NMF.printf_mat(W)
end


text:
"""
    planar push block
        particle with contacts at each corner
"""
struct PlanarPushBlock{I, T} <: Model{I, T}
    n::Int
    m::Int
    d::Int

    mass
    J
    μ
    g

    contact_corner_offset
	control_input_offset

    nq
    nu
    nc
    nf
    nb

    idx_u
    idx_λ
    idx_b
    idx_ψ
    idx_η
    idx_s
end

# Dimensions
nq = 3 # configuration dimension
nu = 2 * 4 + 4 # control dimension
nc = 4 # number of contact points
nf = 4 # number of faces for friction cone pyramid
nb = nc * nf
ns = 1

n = 2 * nq
m = nu + nc + nb + nc + nb + 1
d = 0

idx_u = (1:nu)
idx_λ = nu .+ (1:nc)
idx_b = nu + nc .+ (1:nb)
idx_ψ = nu + nc + nb .+ (1:nc)
idx_η = nu + nc + nb + nc .+ (1:nb)
idx_s = nu + nc + nb + nc + nb .+ (1:ns)

# Kinematics
r = 0.1

# contact corner
cc1 = @SVector [r, r]
cc2 = @SVector [-r, r]
cc3 = @SVector [r, -r]
cc4 = @SVector [-r, -r]

# control location
cu5 = @SVector [r, 0.0]
cu6 = @SVector [0.0, r]
cu7 = @SVector [-r, 0.0]
cu8 = @SVector [0.0, -r]

contact_corner_offset = @SVector [cc1, cc2, cc3, cc4]
control_input_offset = @SVector [cu5, cu6, cu7, cu8]

# Parameters
μ = 1.0  # coefficient of friction
g = 9.81
mass = 1.0   # mass
J = 1.0 / 12.0 * mass * ((2.0 * r)^2 + (2.0 * r)^2)

# Methods
M_func(model::PlanarPushBlock, q) = Diagonal(@SVector [model.mass, model.mass, model.J])

function C_func(model::PlanarPushBlock, q, q̇)
	SVector{3}([0.0, 0.0, 0.0])
end

function rotation_matrix(x)
	SMatrix{2,2}([cos(x) -sin(x); sin(x) cos(x)])
end

function ϕ_func(model::PlanarPushBlock, q)
    p = view(q, 1:2)
	θ = q[3]
    R = rotation_matrix(θ)

    @SVector [0.0, 0.0, 0.0, 0.0]
end

function control_kinematics_func(model::PlanarPushBlock, q)
    p = q[1:2]
	θ = q[3]
    R = rotation_matrix(θ)

    SVector{8}([p + R * model.control_input_offset[1];
              p + R * model.control_input_offset[2];
              p + R * model.control_input_offset[3];
              p + R * model.control_input_offset[4]])
              # p + R * model.control_input_offset[5];
              # p + R * model.control_input_offset[6];
              # p + R * model.control_input_offset[7];
              # p + R * model.control_input_offset[8]])
end

function B_func(model::PlanarPushBlock, q, u)
	tmp(z) = control_kinematics_func(model, z)
	ForwardDiff.jacobian(tmp, q)

	mask = Diagonal([1, 0, 1, 0, 1, 0, 1, 0])

	R1 = rotation_matrix(-1.0 * π)
	R2 = rotation_matrix(-0.5 * π)
	R3 = rotation_matrix(0.0)
	R4 = rotation_matrix(0.5 * π)
	R = cat(R1, R2, R3, R4, dims=(1, 2))

	T1 = rotation_matrix(u[2])
	T2 = rotation_matrix(u[4])
	T3 = rotation_matrix(u[6])
	T4 = rotation_matrix(u[8])
	T = cat(T1, T2, T3, T4, dims=(1, 2))

	transpose(ForwardDiff.jacobian(tmp, q)) * T * R * mask * u
end

function N_func(model::PlanarPushBlock, q)
    tmp(z) = ϕ_func(model, z)
    ForwardDiff.jacobian(tmp, q)
end

function P_func(model::PlanarPushBlock, q)
    map = [1. 0.;
           0. 1.;
           -1. 0.;
           0. -1.]

    function p(x)
        pos = view(x, 1:2)
		θ = x[3]
        R = rotation_matrix(θ)

        [map * (pos + R * model.contact_corner_offset[1])[1:2];
         map * (pos + R * model.contact_corner_offset[2])[1:2];
         map * (pos + R * model.contact_corner_offset[3])[1:2];
         map * (pos + R * model.contact_corner_offset[4])[1:2]]
    end

    ForwardDiff.jacobian(p, q)
end

function friction_cone(model::PlanarPushBlock,u)
    λ = u[model.idx_λ]
    b = u[model.idx_b]

    @SVector [model.μ[1] * λ[1] - sum(b[1:4]),
              model.μ[2] * λ[2] - sum(b[5:8]),
              model.μ[3] * λ[3] - sum(b[9:12]),
              model.μ[4] * λ[4] - sum(b[13:16])]
end

function maximum_dissipation(model::PlanarPushBlock, x⁺, u, h)
    q3 = x⁺[model.nq .+ (1:model.nq)]
	q2 = x⁺[1:model.nq]

    ψ = u[model.idx_ψ]
    ψ_stack = [ψ[1] * ones(4);
               ψ[2] * ones(4);
               ψ[3] * ones(4);
               ψ[4] * ones(4)]

    η = u[model.idx_η]

    P_func(model, q3) * (q3 - q2) / h + ψ_stack - η
end

function lagrangian_derivatives(model, q, v)
	D1L = -1.0 * C_func(model, q, v)
    D2L = M_func(model, q) * v
	return D1L, D2L
end

function fd(model::PlanarPushBlock{Discrete, FixedTime}, x⁺, x, u, w, h, t)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2⁺ = view(x⁺, 1:model.nq)
	q2⁻ = view(x, model.nq .+ (1:model.nq))
	q1 = view(x, 1:model.nq)
	u_ctrl = view(u, model.idx_u)
	λ = view(u, model.idx_λ)
	b = view(u, model.idx_b)

	qm1 = 0.5 * (q1 + q2⁺)
    vm1 = (q2⁺ - q1) / h
    qm2 = 0.5 * (q2⁺ + q3)
    vm2 = (q3 - q2⁺) / h

	D1L1, D2L1 = lagrangian_derivatives(model, qm1, vm1)
	D1L2, D2L2 = lagrangian_derivatives(model, qm2, vm2)

	[q2⁺ - q2⁻;
     (0.5 * h * D1L1 + D2L1 + 0.5 * h * D1L2 - D2L2
     + B_func(model, qm2, u_ctrl[1:8])
     + transpose(P_func(model, q3)) * SVector{16}(b))]
end

model = PlanarPushBlock{Discrete, FixedTime}(n, m, d,
			mass, J, [[μ for i = 1:nc]..., 0.1], g,
			contact_corner_offset,
			control_input_offset,
            nq, nu, nc, nf, nb,
            idx_u,
            idx_λ,
            idx_b,
            idx_ψ,
            idx_η,
            idx_s)

# function control_input(q, u)
# 	idx = [(i - 1) * 2 .+ (1:2) for i = 1:4]
#
#     k = control_kinematics_func(model, q)
# 	k_input = u[model.idx_u][9:10]
#
# 	d1 = norm(k[1:2] - k_input)
# 	d2 = norm(k[3:4] - k_input)
# 	d3 = norm(k[5:6] - k_input)
# 	d4 = norm(k[7:8] - k_input)
# 	# d5 = norm(k[9:10] - k_input)
# 	# d6 = norm(k[11:12] - k_input)
# 	# d7 = norm(k[13:14] - k_input)
# 	# d8 = norm(k[15:16] - k_input)
#
# 	_, min_idx = findmin([d1, d2, d3, d4])#, d5, d6, d7, d8])
#
# 	return u[idx[min_idx]], u[9:10], min_idx
# end

function visualize!(vis, model::PlanarPushBlock, q, u; r = r,
        Δt = 0.1)

	default_background!(vis)

    setobject!(vis["box"], GeometryBasics.Rect(Vec(-1.0 * r,
		-1.0 * r,
		-1.0 * r),
		Vec(2.0 * r, 2.0 * r, 2.0 * r)),
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

    for i = 1:model.nc
        setobject!(vis["contact$i"], GeometryBasics.Sphere(Point3f0(0),
            convert(Float32, 0.02)),
            MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0.0, 1.0)))
    end

	# for i = 1:4
    #     setobject!(vis["control$i"], GeometryBasics.Sphere(Point3f0(0),
    #         convert(Float32, 0.02)),
    #         MeshPhongMaterial(color = RGBA(51.0 / 255.0, 1.0, 1.0, 1.0)))
    # end

	# force_vis = ArrowVisualizer(vis[:force])
	# setobject!(force_vis, MeshPhongMaterial(color=RGBA(1.0, 0.0, 0.0, 1.0)))
	#
	# uf, up, min_idx = control_input(q[2], u[1])
	# uf_norm = uf / 10.0
	# settransform!(force_vis,
	# 			Point(up[1] - uf_norm[1], up[2] - uf_norm[2], 2 * r),
	# 			Vec(uf_norm[1], uf_norm[2], 2 * r),
	# 			shaft_radius=0.01,
	# 			max_head_radius=0.025)

    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

	T = length(q)
    for t = 1:T-1
        MeshCat.atframe(anim, t) do
			if t < T-1
				# uf, up, min_idx = control_input(q[t+1], u[t])
				# println("t = $t, control_input: $min_idx")
				#
				# if norm(uf) < 1.0e-6
				# 	setvisible!(vis[:force], false)
				# else
				# 	setvisible!(vis[:force], true)
				# 	uf_norm = uf ./ 10.0
				# 	settransform!(force_vis,
				# 				Point(up[1] - uf_norm[1], up[2] - uf_norm[2], 2 * r),
				# 				Vec(uf_norm[1], uf_norm[2], 2 * r),
				# 				shaft_radius=0.01,
				# 				max_head_radius=0.025)
				# end
			end

            settransform!(vis["box"],
				compose(Translation(q[t+1][1], q[t+1][2], r), LinearMap(RotZ(q[t+1][3]))))

            for i = 1:model.nc
                settransform!(vis["contact$i"],
                    Translation(([q[t+1][1:2]; 0.0] + RotZ(q[t+1][3]) * [contact_corner_offset[i]; 0.0])...))
            end

			# for i = 1:4
			# 	settransform!(vis["control$i"],
			# 		Translation(([q[t+1][1:2]; r] + RotZ(q[t+1][3]) * [control_input_offset[i]; 0.0])...))
			# end
        end
    end

	settransform!(vis["/Cameras/default"],
		compose(Translation(0.0, 0.0, 50.0), LinearMap(RotZ(0.5 * pi) * RotY(-pi/2.5))))
	setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 25)

    MeshCat.setanimation!(vis, anim)
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
abstract type AbstractWorld{T, ND} end

struct SSDInterval{T <: SSDFloat, L, R, BL, BR} <: IntervalSets.TypedEndpointsInterval{L,R,T}
    left::T
    right::T
end

function get_boundary_types(int::SSDInterval{T, L, R, BL, BR}) where {T <:SSDFloat, L, R, BL, BR}
    return L, R, BL, BR
end

const boundary_condition_mapping = Dict{String, Symbol}(
    "r0" => :r0,
    "inf" => :infinite,
    "infinite" => :infinite,
    "fixed" => :fixed,
    "fix" => :fixed,
    "reflecting" => :reflecting,
    "periodic" => :periodic,
)

struct World{T <: SSDFloat, ND, S} <: AbstractWorld{T, ND}
    intervals::NTuple{ND, SSDInterval{T}}
end

function World{T, ND, S}(args...) where {T <: SSDFloat, ND, S}
    return World{T, ND, S}(args)
end

function World(T, dict::Dict, inputunit_dict::Dict{String, Unitful.Units})::World
    if dict["coordinates"] == "cylindrical"
        CylindricalWorld(T, dict["axes"], inputunit_dict)
    elseif dict["coordinates"] == "cartesian"
        CartesianWorld(T, dict["axes"], inputunit_dict)
    else
        error("Gridtype must be \"cylindrical\" or \"cartesian\"")
    end
end

function get_periodicity(int::SSDInterval{T, L, R, BL, BR}) where {T <:SSDFloat, L, R, BL, BR}
    return int.right - int.left
end

function get_interval_boundary_types(dict::Dict)
    BL, BR = missing, missing
    if haskey(dict, "boundaries")
        if typeof(dict["boundaries"]) == String
            BL = boundary_condition_mapping[dict["boundaries"]]
            BR = boundary_condition_mapping[dict["boundaries"]]
        else
            BL = boundary_condition_mapping[dict["boundaries"]["left"]]
            BR = boundary_condition_mapping[dict["boundaries"]["right"]]
        end
    end
    if !ismissing(BL)
        if BL == :periodic || BR == :periodic
            if !(BL == :periodic && BR == :periodic) && !((BL == :periodic && BR == :reflecting) || (BL == :reflecting && BR == :periodic))
                throw(ConfigFileError("both or none endings must be \"periodic\" for an interval. Or, \"periodic\" and \"reflecting\" (for periodic behaviour plus mirror symmetry)."))
            end
        end
    end
    return BL, BR
end

function is_periodic_plus_mirror_symmetric(BL::Symbol, BR::Symbol)::Bool
    return (BL == :periodic && BR == :reflecting) || (BL == :reflecting && BR == :periodic)
end

function get_r_SSDInterval(T, dict, inputunit_dict::Dict{String, Unitful.Units})
    if isa(dict, Dict)
        from::T = 0
        if "from" in keys(dict) @warn "ConfigFileWarning: \"from\" is not used in r-axis. It is fixed to 0." end
        to::T = "to" in keys(dict) ? geom_round(ustrip(uconvert(internal_length_unit, T(dict["to"]) * inputunit_dict["length"]))) : throw(ConfigFileError("No \"to\" given for r-axis."))
        if from < 0 throw(ConfigFileError("left boundary of r-axis cannot be negative.")) end
        if to < 0 throw(ConfigFileError("right boundary of r-axis cannot be negative.")) end
        L = :closed
        BL = :r0
        R, BR = :closed, :infinite
        if "boundaries" in keys(dict)
            BR = boundary_condition_mapping[dict["boundaries"]]
        end
        return SSDInterval{T, L, R, BL, BR}(from, to)
    else
        SSDInterval{T, :closed, :closed, :r0, :infinite}(0, geom_round(ustrip(uconvert(internal_length_unit, T(dict) * inputunit_dict["length"]))))
    end
end
function get_φ_SSDInterval(T, dict::Dict, inputunit_dict::Dict{String, Unitful.Units})
    if haskey(dict, "phi")
        dp = dict["phi"]
        from::T = "from" in keys(dp) ? geom_round(ustrip(uconvert(internal_angle_unit, T(dp["from"]) * inputunit_dict["angle"]))) : T(0)
        to::T = "to" in keys(dp) ? geom_round(ustrip(uconvert(internal_angle_unit, T(dp["to"]) * inputunit_dict["angle"]))) : T(2π)
        L = :closed
        R = :open
        BL = :periodic
        BR = :periodic
        cfBL, cfBR = get_interval_boundary_types(dp)
        if !ismissing(cfBL) BL = cfBL; BR = cfBR; end
        if is_periodic_plus_mirror_symmetric(BL, BR)
            L = :closed; R = :closed;
            BL = :reflecting; BR = :reflecting
        end
        if from == to # 2D
            L = :closed; R = :closed;
            BL = :reflecting; BR = :reflecting
        end
        return SSDInterval{T, L, R, BL, BR}(from, to)
    else
        return SSDInterval{T, :closed, :open, :periodic, :periodic}(T(0), T(2π))
    end
end

function get_cartesian_SSDInterval(T, dict::Dict, inputunit_dict::Dict{String, Unitful.Units})
    from::T = "from" in keys(dict) ? geom_round(ustrip(uconvert(internal_length_unit, T(dict["from"]) * inputunit_dict["length"]))) : 0
    to = "to" in keys(dict) ? geom_round(ustrip(uconvert(internal_length_unit, T(dict["to"]) * inputunit_dict["length"]))) : throw(ConfigFileError("No \"to\" given for z-axis."))
    L = :closed
    R = :closed
    BL, BR = :infinite, :infinite
    cfBL, cfBR = get_interval_boundary_types(dict)
    if !ismissing(cfBL) BL = cfBL; BR = cfBR; end
    if BL == BR == :periodic
        L = :closed; R = :open
    end
    return SSDInterval{T, L, R, BL, BR}(from, to)
end

function CylindricalWorld(T, dict::Dict, inputunit_dict::Dict{String, Unitful.Units})::World
    r_int = get_r_SSDInterval(T, dict["r"], inputunit_dict)
    φ_int = get_φ_SSDInterval(T, dict, inputunit_dict)
    z_int = get_cartesian_SSDInterval(T, dict["z"], inputunit_dict)
    return World{T, 3, :cylindrical}( r_int, φ_int, z_int )
end

function CartesianWorld(T, dict::Dict, inputunit_dict::Dict{String, Unitful.Units})::World
    x_int = get_cartesian_SSDInterval(T, dict["x"], inputunit_dict)
    y_int = get_cartesian_SSDInterval(T, dict["y"], inputunit_dict)
    z_int = get_cartesian_SSDInterval(T, dict["z"], inputunit_dict)
    return World{T, 3, :cartesian}( x_int, y_int, z_int )
end

function CartesianWorld(xl::T, xr::T, yl::T, yr::T, zl::T, zr::T)::World where {T <: SSDFloat}
    Δx::T = (xr - xl) / 10
    Δy::T = (yr - yl) / 10
    Δz::T = (zr - zl) / 10
    x_int = SSDInterval{T, :closed, :closed, :infinite, :infinite}(xl - Δx, xr + Δx)
    y_int = SSDInterval{T, :closed, :closed, :infinite, :infinite}(yl - Δy, yr + Δy)
    z_int = SSDInterval{T, :closed, :closed, :infinite, :infinite}(zl - Δz, zr + Δz)
    return World{T, 3, :cartesian}( x_int, y_int, z_int )
end

function CylindricalWorld(r_max::T, zl::T, zr::T)::World where {T <: SSDFloat}
    r_int = SSDInterval{T, :closed, :closed, :r0, :infinite}(T(0), abs(r_max * T(1.1)))
    φ_int = SSDInterval{T, :closed, :open, :periodic, :periodic}(T(0), T(2π))
    Δz::T = zr - zl
    z_int = SSDInterval{T, :closed, :closed, :infinite, :infinite}(zl - (Δz / 10), zr + (Δz / 10))
    return World{T, 3, :cylindrical}( r_int, φ_int, z_int )
end

function World(S::Val{:cylindrical}, limits::NTuple{6, T})::World where {T <: SSDFloat}
    return CylindricalWorld(limits[2], limits[5], limits[6])
end
function World(S::Val{:cartesian}, limits::NTuple{6, T})::World where {T <: SSDFloat}
    return CartesianWorld(limits...)
end


text:
using GriddingMachine.Collections
using GriddingMachine.Indexer
using GriddingMachine.Requestor
using Test

# test Collections functions
println();
@testset "GriddingMachine : Library" begin
    canopy_height_collection();    @test true;
    clumping_index_collection();   @test true;
    elevation_collection();        @test true;
    gpp_collection();              @test true;
    lai_collection();              @test true;
    land_mask_collection();        @test true;
    leaf_chlorophyll_collection(); @test true;
    leaf_nitrogen_collection();    @test true;
    leaf_phosphorus_collection();  @test true;
    pft_collection();              @test true;
    sif_collection();              @test true;
    sla_collection();              @test true;
    soil_color_collection();       @test true;
    soil_hydraulics_collection();  @test true;
    surface_area_collection();     @test true;
    tree_density_collection();     @test true;
    vcmax_collection();            @test true;
    wood_density_collection();     @test true;
end;

println();
@testset "GriddingMachine : Query Collections" begin
    # test query_collection function
    query_collection(pft_collection()); @test true;
    query_collection(sla_collection()); @test true;
    query_collection("PFT_2X_1Y_V1"); @test true;

    # clean up artifacts
    clean_collections!("old"); @test true;
    clean_collections!(pft_collection()); @test true;

    # only for high memory and storage cases, e.g., server
    if Sys.islinux() && (Sys.total_memory() / 2^30) > 30
        query_collection(canopy_height_collection(), "20X_1Y_V1"); @test true;
        query_collection(canopy_height_collection(), "2X_1Y_V2" ); @test true;

        query_collection(elevation_collection(), "4X_1Y_V1" ); @test true;

        query_collection(clumping_index_collection(), "240X_1Y_V1"); @test true;
        query_collection(clumping_index_collection(), "2X_1Y_V1"  ); @test true;
        query_collection(clumping_index_collection(), "2X_1Y_V2"  ); @test true;

        for year in 2001:2019
            query_collection(gpp_collection(), "MPI_RS_2X_1M_$(year)_V1"); @test true;
            query_collection(gpp_collection(), "MPI_RS_2X_8D_$(year)_V1"); @test true;
        end;
        for year in 2000:2019
            query_collection(gpp_collection(), "VPM_5X_8D_$(year)_V2" ); @test true;
            query_collection(gpp_collection(), "VPM_12X_8D_$(year)_V2"); @test true;
        end;

        for year in 2000:2020
            query_collection(lai_collection(), "MODIS_2X_1M_$(year)_V1" ); @test true;
            query_collection(lai_collection(), "MODIS_2X_8D_$(year)_V1" ); @test true;
            query_collection(lai_collection(), "MODIS_10X_1M_$(year)_V1"); @test true;
            query_collection(lai_collection(), "MODIS_10X_8D_$(year)_V1"); @test true;
            query_collection(lai_collection(), "MODIS_20X_1M_$(year)_V1"); @test true;
            query_collection(lai_collection(), "MODIS_20X_8D_$(year)_V1"); @test true;
        end;

        query_collection(land_mask_collection(), "4X_1Y_V1"); @test true;

        # expect warning here
        query_collection(leaf_chlorophyll_collection(), "2X_7D_V1"); @test true;

        query_collection(leaf_nitrogen_collection(), "2X_1Y_V1"); @test true;
        query_collection(leaf_nitrogen_collection(), "2X_1Y_V2"); @test true;

        query_collection(leaf_phosphorus_collection(), "2X_1Y_V1"); @test true;

        query_collection(pft_collection(), "2X_1Y_V1"); @test true;

        for year in 2018:2019
            query_collection(sif_collection(), "TROPOMI_740_1X_1M_$(year)_V1"   ); @test true;
            query_collection(sif_collection(), "TROPOMI_740_12X_8D_$(year)_V1"  ); @test true;
            query_collection(sif_collection(), "TROPOMI_740DC_1X_1M_$(year)_V1" ); @test true;
            query_collection(sif_collection(), "TROPOMI_740DC_12X_8D_$(year)_V1"); @test true;
        end;

        query_collection(sla_collection(), "2X_1Y_V1"); @test true;
        query_collection(sla_collection(), "2X_1Y_V2"); @test true;

        query_collection(soil_color_collection(), "2X_1Y_V1"); @test true;

        query_collection(soil_hydraulics_collection(), "SWCR_120X_1Y_V1"); @test true;
        query_collection(soil_hydraulics_collection(), "SWCR_12X_1Y_V1" ); @test true;
        query_collection(soil_hydraulics_collection(), "SWCS_120X_1Y_V1"); @test true;
        query_collection(soil_hydraulics_collection(), "SWCS_12X_1Y_V1" ); @test true;
        query_collection(soil_hydraulics_collection(), "VGA_120X_1Y_V1" ); @test true;
        query_collection(soil_hydraulics_collection(), "VGA_12X_1Y_V1"  ); @test true;
        query_collection(soil_hydraulics_collection(), "VGN_120X_1Y_V1" ); @test true;
        query_collection(soil_hydraulics_collection(), "VGN_12X_1Y_V1"  ); @test true;

        query_collection(surface_area_collection(), "2X_1Y_V1"); @test true;
        query_collection(surface_area_collection(), "1X_1Y_V1"); @test true;

        query_collection(tree_density_collection(), "120X_1Y_V1"); @test true;
        query_collection(tree_density_collection(), "2X_1Y_V1"  ); @test true;

        query_collection(vcmax_collection(), "2X_1Y_V1"); @test true;
        query_collection(vcmax_collection(), "2X_1Y_V2"); @test true;

        query_collection(wood_density_collection(), "2X_1Y_V1"); @test true;
    end;
end;

# test Indexer functions
println();
@testset "GriddingMachine : Indexer" begin
    # read the full dataset
    read_LUT(query_collection(vcmax_collection())); @test true;

    # read the global map at a given cycle index
    read_LUT(query_collection(gpp_collection()), 8); @test true;

    # read the data at given lat and lon
    read_LUT(query_collection(vcmax_collection()), 30, 116); @test true;
    read_LUT(query_collection(vcmax_collection()), 30, 116, 0.5); @test true;
    read_LUT(query_collection(vcmax_collection()), 30, 116; interpolation=true); @test true;
    read_LUT(query_collection(vcmax_collection()), 30, 116, 0.5; interpolation=true); @test true;

    # read the data at given lat, lon, and cycle index
    read_LUT(query_collection(gpp_collection()), 30, 116, 8); @test true;
    read_LUT(query_collection(gpp_collection()), 30, 116, 8, 0.5); @test true;
    read_LUT(query_collection(gpp_collection()), 30, 116, 8; interpolation=true); @test true;
    read_LUT(query_collection(gpp_collection()), 30, 116, 8, 0.5; interpolation=true); @test true;
end;

# test Requestor functions
println();
@testset "GriddingMachine : Indexer" begin
    request_LUT("LAI_MODIS_2X_8D_2017_V1", 30.5, 115.5); @test true;
    request_LUT("LAI_MODIS_2X_8D_2017_V1", 30.5, 115.5; interpolation=true); @test true;
    request_LUT("LAI_MODIS_2X_8D_2017_V1", 30.5, 115.5, 8); @test true;
    request_LUT("LAI_MODIS_2X_8D_2017_V1", 30.5, 115.5, 8; interpolation=true); @test true;
end;


text:
using HypothesisTests
using Distances
using DelayEmbeddings

function normalise_minmax(x, vmin, vmax)
    (x - vmin)/(vmax - vmin)
end

"""
    jdd(source, target; distance_metric = SqEuclidean(),
        B::Int = 10, D::Int = 2, τ::Int = 1) → Vector{Float64}

Compute the joint distance distribution [^1] from `source` to `target` using
the provided `distance_metric`, with `B` controlling the number of subintervals,
`D` the embedding dimension and `τ` the embedding lag.

## Example

```julia
using CausalityTools
x, y = rand(1000), rand(1000)

jdd(x, y)
```

## Keyword arguments

- **`distance_metric::Metric`**: An instance of a valid distance metric from `Distances.jl`.
    Defaults to `SqEuclidean()`.
- **`B`**: The number of equidistant subintervals to divide the interval `[0, 1]` into
    when comparing the normalised distances.
- **`D`**: Embedding dimension.
- **`τ`**: Embedding delay.

## References
[^1] Amigó, José M., and Yoshito Hirata. "Detecting directional couplings from multivariate flows by the joint distance distribution." Chaos: An Interdisciplinary Journal of Nonlinear Science 28.7 (2018): 075302.
"""
function jdd(source, target;
        distance_metric = Euclidean(),
        B::Int = 10,
        D::Int = 2,
        τ::Int = 1)

    js = ([1 for i = 1:D]...,)
    τs = (collect(0:-τ:-(D-1)*τ)...,)
    Ex = DelayEmbeddings.genembed(source, τs, js)
    Ey = DelayEmbeddings.genembed(target, τs, js)
    Mx = DelayEmbeddings.Matrix(Ex)
    My = DelayEmbeddings.Matrix(Ey)

    npts = length(Ex)
    Dx = pairwise(distance_metric, Mx, Mx, dims = 1)
    Dy = pairwise(distance_metric, My, My, dims = 1)

    # Normalise the distances to the interval [0, 1]
    Dx_min = minimum(Dx[Dx .> 0])
    Dy_min = minimum(Dy[Dy .> 0])
    Dx_max = maximum(Dx[Dx .> 0])
    Dy_max = maximum(Dy[Dy .> 0])

    Dx_norm = zeros(Float64, size(Dx))
    Dy_norm = zeros(Float64, size(Dy))

    for i in LinearIndices(Dx[Dx .> 0])
        Dx_norm[i] = normalise_minmax(Dx[i], Dx_min, Dx_max)
    end

    for i in LinearIndices(Dy[Dy .> 0])
        Dy_norm[i] = normalise_minmax(Dy[i], Dy_min, Dy_max)
    end

    mins_δ_yi_yj = fill(2.0, 2*B)

    for (k, b) in enumerate(1:2*B)
        bmin = (b-1)/(2*B)
        bmax = b/(2*B)

        # Find the indices of all pairs (i, j) in Dx whose distances fall inside the interval Ib
        #idxs_Dxs_in_Ib = findall(Dx[])

        # We don't need to store any indices or distances explicitly, but only
        # keep track of whether a smaller distance has has been detected.
        # The maximum possible distance after normalisation is 1.0, so this
        # value can only decrease as we update.
        min_δ_yi_yj = 1.0

        for i = 1:npts
            for j = (i+1):npts
                δ_xi_xj = Dx_norm[i, j]
                if bmin < δ_xi_xj <= bmax
                    δ_yi_yj = Dy_norm[i, j]
                    if δ_yi_yj < min_δ_yi_yj
                        min_δ_yi_yj = δ_yi_yj
                    end
                end
            end
        end
        mins_δ_yi_yj[k] = min_δ_yi_yj
    end

    Δjdd = [mins_δ_yi_yj[B + i] - mins_δ_yi_yj[i] for i in 1:B]

    return Δjdd
end

"""
    jdd(test::OneSampleTTest, source, target;
        distance_metric = SqEuclidean(), B::Int = 10, D::Int = 2, τ::Int = 1,
        μ0 = 0.0) → OneSampleTTest

Perform a one sample t-test to check that the joint distance distribution[^Amigo]
computed from `source` to `target` is biased towards positive values, using the null
hypothesis that the mean of the distribution is `μ0`.

The interpretation of the t-test is that if we can reject the null, then the
joint distance distribution is biased towards positive values, and then there
exists an underlying coupling from `source` to `target`.

## Example

```julia
using CausalityTools, HypothesisTests
x, y = rand(1000), rand(1000)

jdd(OneSampleTTest, x, y)
```

which gives

```julia
One sample t-test
-----------------
Population details:
    parameter of interest:   Mean
    value under h_0:         0.0
    point estimate:          0.06361857324022721
    95% confidence interval: (0.0185, 0.1087)

Test summary:
    outcome with 95% confidence: reject h_0
    two-sided p-value:           0.0082

Details:
    number of observations:   20
    t-statistic:              2.9517208721082873
    degrees of freedom:       19
    empirical standard error: 0.0215530451545668
```

The lower bound of the confidence interval for the mean of the joint
distance distribution is `0.0185` at confidence level `α = 0.05`. The
meaning that the test falsely detected causality from `x` to `y`
between these two random time series.

To get the confidence intervals at confidence level `α`, run

```confinf(jdd(OneSampleTTest, x, y), α)```.

If you just want the p-value at 95% confidence, run

```pvalue(jdd(OneSampleTTest, x, y), tail = :right)```

## Keyword arguments

- **`distance_metric::Metric`**: An instance of a valid distance metric from `Distances.jl`.
    Defaults to `SqEuclidean()`.
- **`B`**: The number of equidistant subintervals to divide the interval `[0, 1]` into
    when comparing the normalised distances.
- **`D`**: Embedding dimension.
- **`τ`**: Embedding delay.
- **`μ0`**: The hypothetical mean value of the joint distance distribution if there
    is no coupling between `x` and `y` (default is `μ0 = 0.0`).

## References
[^Amigo]: Amigó, José M., and Yoshito Hirata. "Detecting directional couplings from multivariate flows by the joint distance distribution." Chaos: An Interdisciplinary Journal of Nonlinear Science 28.7 (2018): 075302.
"""
function jdd(test::Type{OneSampleTTest}, source, target;
        distance_metric = SqEuclidean(),
        B::Int = 10,
        D::Int = 2, τ::Int = 1,
        μ0 = 0.0)

    Δjdd = jdd(source, target,
        distance_metric = distance_metric,
        B = B,
        D = D, τ = τ)

    OneSampleTTest(Δjdd, μ0)
end

export jdd


text:
using ITensors,
      Test

import ITensors: In,Out,Neither

@testset "Index" begin
  @testset "Default Index" begin
    i = Index()
    @test id(i) == 0
    @test dim(i) == 1
    @test dir(i) == Neither
    @test -dir(i) == Neither
    @test plev(i) == 0
    @test tags(i) == TagSet(("",0))
  end
  @testset "Index with dim" begin
    i = Index(2)
    @test id(i) != 0
    @test dim(i) == 2
    @test dir(i) == Out
    @test plev(i) == 0
    @test tags(i) == TagSet(("",0))
  end
  @testset "Index with all args" begin
    i = Index(UInt64(1), 2, In, TagSet(("Link",1)))
    @test id(i) == 1
    @test dim(i) == 2
    @test dir(i) == In
    @test plev(i) == 1
    @test tags(i) == TagSet(("Link",1))
    j = copy(i)
    @test id(j) == 1
    @test dim(j) == 2
    @test dir(j) == In
    @test plev(j) == 1
    @test tags(j) == TagSet(("Link",1))
    @test j == i
  end
  @testset "prime" begin
    i = Index(2)
    @test plev(i) == 0
    i2 = prime(i, 2)
    @test plev(i2) == 2
    i1 = i'
    @test plev(i1) == 1
    i2 = i''
    @test plev(i2) == 2
    i3 = i'''
    @test plev(i3) == 3
    i6 = i^6
    @test plev(i6) == 6
    i0 = noprime(i)
    @test plev(i0) == 0
  end
  @testset "IndexVal" begin
    i = Index(2)
    @test_throws ErrorException IndexVal(i, 4)
    @test_throws ErrorException IndexVal(i, 0)
    @test i(2) == IndexVal(i, 2)
    @test val(IndexVal(i, 1)) == 1
    @test ind(IndexVal(i, 1)) == i
    @test i == IndexVal(i, 2)
    @test IndexVal(i, 2) == i
    @test plev(i(2)') == 1
    @test val(i(2)') == 2
    @test plev(prime(i(2),4)) == 4
    #@test i[:] == [i(1); i(2)]
    @test sprint(show, i(2)) == sprint(show, i)*"=2"

    @test IndexVal() == IndexVal(Index(),1)
  end
  @testset "Iteration" begin
    i = Index(10)
    c = 1
    for n in i
      @test n == c
      c += 1
    end
  end
end