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

using Base.Order: Forward

@test sort([2,3,1]) == [1,2,3]
@test sort([2,3,1], rev=true) == [3,2,1]
@test sort(['z':-1:'a';]) == ['a':'z';]
@test sort(['a':'z';], rev=true) == ['z':-1:'a';]
@test sortperm([2,3,1]) == [3,1,2]
@test sortperm!([1,2,3], [2,3,1]) == [3,1,2]
let s = view([1,2,3,4], 1:3)
    r = sortperm!(s, [2,3,1])
    @test r == [3,1,2]
    @test r === s
end
@test_throws ArgumentError sortperm!(view([1,2,3,4], 1:4), [2,3,1])
@test !issorted([2,3,1])
@test issorted([1,2,3])
@test reverse([2,3,1]) == [1,3,2]
@test select([3,6,30,1,9],3) == 6
@test select([3,6,30,1,9],3:4) == [6,9]
@test selectperm([3,6,30,1,9], 3:4) == [2,5]
@test selectperm!(collect(1:5), [3,6,30,1,9], 3:4) == [2,5]
let a=[1:10;]
    for r in Any[2:4, 1:2, 10:10, 4:2, 2:1, 4:-1:2, 2:-1:1, 10:-1:10, 4:1:3, 1:2:8, 10:-3:1]
        @test select(a, r) == [r;]
        @test selectperm(a, r) == [r;]
        @test select(a, r, rev=true) == (11 .- [r;])
        @test selectperm(a, r, rev=true) == (11 .- [r;])
    end
end
@test sum(randperm(6)) == 21

numTypes = [ Int8,  Int16,  Int32,  Int64,  Int128,
            UInt8, UInt16, UInt32, UInt64, UInt128,
            Float16, Float32, Float64, BigInt, BigFloat]

@test searchsorted([1:10;], 1, by=(x -> x >= 5)) == 1:4
@test searchsorted([1:10;], 10, by=(x -> x >= 5)) == 5:10
@test searchsorted([1:5; 1:5; 1:5], 1, 6, 10, Forward) == 6:6
@test searchsorted(ones(15), 1, 6, 10, Forward) == 6:10

for R in numTypes, T in numTypes
    @test searchsorted(R[1, 1, 2, 2, 3, 3], T(0)) == 1:0
    @test searchsorted(R[1, 1, 2, 2, 3, 3], T(1)) == 1:2
    @test searchsorted(R[1, 1, 2, 2, 3, 3], T(2)) == 3:4
    @test searchsorted(R[1, 1, 2, 2, 3, 3], T(4)) == 7:6
    @test searchsorted(R[1, 1, 2, 2, 3, 3], 2.5) == 5:4

    @test searchsorted(1:3, T(0)) == 1:0
    @test searchsorted(1:3, T(1)) == 1:1
    @test searchsorted(1:3, T(2)) == 2:2
    @test searchsorted(1:3, T(4)) == 4:3

    @test searchsorted(R[1:10;], T(1), by=(x -> x >= 5)) == 1:4
    @test searchsorted(R[1:10;], T(10), by=(x -> x >= 5)) == 5:10
    @test searchsorted(R[1:5; 1:5; 1:5], T(1), 6, 10, Forward) == 6:6
    @test searchsorted(ones(R, 15), T(1), 6, 10, Forward) == 6:10
end

for (rg,I) in [(49:57,47:59), (1:2:17,-1:19), (-3:0.5:2,-5:.5:4)]
    rg_r = reverse(rg)
    rgv, rgv_r = [rg;], [rg_r;]
    for i = I
        @test searchsorted(rg,i) == searchsorted(rgv,i)
        @test searchsorted(rg_r,i,rev=true) == searchsorted(rgv_r,i,rev=true)
    end
end

rg = 0.0:0.01:1.0
for i = 2:101
    @test searchsorted(rg, rg[i]) == i:i
    @test searchsorted(rg, prevfloat(rg[i])) == i:i-1
    @test searchsorted(rg, nextfloat(rg[i])) == i+1:i
end

rg_r = reverse(rg)
for i = 1:100
    @test searchsorted(rg_r, rg_r[i], rev=true) == i:i
    @test searchsorted(rg_r, prevfloat(rg_r[i]), rev=true) == i+1:i
    @test searchsorted(rg_r, nextfloat(rg_r[i]), rev=true) == i:i-1
end

@test searchsorted(1:10, 1, by=(x -> x >= 5)) == searchsorted([1:10;], 1, by=(x -> x >= 5))
@test searchsorted(1:10, 10, by=(x -> x >= 5)) == searchsorted([1:10;], 10, by=(x -> x >= 5))

@test searchsorted([], 0) == 1:0
@test searchsorted([1,2,3], 0) == 1:0
@test searchsorted([1,2,3], 4) == 4:3

# issue 8866
@test searchsortedfirst(500:1.0:600, -1.0e20) == 1
@test searchsortedfirst(500:1.0:600, 1.0e20) == 102
@test searchsortedlast(500:1.0:600, -1.0e20) == 0
@test searchsortedlast(500:1.0:600, 1.0e20) == 101

# exercise the codepath in searchsorted* methods for ranges that check for zero step range
immutable ConstantRange{T} <: Range{T}
   val::T
   len::Int
end

Base.length(r::ConstantRange) = r.len
Base.getindex(r::ConstantRange, i::Int) = (1 <= i <= r.len || throw(BoundsError(r,i)); r.val)
Base.step(r::ConstantRange) = 0

r = ConstantRange(1, 5)

@test searchsortedfirst(r, 1.0, Forward) == 1
@test searchsortedfirst(r, 1, Forward) == 1
@test searchsortedfirst(r, UInt(1), Forward) == 1

@test searchsortedlast(r, 1.0, Forward) == 5
@test searchsortedlast(r, 1, Forward) == 5
@test searchsortedlast(r, UInt(1), Forward) == 5

a = rand(1:10000, 1000)

for alg in [InsertionSort, MergeSort]
    b = sort(a, alg=alg)
    @test issorted(b)

    ix = sortperm(a, alg=alg)
    b = a[ix]
    @test issorted(b)
    @test a[ix] == b

    sortperm!(view(ix, 1:100), view(a, 1:100), alg=alg)
    b = a[ix][1:100]
    @test issorted(b)

    sortperm!(ix, a, alg=alg)
    b = a[ix]
    @test issorted(b)
    @test a[ix] == b

    b = sort(a, alg=alg, rev=true)
    @test issorted(b, rev=true)
    ix = sortperm(a, alg=alg, rev=true)
    b = a[ix]
    @test issorted(b, rev=true)
    @test a[ix] == b

    sortperm!(view(ix, 1:100), view(a, 1:100), alg=alg, rev=true)
    b = a[ix][1:100]
    @test issorted(b, rev=true)

    sortperm!(ix, a, alg=alg, rev=true)
    b = a[ix]
    @test issorted(b, rev=true)
    @test a[ix] == b

    b = sort(a, alg=alg, by=x->1/x)
    @test issorted(b, by=x->1/x)
    ix = sortperm(a, alg=alg, by=x->1/x)
    b = a[ix]
    @test issorted(b, by=x->1/x)
    @test a[ix] == b

    sortperm!(view(ix, 1:100), view(a, 1:100), by=x->1/x)
    b = a[ix][1:100]
    @test issorted(b, by=x->1/x)

    sortperm!(ix, a, alg=alg, by=x->1/x)
    b = a[ix]
    @test issorted(b, by=x->1/x)
    @test a[ix] == b

    c = copy(a)
    permute!(c, ix)
    @test c == b

    ipermute!(c, ix)
    @test c == a

    c = sort(a, alg=alg, lt=(>))
    @test b == c

    c = sort(a, alg=alg, by=x->1/x)
    @test b == c
end

# unstable algorithms
for alg in [QuickSort, PartialQuickSort(length(a))]
    b = sort(a, alg=alg)
    @test issorted(b)
    b = sort(a, alg=alg, rev=true)
    @test issorted(b, rev=true)
    b = sort(a, alg=alg, by=x->1/x)
    @test issorted(b, by=x->1/x)
end

# test PartialQuickSort only does a partial sort
let alg = PartialQuickSort(div(length(a), 10))
    k = alg.k
    b = sort(a, alg=alg)
    c = sort(a, alg=alg, by=x->1/x)
    d = sort(a, alg=alg, rev=true)
    @test issorted(b[1:k])
    @test issorted(c[1:k], by=x->1/x)
    @test issorted(d[1:k], rev=true)
    @test !issorted(b)
    @test !issorted(c, by=x->1/x)
    @test !issorted(d, rev=true)
end

@test select([3,6,30,1,9], 2, rev=true) == 9
@test select([3,6,30,1,9], 2, by=x->1/x) == 9
@test selectperm([3,6,30,1,9], 2, rev=true) == 5
@test selectperm([3,6,30,1,9], 2, by=x->1/x) == 5

## more advanced sorting tests ##

randnans(n) = reinterpret(Float64,[rand(UInt64)|0x7ff8000000000000 for i=1:n])

function randn_with_nans(n,p)
    v = randn(n)
    x = find(rand(n).<p)
    v[x] = randnans(length(x))
    return v
end

srand(0xdeadbeef)

for n in [0:10; 100; 101; 1000; 1001]
    r = -5:5
    v = rand(r,n)
    h = [sum(v .== x) for x in r]

    for rev in [false,true]
        # insertion sort (stable) as reference
        pi = sortperm(v, alg=InsertionSort, rev=rev)
        @test pi == sortperm(float(v), alg=InsertionSort, rev=rev)
        @test isperm(pi)
        si = v[pi]
        @test [sum(si .== x) for x in r] == h
        @test issorted(si, rev=rev)
        @test all(issorted,[pi[si.==x] for x in r])
        c = copy(v)
        permute!(c, pi)
        @test c == si
        ipermute!(c, pi)
        @test c == v

        # stable algorithms
        for alg in [MergeSort]
            p = sortperm(v, alg=alg, rev=rev)
            @test p == sortperm(float(v), alg=alg, rev=rev)
            @test p == pi
            s = copy(v)
            permute!(s, p)
            @test s == si
            ipermute!(s, p)
            @test s == v
        end

        # unstable algorithms
        for alg in [QuickSort, PartialQuickSort(n)]
            p = sortperm(v, alg=alg, rev=rev)
            @test p == sortperm(float(v), alg=alg, rev=rev)
            @test isperm(p)
            @test v[p] == si
            s = copy(v)
            permute!(s, p)
            @test s == si
            ipermute!(s, p)
            @test s == v
        end
    end

    v = randn_with_nans(n,0.1)
    # TODO: alg = PartialQuickSort(n) fails here
    for alg in [InsertionSort, QuickSort, MergeSort],
        rev in [false,true]
        # test float sorting with NaNs
        s = sort(v, alg=alg, rev=rev)
        @test issorted(s, rev=rev)
        @test reinterpret(UInt64,v[isnan.(v)]) == reinterpret(UInt64,s[isnan.(s)])

        # test float permutation with NaNs
        p = sortperm(v, alg=alg, rev=rev)
        @test isperm(p)
        vp = v[p]
        @test isequal(vp,s)
        @test reinterpret(UInt64,vp) == reinterpret(UInt64,s)
    end
end

inds = [
    1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6,7,7,7,8,8,8,9,9,9,10,
    10,10,11,11,11,12,12,12,13,13,13,14,14,14,15,15,15,16,16,
    16,17,17,17,18,18,18,19,19,19,20,20,20,21,21,22,22,22,23,
    23,24,24,24,25,25,25,26,26,26,27,27,27,28,28,28,29,29,29,
    30,30,30,31,31,32,32,32,33,33,33,34,34,34,35,35,35,36,36,
    36,37,37,37,38,38,38,39,39,39,40,40,40,41,41,41,42,42,42,
    43,43,43,44,44,44,45,45,45,46,46,46,47,47,47,48,48,48,49,
    49,49,50,50,50,51,51,52,52,52,53,53,53,54,54,54,55,55,55,
    56,56,56,57,57,57,58,58,58,59,60,60,60,61,61,61,62,62,63,
    64,64,64,65,65,65,66,66,66,67,67,67,68,68,69,69,69,70,70,
    70,71,71,71,72,72,72,73,73,73,74,75,75,75,76,76,76,77,77,
    77,78,78,78,79,79,79,80,80,80,81,81,82,82,82,83,83,83,84,
    84,84,85,85,85,86,86,86,87,87,87,88,88,88,89,89,89,90,90,
    90,91,91,91,92,92,93,93,93,94,94,94,95,95,95,96,96,96,97,
    97,98,98,98,99,99,99,100,100,100,101,101,101,102,102,102,
    103,103,103,104,105,105,105,106,106,106,107,107,107,108,
    108,108,109,109,109,110,110,110,111,111,111,112,112,112,
    113,113,113,114,114,115,115,115,116,116,116,117,117,117,
    118,118,118,119,119,119,120,120,120,121,121,121,122,122,
    122,123,123,123,124,124,124,125,125,125,126,126,126,127,
    127,127,128,128,128,129,129,129,130,130,130,131,131,131,
    132,132,132,133,133,133,134,134,134,135,135,135,136,136,
    136,137,137,137,138,138,138,139,139,139,140,140,140,141,
    141,142,142,142,143,143,143,144,144,144,145,145,145,146,
    146,146,147,147,147,148,148,148,149,149,149,150,150,150,
    151,151,151,152,152,152,153,153,153,154,154,154,155,155,
    155,156,156,156,157,157,157,158,158,158,159,159,159,160,
    160,160,161,161,161,162,162,162,163,163,163,164,164,164,
    165,165,165,166,166,166,167,167,167,168,168,168,169,169,
    169,170,170,170,171,171,171,172,172,172,173,173,173,174,
    174,174,175,175,175,176,176,176,177,177,177,178,178,178,
    179,179,179,180,180,180,181,181,181,182,182,182,183,183,
    183,184,184,184,185,185,185,186,186,186,187,187,187,188,
    188,188,189,189,189,190,190,190,191,191,191,192,192,192,
    193,193,193,194,194,194,195,195,195,196,196,197,197,197,
    198,198,198,199,199,199,200,200,200
]
sp = sortperm(inds)
@test all(issorted, [sp[inds.==x] for x in 1:200])
for alg in [InsertionSort, MergeSort]
    sp = sortperm(inds, alg=alg)
    @test all(issorted, [sp[inds.==x] for x in 1:200])
end

# issue #6177
@test sortperm([ 0.0, 1.0, 1.0], rev=true) == [2, 3, 1]
@test sortperm([-0.0, 1.0, 1.0], rev=true) == [2, 3, 1]
@test sortperm([-1.0, 1.0, 1.0], rev=true) == [2, 3, 1]

# issue #8825 - stability of min/max
type Twain
    a :: Int
    b :: Int
end
Base.isless(x :: Twain, y :: Twain) = x.a < y.a
let x = Twain(2,3), y = Twain(2,4)
    @test (min(x,y), max(x,y)) == (x,y) == minmax(x,y)
end

# issue #12833 - type stability of sort
@test Base.return_types(sort, (Vector{Int},)) == [Vector{Int}]

# PR #18791
@test sort([typemax(Int),typemin(Int)]) == [typemin(Int),typemax(Int)]
@test sort([typemax(UInt),0]) == [0,typemax(UInt)]

# issue #19005
@test searchsortedfirst(0:256, 0x80) == 129
@test searchsortedlast(0:256, 0x80) == 129

max_stars_repo_name: MichaelHatherly/julia

max_issues_repo_name: MichaelHatherly/julia

max_forks_repo_name: MichaelHatherly/julia


text:
using QuantumControl
using Test

@testset "QuantumControl.jl" begin
    # Write your tests here.
end

ext: jl


text:
using ITensors
using ITensorInfiniteMPS
using Test
using Random

@testset "vumps_extended_ising" begin
  Random.seed!(1234)

  model = Model"ising_extended"()
  model_kwargs = (J=1.0, h=1.1, J₂=0.2)
  initstate(n) = "↑"

  function space_shifted(::Model"ising_extended", q̃sz; conserve_qns=true)
    if conserve_qns
      return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1]
    else
      return [QN() => 2]
    end
  end

  # Compare to DMRG
  Nfinite = 100
  sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
  Hfinite = MPO(model, sfinite; model_kwargs...)
  ψfinite = randomMPS(sfinite, initstate)
  sweeps = Sweeps(20)
  setmaxdim!(sweeps, 30)
  setcutoff!(sweeps, 1E-10)
  energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps; outputlevel=0)
  Szs_finite = expect(ψfinite, "Sz")

  function energy(ψ, h, n)
    ϕ = ψ[n] * ψ[n + 1] * ψ[n + 2]
    return (noprime(ϕ * h) * dag(ϕ))[]
  end

  nfinite = Nfinite ÷ 2
  hnfinite = ITensor(
    model, sfinite[nfinite], sfinite[nfinite + 1], sfinite[nfinite + 2]; model_kwargs...
  )
  orthogonalize!(ψfinite, nfinite)
  energy_finite = energy(ψfinite, hnfinite, nfinite)

  cutoff = 1e-8
  maxdim = 100
  tol = 1e-8
  maxiter = 20
  outer_iters = 4
  for multisite_update_alg in ["sequential", "parallel"],
    conserve_qns in [true, false],
    nsites in [1, 2],
    time_step in [-Inf, -0.5]

    vumps_kwargs = (
      multisite_update_alg=multisite_update_alg,
      tol=tol,
      maxiter=maxiter,
      outputlevel=0,
      time_step=time_step,
    )
    subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)

    space_ = fill(space_shifted(model, 1; conserve_qns=conserve_qns), nsites)
    s = infsiteinds("S=1/2", nsites; space=space_)
    H = InfiniteSum{MPO}(model, s; model_kwargs...)
    ψ = InfMPS(s, initstate)

    # Check translational invariance
    @test contract(ψ.AL[1:nsites]..., ψ.C[nsites]) ≈ contract(ψ.C[0], ψ.AR[1:nsites]...)

    ψ = tdvp(H, ψ; vumps_kwargs...)
    for _ in 1:outer_iters
      ψ = subspace_expansion(ψ, H; subspace_expansion_kwargs...)
      ψ = tdvp(H, ψ; vumps_kwargs...)
    end
    # Check translational invariance
    @test norm(
      contract(ψ.AL[1:nsites]..., ψ.C[nsites]) - contract(ψ.C[0], ψ.AR[1:nsites]...)
    ) ≈ 0 atol = 1e-5

    energy_infinite = expect(ψ, H)
    Szs_infinite = [expect(ψ, "Sz", n) for n in 1:nsites]

    @test energy_finite ≈ sum(energy_infinite) / nsites rtol = 1e-4
    @test Szs_finite[nfinite:(nfinite + nsites - 1)] ≈ Szs_infinite rtol = 1e-3
  end
end

# XXX: orthogonalize is broken right now
## @testset "ITensorInfiniteMPS.jl" begin
##   @testset "Mixed canonical gauge" begin
##     N = 10
##     s = siteinds("S=1/2", N; conserve_szparity=true)
##     χ = 6
##     @test iseven(χ)
##     space = (("SzParity", 1, 2) => χ ÷ 2) ⊕ (("SzParity", 0, 2) => χ ÷ 2)
##     ψ = InfiniteMPS(ComplexF64, s; space=space)
##     randn!.(ψ)
##
##     ψ = orthogonalize(ψ, :)
##     @test prod(ψ.AL[1:N]) * ψ.C[N] ≈ ψ.C[0] * prod(ψ.AR[1:N])
##   end
## end


text:
using Flux
using KitBase.JLD2
using Flux: @epochs

cd(@__DIR__)
include("../common.jl")
@load "data/dataset.jld2" X Y
#@load "model/nn_pro.jld2" nn
@load "model/nn_layer.jld2" nn

data = Flux.Data.DataLoader((X, Y), shuffle = true)
ps = params(nn)
loss(x, y) = Flux.binarycrossentropy(nn(x), y)
cb = () -> println("loss: $(loss(X, Y))")
opt = ADAM()

@epochs 5 Flux.train!(loss, ps, data, opt, cb = Flux.throttle(cb, 1))

accuracy(nn, X, Y)

@save "model/nn_layer.jld2" nn

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
twoband(ρ) = UnionDomain(-1..(-ρ), ρ..1)

"""
    TwoBandWeight(a, b)

is a quasi-vector representing `|x|^(2c) * (x^2-ρ)^b * (1-x^2)^a`
"""
struct TwoBandWeight{T} <: Weight{T}
    ρ::T
    a::T
    b::T
    c::T
end

TwoBandWeight(ρ::R, a::A, b::B, c::C) where {R,A,B,C} = TwoBandWeight{promote_type(R,A,B,C)}(ρ, a, b, c)

copy(w::TwoBandWeight) = w

axes(w::TwoBandWeight{T}) where T = (Inclusion(twoband(w.ρ)),)

==(w::TwoBandWeight, v::TwoBandWeight) = w.a == v.a && w.b == v.b && w.ρ == v.ρ && w.c == v.c

function getindex(w::TwoBandWeight, x::Real)
    @boundscheck checkbounds(w, x)
    abs(x)^(2w.c) * (x^2- w.ρ^2)^w.b * (1-x^2)^w.a
end

function sum(w::TwoBandWeight{T}) where T
    if 2w.a == 2w.b == -1 && 2w.c == 1
        convert(T,π)/2
    else
        error("not implemented")
#         1/2 \[Pi] (-((\[Rho]^(1 + 2 b + 2 c)
#       Gamma[1 + b] Hypergeometric2F1Regularized[-a, 1/2 + c,
#       3/2 + b + c, \[Rho]^2])/Gamma[1/2 - c]) + (
#    Gamma[1 + a] Hypergeometric2F1Regularized[-b, -(1/2) - a - b - c,
#      1/2 - b - c, \[Rho]^2])/
#    Gamma[3/2 + a + b + c]) Sec[(b + c) \[Pi]]
    end
end

"""
    TwoBandJacobi(a, b)

is a quasi-matrix orthogonal `|x|^(2c) * (x^2 - ρ^2)^b * (1-x^2)^a`.
"""
struct TwoBandJacobi{T} <: OrthogonalPolynomial{T}
    ρ::T
    a::T
    b::T
    c::T
    P::SemiclassicalJacobi{T}
    Q::SemiclassicalJacobi{T}
    TwoBandJacobi{T}(ρ::T, a::T, b::T, c::T, P::SemiclassicalJacobi{T}, Q::SemiclassicalJacobi{T}) where T = new{T}(ρ, a, b, c, P, Q)
end
function TwoBandJacobi{T}(ρ, a, b, c) where T
    t = inv(1-ρ^2)
    P = SemiclassicalJacobi(t, a, b, c-one(T)/2)
    TwoBandJacobi{T}(convert(T,ρ), convert(T,a), convert(T,b), convert(T,c), P, SemiclassicalJacobi(t, a, b, c+one(T)/2, P))
end
TwoBandJacobi(ρ, a, b, c) = TwoBandJacobi{float(promote_type(eltype(ρ),eltype(a),eltype(b),eltype(c)))}(ρ, a, b, c)

axes(P::TwoBandJacobi{T}) where T = (Inclusion(twoband(P.ρ)),oneto(∞))

==(w::TwoBandJacobi, v::TwoBandJacobi) = w.ρ == v.ρ && w.a == v.a && w.b == v.b && w.c == v.c

copy(A::TwoBandJacobi) = A

orthogonalityweight(Z::TwoBandJacobi) = TwoBandWeight(Z.ρ, Z.a, Z.b, Z.c)

function  getindex(R::TwoBandJacobi, x::Real, j::Integer)
    ρ = R.ρ
    if isodd(j)
        R.P[(1-x^2)/(1-ρ^2), (j+1)÷2]
    else
        x * R.Q[(1-x^2)/(1-ρ^2), j÷2]
    end
end

struct Interlace{T,AA,BB} <: LazyVector{T}
    a::AA
    b::BB
end

Interlace{T}(a::AbstractVector{T}, b::AbstractVector{T}) where T = Interlace{T,typeof(a),typeof(b)}(a,b)
Interlace(a::AbstractVector{T}, b::AbstractVector{V}) where {T,V} = Interlace{promote_type(T,V)}(a, b)

size(::Interlace) = (ℵ₀,)

getindex(A::Interlace{T}, k::Int) where T = convert(T, isodd(k) ? A.a[(k+1)÷2] : A.b[k÷2])::T

function jacobimatrix(R::TwoBandJacobi{T}) where T
    ρ = R.ρ
    L = R.P \ (SemiclassicalJacobiWeight(R.P.t,0,0,1) .* R.Q)
    Tridiagonal(Interlace(L.dv/L[1,1], (1-ρ^2) * L.ev), Zeros{T}(∞), Interlace((1-ρ^2) * L.dv,L.ev/L[1,1]))
end

@simplify function *(H::Hilbert, w::TwoBandWeight)
    if 2w.a == 2w.b == -1 && 2w.c == 1
        zeros(promote_type(eltype(H),eltype(w)), axes(w,1))
    else
       error("Not Implemented")
    end
end

function grid(L::SubQuasiArray{T,2,<:Associated{<:Any,<:TwoBandJacobi},<:Tuple{Inclusion,OneTo}}) where T
    g = grid(legendre(parent(L).P.ρ .. 1)[:,parentindices(L)[2]])
    [-g; g]
end

LinearAlgebra.factorize(L::SubQuasiArray{<:Any,2,<:Associated{<:Any,<:TwoBandJacobi},<:Tuple{Inclusion,OneTo}}) =
    ContinuumArrays._factorize(BasisLayout(), L)


text:
@testset "Cross validating LMHMM on ZX..." begin
        @time ctb = ZhuXian()

        nb = 10
        batches = KongYiji.foldbatch(ctb, nb)
        tbs = KongYiji.HmmScoreTable[]
        @showprogress 1 "Cross Validating ooHMM..." for (tr, te) in batches
                hmm = KongYiji.HMM(tr, "", 0.)
                lm = KongYiji.LM(tr; epochs=20, rnnType=:lstm, numLayers=1, bidirectional=false)
                m = hmm + lm
                x = KongYiji.rawsents(te)
                z = KongYiji.wordsents(te)
                @time y = m(x; beams=500, recover=false, withpos=false)
                push!(tbs, KongYiji.HmmScoreTable(z, y))
                println(tbs[end])
        end
        hmmscoretable = sum(tbs)
        println(hmmscoretable)
end


text:
#= ForwardDiff is removed as a dependency from this version
using Revise
import MechanicalSketch: @import_expand, empty_figure, WI, HE, EM, O, PT, FS, finish
import MechanicalSketch: settext, place_image, PALETTE, color_with_lumin, circle, line
import MechanicalSketch: circle, ∙, arrow_nofill, Point, @layer, sethue, draw_expr
import MechanicalSketch: latexify, arrow, Quantity, readsvg, unit, ustrip
import MechanicalSketch: ForwardDiff, box_fill_outline, brush, setopacity
import MechanicalSketch: @ev_draw, rotate_hue, setline
import MechanicalSketch: Drawing, configure_mechanical, preview
import MechanicalSketch: Luxor
import Luxor:            CairoContext, CairoRecordingSurface, flush, set_source
import Luxor:            paint, getmatrix, setmatrix, currentpoint, get_current_cr
import Luxor:            current_surface_type, CURRENTDRAWING
import Luxor.Cairo:      device_get_type, flush

if !@isdefined N
    @import_expand(~m, ~s, °, N)
end

"""
    make_recording_surface_drawing()
    --> Luxor.Drawing

This makes, then mutates a Luxor drawing to allow
a special type of Cairo surface.
"""
function make_recording_surface_drawing()
    d = Drawing(WI, HE, :svg)
    s1 = CairoRecordingSurface()
    c1 = CairoContext(s1)
    d.surface = s1
    d.cr = c1
    configure_mechanical()
    d
end

"""
    snapshot_of_working_surface(fnam)

Takes a snapshot and saves to file with 'fnam' name and suffix.
One could continue drawing, or do the other things.
"""
function snapshot_of_working_surface(fnam)
    drec = currentdrawing()
    @assert typeof(drec.surface)
    flush(drec.surface)
    d = empty_figure(fnam)
    d.width = drec.width
    d.height = drec.height
    d.redvalue = drec.redvalue
    d.greenvalue = drec.greenvalue
    d.bluevalue = drec.bluevalue
    d.alpha = drec.alpha
    set_source(d.cr, drec.surface, -WI / 2, -HE / 2)
    paint()
    finish()
    CURRENTDRAWING[1] = drec
end

make_recording_surface_drawing()

for r = 1m:0.2m:60m
   α = 360.0° * (r / 20m)
   circle(O + 100 .* (cos(α), sin(α)); r)
   setline(8)
   sethue(rotate_hue(PALETTE[5], α))
   settext(string(r), Point(0m, r))
end
snapshotrecordingsurface("snap1.png")
for r = 1m:0.2m:60m
    α = 360.0° * (r / 20m)
    circle(O + 100 .* (cos(-α), sin(α)); r)
    setline(8)
    sethue(rotate_hue(PALETTE[3], α))
    settext(string(r), Point(0m, r))
 end
 snap2 = snapshotrecordingsurface(drec)
#snap1 = snapshotrecordingsurface(drec)
#circle(O + (3m, 0m), r = 5m)
#snap2 = snapshotrecordingsurface(drec)

=#


text:
function PolynomialNutation(arg0::Vector{jdouble})
    return PolynomialNutation((Vector{jdouble},), arg0)
end

function derivative(obj::PolynomialNutation, arg0::RealFieldElement)
    return jcall(obj, "derivative", RealFieldElement, (RealFieldElement,), arg0)
end

function derivative(obj::PolynomialNutation, arg0::jdouble)
    return jcall(obj, "derivative", jdouble, (jdouble,), arg0)
end

function value(obj::PolynomialNutation, arg0::RealFieldElement)
    return jcall(obj, "value", RealFieldElement, (RealFieldElement,), arg0)
end

function value(obj::PolynomialNutation, arg0::jdouble)
    return jcall(obj, "value", jdouble, (jdouble,), arg0)
end


text:
#Define structure for testing
Model = HA_solver.Ayiagari_params()

Model.a_min = 0.0 #minimal number of points
Model.MC = "T" #use Tauchen method

Model.start_r =0.02

@testset "UnpackAyiagari_params" begin
  #test for Unpacking Ayiagari, check if params are well written to the tuple and if they make sense
  Model.β = 2.0
  @test_throws AssertionError HA_solver.UnpackAyiagari_params(Model)
  Model.β = 0.96

  Model.σ = 0.0
  @test_throws AssertionError HA_solver.UnpackAyiagari_params(Model)
  Model.σ = 2.0

  Model.rho = -1.0
  @test_throws AssertionError HA_solver.UnpackAyiagari_params(Model)
  Model.rho = 0.9

  Model.sigma_e = 0.0
  @test_throws AssertionError HA_solver.UnpackAyiagari_params(Model)
  Model.sigma_e = 0.2

  Model.Δ = -0.01
  @test_throws AssertionError HA_solver.UnpackAyiagari_params(Model)
  Model.Δ = 0.01

  Model.g = -0.01
  @test_throws AssertionError HA_solver.UnpackAyiagari_params(Model)
  Model.g = 0.015

  Model.a_n = 2
  @test_throws AssertionError HA_solver.UnpackAyiagari_params(Model)
  Model.a_n = 500

  Model.ϵ_n = -2
  @test_throws AssertionError HA_solver.UnpackAyiagari_params(Model)
  Model.ϵ_n = 9

  Model.α = -0.5
  @test_throws AssertionError HA_solver.UnpackAyiagari_params(Model)
  Model.α = 0.36

  Model.δ = -0.08
  @test_throws AssertionError HA_solver.UnpackAyiagari_params(Model)
  Model.δ = 0.08

  Model.γ = -0.95
  @test_throws AssertionError HA_solver.UnpackAyiagari_params(Model)
  Model.γ = 0.95
#write down params for the future test
end

@testset "UpdateGEq_Ayiagari_r" begin
    NumParams, Params = HA_solver.UnpackAyiagari_params(Model)
    SolParam_old = []
    #check if the negative capital throw the error
    @test_throws AssertionError HA_solver.UpdateGEq_Ayiagari_r(Params,SolParam_old, -1.0, nothing, 0)

end

@testset "ConvSolParam_Ayiagari_r" begin
    SolParam_old= nothing

    @test_throws AssertionError HA_solver.ConvSolParam_Ayiagari_r(SolParam_old, 0.01, 1.0, 2000, 2)

end

@testset "SaveHAEq_Ayiagari" begin
    # I don't test this savings, as it's just some savings
    @test 1==1
end

@testset "SolveAgP_Ayiagari_EGM and SolveDistr_Ayiagari_Iter" begin
      NumParams, Params = HA_solver.UnpackAyiagari_params(Model)
      test_r = 0.0:0.001:0.03

    #check if the problem gest solved
    for i in eachindex(test_r)
        Policy = HA_solver.SolveAgP_Ayiagari_EGM(Params, NumParams, test_r[i])
        Measure, k = HA_solver.SolveDistr_Ayiagari_Iter(Params, NumParams,test_r[i], Policy )
        @test Policy[4] < 1e-7
        @test abs(sum(Measure) -1.0) <=1e-9

    end

end

#Next test if results of Ayiagari model are the same as in the paper. Is it not possible to  exactly replicate the papaer, due to the different grid points, simulation technique etc.
# we check if the fixed point interest rate is not different than 0.5 percent points
@testset "Find_eq_Ayiagari" begin
    σ = 1.0
    ρ_set = [0.0, 0.3, 0.6, 0.9]
    #orginal solution from Ayiagari
    r_computed = zeros(4)
    orginal_solution = [0.04166, 0.041365, 0.040912, 0.03905]
        for j in eachindex(ρ_set)
            #get decent initial guies for this test to not take forever
            Model.γ = 0.99
            Model.start_r =  0.04
            Model.σ = 1.0
            Model.rho = ρ_set[j]
            r_computed[j] = HA_solver.Find_eq_Ayiagari(Model)
            @test abs(r_computed[j] -orginal_solution[j]) <=0.005
        end
end


text:
module Bridge

using ..Ahorn, Maple

function bridgeFinalizer(entity::Maple.Bridge)
    x, y = Ahorn.position(entity)

    entity.data["nodes"] = [
        (x + 96, y),
        (x + 128, y),
    ]
end

const placements = Ahorn.PlacementDict(
    "Bridge" => Ahorn.EntityPlacement(
        Maple.Bridge,
        "rectangle",
        Dict{String, Any}(),
        bridgeFinalizer
    ),
    "Bridge (Fixed)" => Ahorn.EntityPlacement(
        Maple.BridgeFixed,
        "rectangle"
    )
)

const bridgeSprite = "scenery/bridge"
const bridgeFixedSprite = "scenery/bridge_fixed"

Ahorn.minimumSize(entity::Maple.Bridge) = 32, 0
Ahorn.minimumSize(entity::Maple.BridgeFixed) = Ahorn.getSprite(bridgeFixedSprite, "Gameplay").width, 0

Ahorn.resizable(entity::Maple.Bridge) = true, false
Ahorn.resizable(entity::Maple.BridgeFixed) = true, false

Ahorn.nodeLimits(entity::Maple.Bridge) = 2, 2

function Ahorn.selection(entity::Maple.Bridge)
    x, y = Ahorn.position(entity)

    sprite = Ahorn.getSprite(bridgeSprite, "Gameplay")
    nodes = get(entity.data, "nodes", ())

    if length(nodes) == 2
        gapStartX = nodes[1][1]
        gapEndX = nodes[2][1]

        return [
            Ahorn.Rectangle(x, y - 3, bridgeSelectionWidth(entity), sprite.height),
            Ahorn.Rectangle(gapStartX - 4, y - 16, 8, 32),
            Ahorn.Rectangle(gapEndX - 4, y - 16, 8, 32)
        ]
    end
end

function Ahorn.selection(entity::Maple.BridgeFixed)
    x, y = Ahorn.position(entity)

    sprite = Ahorn.getSprite(bridgeFixedSprite, "Gameplay")
    width = get(entity.data, "width", 32)

    selectionWidth = ceil(Int, (width / sprite.width)) * sprite.width

    return Ahorn.Rectangle(x, y - 4, selectionWidth, sprite.height)
end

const sizes = Ahorn.Rectangle[
    Ahorn.Rectangle(0, 0, 16, 55),
    Ahorn.Rectangle(16, 0, 8, 55),
    Ahorn.Rectangle(24, 0, 8, 55),
    Ahorn.Rectangle(32, 0, 8, 55),
    Ahorn.Rectangle(40, 0, 8, 55),
    Ahorn.Rectangle(48, 0, 8, 55),
    Ahorn.Rectangle(56, 0, 8, 55),
    Ahorn.Rectangle(64, 0, 8, 55),
    Ahorn.Rectangle(72, 0, 8, 55),
    Ahorn.Rectangle(80, 0, 16, 55),
    Ahorn.Rectangle(96, 0, 8, 55)
]

function renderBridgeTile(ctx::Ahorn.Cairo.CairoContext, x::Integer, y::Integer, size::Ahorn.Rectangle)
    tileWidth = size.w
    tileHeight = size.h

    if tileWidth == 16
        height = 24
        py = 0

        while py < tileHeight
            Ahorn.drawImage(ctx, bridgeSprite, x, y + py - 3, size.x, py, tileWidth, min(height, tileHeight - py))

            py += height
            height = 12
        end

    else
        Ahorn.drawImage(ctx, bridgeSprite, x, y - 3, size.x, size.y, tileWidth, tileHeight)
    end
end

# Helper for selections
function bridgeSelectionWidth(entity::Maple.Bridge)
    x, y = Ahorn.position(entity)

    width = get(entity.data, "width", 32)
    nodes = get(entity.data, "nodes", ())

    if length(nodes) == 2
        rng = Ahorn.getSimpleEntityRng(entity)

        gapStartX = nodes[1][1]
        gapEndX = nodes[2][1]

        index = 1
        px, py = x, y

        while px < x + width
            tileSize = index < 3 || index > 8 ? sizes[index] : sizes[3 + rand(rng, 0:6)]

            px += tileSize.w
            index = mod1(index + 1, length(sizes))
        end
    end

    return px - x
end

function renderBridge(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Bridge)
    x, y = Ahorn.position(entity)

    width = get(entity.data, "width", 32)
    nodes = get(entity.data, "nodes", ())

    if length(nodes) == 2
        rng = Ahorn.getSimpleEntityRng(entity)

        gapStartX = nodes[1][1]
        gapEndX = nodes[2][1]

        index = 1
        px, py = x, y

        while px < x + width
            tileSize = index < 3 || index > 8 ? sizes[index] : sizes[3 + rand(rng, 0:6)]

            if px < gapStartX || px >= gapEndX
                renderBridgeTile(ctx, px, py, tileSize)
            end

            px += tileSize.w
            index = mod1(index + 1, length(sizes))
        end
    end
end

function renderFixedBridge(ctx::Ahorn.Cairo.CairoContext, entity::Maple.BridgeFixed)
    x, y = Ahorn.position(entity)

    width = get(entity.data, "width", 32)

    px = x
    sprite = Ahorn.getSprite(bridgeFixedSprite, "Gameplay")

    while px < x + width
        Ahorn.drawSprite(ctx, bridgeFixedSprite, px, y - 8, jx=0.0, jy=0.0)

        px += sprite.width
    end
end

Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Bridge, room::Maple.Room) = renderBridge(ctx, entity)
Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.BridgeFixed, room::Maple.Room) = renderFixedBridge(ctx, entity)

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Bridge, room::Maple.Room)
    x, y = Ahorn.position(entity)
    nodes = get(entity.data, "nodes", ())

    width = Int(get(entity.data, "width", 8))

    if length(nodes) == 2
        gapStartX = nodes[1][1]
        gapEndX = nodes[2][1]

        Ahorn.drawRectangle(ctx, gapStartX, y - 16, 1, 32, (1.0, 0.0, 0.0, 1.0), (0.0, 0.0, 0.0, 0.0))
        Ahorn.drawRectangle(ctx, gapEndX, y - 16, 1, 32, (1.0, 0.0, 0.0, 1.0), (0.0, 0.0, 0.0, 0.0))
    end
end

end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text: # TODO (max): how to best handle io: HDF5? or something else

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
function build_op_problem!(op_problem::OperationsProblem{M}; kwargs...) where M<:AbstractOperationsProblem
    sys = get_system(op_problem)
    _build!(op_problem.canonical, op_problem.template, sys; kwargs...)
    return
end

function _build!(canonical::Canonical, template::OperationsTemplate, sys::PSY.System; kwargs...)
    verbose = get(kwargs, :verbose, true)
    transmission = template.transmission
    # Build Injection devices
    for device_model in values(template.devices)
        verbose && @info "Building $(device_model.device_type) with $(device_model.formulation) formulation"
        construct_device!(canonical, sys, device_model, transmission; kwargs...)
    end

    # Build Network
    verbose && @info "Building $(transmission) network formulation"
    construct_network!(canonical, sys, transmission; kwargs...)

    # Build Branches
    for branch_model in values(template.branches)
        verbose && @info "Building $(branch_model.device_type) with $(branch_model.formulation) formulation"
        construct_device!(canonical, sys, branch_model, transmission; kwargs...)
    end

    #Build Services
    construct_services!(canonical, sys, template.services; kwargs...)

    # Objective Function
    verbose && @info "Building Objective"
    JuMP.@objective(canonical.JuMPmodel, MOI.MIN_SENSE, canonical.cost_function)

    return

end


text:
@testset "symmetric" begin
    @testset "Symmetric" begin
        rng, N = MersenneTwister(123456), 3
        rrule_test(Symmetric, randn(rng, N, N), (randn(rng, N, N), randn(rng, N, N)))
    end
end


text:
"""
    OpConnectionIdentity

An OpConnection to store the operatorial identity.
It always returns one entry with value `true` and no changes.
"""
struct OpConnectionIdentity{B,C} <: AbsOpConnection
    c::StateChanges{B,C}
end

OpConnectionIdentity(a::Type, b::Type) =
    OpConnectionIdentity(StateChanges{a,b}())

OpConnectionIdentity{B,C}() where {B,C} = OpConnectionIdentity{B,C}(StateChanges{B,C}())

function OpConnectionIdentity(op::AbsLinearOperator)
    t = conn_type(op)
    v1, v2, v3 = t.parameters
    return OpConnectionIdentity(v2, v3)
end

@inline Base.length(c::OpConnectionIdentity) = true
@inline Base.resize!(c::OpConnectionIdentity, i) = c
@inline Base.append!(c::OpConnectionIdentity, a::Nothing) = c
@inline Base.eltype(::Type{OpConnectionIdentity{B,C}}) where {B,C} =
    Tuple{Bool, StateChanges{B,C}}

Base.:(==)(::OpConnectionIdentity, ::OpConnectionIdentity) = true

@inline function Base.iterate(iter::OpConnectionIdentity, state=(1))
    if state == 1
        return (true, iter.c)
    else
        return nothing
    end
end

@inline function Base.getindex(c::OpConnectionIdentity, i)
    @assert i == 1
    return (true, c.c)
end

#showing
function Base.show(io::IO, c::OpConnectionIdentity{B,C}) where {B,C}
    print(io, "$(length(c))-elements - Identity-OpConnection{bool,StateChanges{$B,$C}}:\n")
    print(io, " (true, [])")
    return io
end


text:
function Pressure(RhoTh,Rho,KE,Param)
    (;Equation,
        Thermo,
        Rd,
        Cvd,
        Grav,
        Grid,
        p0,
        kappa) = Param
    str = Equation
    if str == "Compressible"
        if strcmp(Thermo,"Energy")
          p=(Rd/Cvd)*(RhoTh-Rho.*(KE+Grav*repmat(Grid.zP,1,size(Rho,1))'));
        else
          p=p0*(Rd*RhoTh/p0).^(1.0e0/(1.0e0-kappa));
        end
    elseif str == "Shallow"
        p=0.5*Grav*RhoTh.^2;
    end
    return p
end


text:
# Autogenerated wrapper script for CSPICE_jll for aarch64-linux-gnu
export brief, chronos, ckbrief, commnt, dskbrief, dskexp, frmdiff, inspekt, libcspice, libcsupport, mkdsk, mkspk, msopck, spacit, spkdiff, spkmerge, tobin, toxfr, version

JLLWrappers.@generate_wrapper_header("CSPICE")
JLLWrappers.@declare_executable_product(brief)
JLLWrappers.@declare_executable_product(chronos)
JLLWrappers.@declare_executable_product(ckbrief)
JLLWrappers.@declare_executable_product(commnt)
JLLWrappers.@declare_executable_product(dskbrief)
JLLWrappers.@declare_executable_product(dskexp)
JLLWrappers.@declare_executable_product(frmdiff)
JLLWrappers.@declare_executable_product(inspekt)
JLLWrappers.@declare_library_product(libcspice, "libcspice.so")
JLLWrappers.@declare_library_product(libcsupport, "libcsupport.so")
JLLWrappers.@declare_executable_product(mkdsk)
JLLWrappers.@declare_executable_product(mkspk)
JLLWrappers.@declare_executable_product(msopck)
JLLWrappers.@declare_executable_product(spacit)
JLLWrappers.@declare_executable_product(spkdiff)
JLLWrappers.@declare_executable_product(spkmerge)
JLLWrappers.@declare_executable_product(tobin)
JLLWrappers.@declare_executable_product(toxfr)
JLLWrappers.@declare_executable_product(version)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        brief,
        "bin/brief",
    )

    JLLWrappers.@init_executable_product(
        chronos,
        "bin/chronos",
    )

    JLLWrappers.@init_executable_product(
        ckbrief,
        "bin/ckbrief",
    )

    JLLWrappers.@init_executable_product(
        commnt,
        "bin/commnt",
    )

    JLLWrappers.@init_executable_product(
        dskbrief,
        "bin/dskbrief",
    )

    JLLWrappers.@init_executable_product(
        dskexp,
        "bin/dskexp",
    )

    JLLWrappers.@init_executable_product(
        frmdiff,
        "bin/frmdiff",
    )

    JLLWrappers.@init_executable_product(
        inspekt,
        "bin/inspekt",
    )

    JLLWrappers.@init_library_product(
        libcspice,
        "lib/libcspice.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcsupport,
        "lib/libcsupport.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        mkdsk,
        "bin/mkdsk",
    )

    JLLWrappers.@init_executable_product(
        mkspk,
        "bin/mkspk",
    )

    JLLWrappers.@init_executable_product(
        msopck,
        "bin/msopck",
    )

    JLLWrappers.@init_executable_product(
        spacit,
        "bin/spacit",
    )

    JLLWrappers.@init_executable_product(
        spkdiff,
        "bin/spkdiff",
    )

    JLLWrappers.@init_executable_product(
        spkmerge,
        "bin/spkmerge",
    )

    JLLWrappers.@init_executable_product(
        tobin,
        "bin/tobin",
    )

    JLLWrappers.@init_executable_product(
        toxfr,
        "bin/toxfr",
    )

    JLLWrappers.@init_executable_product(
        version,
        "bin/version",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
# module Commands

# Examples

```jldoctest
julia>
```
"""
module Commands

include("sbatch.jl")

end


text:
import TileMaps as TM
import Test

Test.@testset "TileMaps.jl" begin
    objects = (TM.EXAMPLE_OBJECT_1, TM.EXAMPLE_OBJECT_2, TM.EXAMPLE_OBJECT_3)

    layer_1 = BitArray([0 1 0 0 0
                        0 1 0 0 0
                        0 1 0 0 0
                        0 1 0 0 0])

    layer_2 = BitArray([0 0 0 0 0
                        0 0 0 0 0
                        1 1 0 0 0
                        1 1 0 0 0])

    layer_3 = BitArray([0 0 0 0 0
                        1 1 0 0 0
                        0 0 0 0 0
                        1 1 0 0 0])

    grid = BitArray(undef, 3, 4, 5)

    grid[1, :, :] .= layer_1
    grid[2, :, :] .= layer_2
    grid[3, :, :] .= layer_3

    tile_map = TM.TileMap(objects, grid)

    Test.@test TM.get_objects_type(tile_map) == typeof(objects)
    Test.@test TM.get_objects(tile_map) == objects

    Test.@test TM.get_num_objects(tile_map) == 3
    Test.@test TM.get_height(tile_map) == 4
    Test.@test TM.get_width(tile_map) == 5

    # regular indexing (indexing without using objects)
    Test.@test tile_map[1, 2, 3] == false
    Test.@test tile_map[1:3, 2, 2] == BitArray([1, 0, 1])
    Test.@test tile_map[[1, 3], 2, 2] == BitArray([1, 1])

    # indexing using one object
    Test.@test tile_map[TM.EXAMPLE_OBJECT_2, 1, 1] == tile_map[2, 1, 1] == false
    Test.@test tile_map[TM.EXAMPLE_OBJECT_2, end, end] == tile_map[2, end, end] == false
    Test.@test tile_map[TM.EXAMPLE_OBJECT_2, :, 1] == tile_map[2, :, 1] == BitArray([0, 0, 1, 1])
    Test.@test tile_map[TM.EXAMPLE_OBJECT_2, end, :] == tile_map[2, end, :] == BitArray([1, 1, 0, 0, 0])
    Test.@test tile_map[TM.EXAMPLE_OBJECT_2, 2:3, 2:4] == tile_map[2, 2:3, 2:4] == BitArray([0 0 0
                                                                                1 0 0])

    # indexing using more than one object
    Test.@test tile_map[[TM.EXAMPLE_OBJECT_2], 1, 1] == tile_map[[2], 1, 1] == BitArray([0])
    Test.@test tile_map[[TM.EXAMPLE_OBJECT_3, TM.EXAMPLE_OBJECT_1], 2, 1] == tile_map[[3, 1], 2, 1] == BitArray([1, 0])
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Julia API functions to implement conversions between representations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Conversions to DCM
# ==============================================================================

Base.convert(::Type{<:DCM}, a::EulerAngles)    = angle_to_dcm(a)
Base.convert(::Type{<:DCM}, a::Quaternion)     = quat_to_dcm(a)
Base.convert(::Type{<:DCM}, a::EulerAngleAxis) = angleaxis_to_dcm(a)

# Conversion to Euler angles
# ==============================================================================

function Base.convert(::Type{<:_EulerAngleConversion{R}}, a::DCM) where R
    return dcm_to_angle(a, R)
end

function Base.convert(::Type{<:_EulerAngleConversion{R}}, a::EulerAngles) where R
    return angle_to_angle(a, R)
end

function Base.convert(::Type{<:_EulerAngleConversion{R}}, a::EulerAngleAxis) where R
    return angleaxis_to_angle(a, R)
end

function Base.convert(::Type{<:_EulerAngleConversion{R}}, a::Quaternion) where R
    return quat_to_angle(a, R)
end

Base.convert(::Type{<:EulerAngles}, a::DCM)            = dcm_to_angle(a, :ZYX)
Base.convert(::Type{<:EulerAngles}, a::EulerAngleAxis) = angleaxis_to_angle(a, :ZYX)
Base.convert(::Type{<:EulerAngles}, a::Quaternion)     = quat_to_angle(a, :ZYX)

# Conversions to Euler angle and axis
# ==============================================================================

Base.convert(::Type{<:EulerAngleAxis}, a::DCM)         = dcm_to_angleaxis(a)
Base.convert(::Type{<:EulerAngleAxis}, a::EulerAngles) = angle_to_angleaxis(a)
Base.convert(::Type{<:EulerAngleAxis}, a::Quaternion)  = quat_to_angleaxis(a)

# Conversions to quaternions
# ==============================================================================

Base.convert(::Type{<:Quaternion}, a::DCM)            = dcm_to_quat(a)
Base.convert(::Type{<:Quaternion}, a::EulerAngles)    = angle_to_quat(a)
Base.convert(::Type{<:Quaternion}, a::EulerAngleAxis) = angleaxis_to_quat(a)


text:
function parse_input(file)
    inside_block = false
    block_type = Symbol()
    global_p = Dict{Symbol,Any}()
    current_p = Dict{Symbol,Any}()
    prev_line = ""

    for line in eachline(file)
        l = strip(line)
        @debug "Parsing line $l"
        startswith(l, '#') || isempty(l) && continue

        if !inside_block
            directive = strip.(split(l, ':'))
            @assert length(directive) == 2 "Unknown input directive $directive on line $line"

            if directive[1] == "begin"
                inside_block = true
                block_type = Symbol(directive[2])
            elseif directive[1] == "import"
                parse_input(directive[2])
            else
                @error "Unknown input command $(directive[1])"
            end
        else
            if startswith(l, "end")
                ending = strip(split(l, ':')[2])
                @assert ending == string(block_type) "Invalid block ending on line $l expected $block_type got $ending"
                inside_block = false
                # handle repeted blocks with same type
                if block_type == :species
                    block_name = Symbol(string(block_type) * "_" * get(current_p, :name, "invalid"))
                    global_p = push!!(global_p, block_name=>current_p)
                else
                    global_p = push!!(global_p, block_type=>current_p)
                end
                current_p = Dict{Symbol,Any}()
                block_type = Symbol()
                continue
            end
            # trim comments
            ln = strip(split(l, '#')[1])
            isempty(ln) && continue

            if endswith(ln, "\\")
                prev_line = ln
                continue
            end
            if !isempty(prev_line)
                ln *= prev_line
                prev_line = ""
            end

            @debug "Splitting line $ln"
            if occursin('=', ln)
                key, value = split(ln, '=')
            elseif occursin(':', ln)
                key, value = split(ln, ':')
            else
                @error "Unable to split key value pair on line $ln"
            end
            key = strip(key)
            value = strip(value)

            # replace constants
            for (i, j) in input_deck_constants
                value = replace(value, i=>j)
            end
            @debug "Got key-value pair: $key = $value"
            global_p, current_p = process_value(
                Val(block_type),
                key,
                value,
                global_p, current_p)
        end
    end

    all_species = ()
    @debug "Adding species list"
    for k in keys(global_p)
        str = string(k)
        if startswith(str, "species")
            @debug "Processing species $str"
            all_species = push!!(all_species, split(str, '_', limit=2)[2])
        end
    end
    push!!(global_p, :species=>all_species)

    return global_p
end

function process_value(key, val, global_p, current_p, block_parameters)
    k = Symbol(key)
    v = parse_value(k, val, merge!!(global_p, current_p))
    if haskey(current_p, k)
        @debug "Key $k already exists in $current_p. Replacing with $v."
        current_p = delete!!(current_p, k)
    end
    current_p = push!!(current_p, k=>v)
    if is_global(key, block_parameters)
        global_p = push!!(global_p, k=>v)
    end
    return global_p, current_p
end

function is_global(key, block_p)
    any(occursin.(block_p, (key,)))
end

function parse_value(k, v, existing_params)
    str, val_unit = replace_existing_params(v, existing_params)
    @debug "Trying to parse units in $str. Units from replacement: $val_unit"
    val = try
        uparse(str)
    catch err
        @debug "Could not parse $str. Error: $err"
        str
    end
    @debug "Tried to parse $k with unitful and got $val"
    # try to auto-add units
    if !isa(val, String) &&
            k in keys(input_unitful_entries) &&
            unit(val) == NoUnits &&
            val_unit == NoUnits

        val *= input_unitful_entries[k]
        @debug "Added unit: $val"
    end

    if val isa String
        @debug "Could not fully parse $val."
        v
    else
        @debug "Successfully parsed $val. Added additional units $val_unit"
        val * val_unit
    end
end

function replace_existing_params(str, existing_params)
    # simple replacements with previously defined keys
    @debug "Looking for known values in $str"
    new_str = deepcopy(str)
    val_units = NoUnits
    for m in eachmatch(r"\w+\d?", str)
        val = m.match
        if val in string.(keys(existing_params))
            known_val = getindex(existing_params, Symbol(val))
            if known_val isa String
                continue
            end
            val_units = unit(known_val)
            @debug "Replacing $val with $known_val"
            # Workaround https://github.com/PainterQubits/Unitful.jl/issues/391
            if val_units ≠ NoUnits
                new_val = string(ustrip(known_val))
            else
                new_val = string(known_val)
            end
            new_str = replace(new_str, r"\b"*val*r"\b"=>"($new_val)")
            @debug "Updated string with new value: $new_str"
        end
    end
    new_str, val_units
end

const input_deck_constants = Dict(
    r"\bT\b" => "true",
    r"\bF\b" => "false",
    r"\bqe\b" => "q",
    r"\bepsilon0\b" => "ϵ0",
    r"\bmu0\b" => "μ0",
    r"\bev\b" => "eV",
    r"\bkev\b" => "keV",
    r"\bmev\b" => "MeV",
    r"\bmicron\b" => "μm",
    r"\bmilli\b" => "1e-3",
    r"\bmicro\b" => "1e-6",
    r"\bnano\b" => "1e-9",
    r"\bpico\b" => "1e-12",
    r"\bfemto\b" => "1e-15",
    r"\batto\b" => "1e-18",
    r"\bcc\b" => "1e-6",
)

const input_unitful_entries = Dict(
    :t_end => u"s",
    :x_min => u"m",
    :y_min => u"m",
    :z_min => u"m",
    :x_max => u"m",
    :y_max => u"m",
    :z_max => u"m",
    :charge => u"q",
    :mass => m_e,
    :number_density => u"m^-3",
    :temperature => u"K",
    :temp => u"K",
    :temperature_ev => u"eV",
    :temp_ev => u"eV",
    :amp => u"V/m",
    :intensity_w_cm2 => u"W/cm^2",
    :omega => u"rad/s",
    :frequency => u"Hz",
    :lambda => u"m",
    :polarisation_angle => u"rad",
    :pol_angle => u"rad",
    :polarisation => u"°",
    :pol => u"°",
    :phase => u"rad",
    :t_start => u"s",
    :t_end => u"s",
    :dt_snapshot => u"s",
)


text:
using ProbNumDiffEq
using Test
using LinearAlgebra

@testset "destats.nf testing $alg" for q in (1, 2, 3, 5),
    init in (TaylorModeInit(), ClassicSolverInit()),
    alg in (
        EK0(order=q, smooth=false, initialization=init),
        EK1(order=q, smooth=false, initialization=init),
        EK1(order=q, smooth=false, initialization=init, autodiff=false),
        # EK1FDB(order=q, smooth=false, initialization=init, jac_quality=1),
    )

    f_counter = [0]
    function f(du, u, p, t; f_counter=f_counter)
        f_counter .+= 1
        mul!(du, p, u)
        return nothing
    end
    u0 = [1]
    p = [-1]
    tspan = (0.0, 3.0)
    prob = ODEProblem(f, u0, tspan, p)
    sol = solve(prob, alg, save_everystep=false)
    # @info alg sol.destats.nf f_counter[1]
    # @info sol.destats f_counter
    @test sol.destats.nf == f_counter[1]
end


text:
using NiLang, NiLang.AD
using SymbolicUtils
using SymbolicUtils: Term, Sym
using LinearAlgebra

const SymReal = Sym{Real}
const TermReal = Term{Real}
const SReals = Union{Term{Real}, Sym{Real}}

import NiLang: INC, DEC, ROT, IROT, FLIP
@inline FLIP(b::Sym{Bool}) = !b

@inline function INC(a!::SReals)
    a! + one(a!)
end

@inline function DEC(a!::SReals)
    a! - one(a!)
end

@inline function ROT(i::SReals, j::SReals, θ::SReals)
    a, b = rot(i, j, θ)
    a, b, θ
end

@inline function IROT(i::SReals, j::SReals, θ::SReals)
    i, j, _ = ROT(i, j, -θ)
    i, j, θ
end

NiLang.AD.GVar(x::SReals) = NiLang.AD.GVar(x, zero(x))
Base.convert(::Type{SymReal}, x::Integer) = SymReal(Symbol(x))
Base.convert(::Type{Term{Real}}, x::Integer) = TermReal(Symbol(x))

Base.zero(x::Sym{T}) where T = zero(Sym{T})
Base.one(x::Sym{T}) where T = one(Sym{T})
Base.zero(::Type{<:Sym{T}}) where T = Sym{T}(Symbol(0))
Base.zero(::Type{<:Term{T}}) where T = Term{T}(Symbol(0))
Base.one(::Type{<:Sym{T}}) where T = Sym{T}(Symbol(1))
Base.one(::Type{<:Term{T}}) where T = Term{T}(Symbol(1))
Base.iszero(x::Sym{T}) where T = x === zero(x)
Base.adjoint(x::SReals) = x
SymbolicUtils.Term{T}(x::Sym{T}) where T = Term{T}(x.name)

LinearAlgebra.dot(a::T, b::T) where T<:SReals = a * b

include("sparse.jl")

using BenchmarkTools, Random
syms = @syms a::Real b::Real c::Real d::Real e::Real f::Real g::Real
Base.rand(r::Random.AbstractRNG, ::Type{SymReal}, i::Integer) = rand(r, syms, i)
Base.rand(r::Random.AbstractRNG, ::Type{TermReal}, i::Integer) = rand(r, TermReal.(syms), i)
a = sprand(TermReal, 100, 100, 0.05);
b = sprand(TermReal, 100, 100, 0.05);
@benchmark SparseArrays.dot($a, $b)
@benchmark idot(TermReal(Symbol(0)), $a, $b)
@benchmark Grad(idot)(Val(1), TermReal(Symbol(0)), $a, $b)
GVar(a)

include("Symbolics/symlib.jl")
syms = @vars a b c d e f g
Base.rand(r::Random.AbstractRNG, ::Type{<:Basic}, i::Integer) = rand(r, syms, i)
a = sprand(Basic, 100, 100, 0.05);
b = sprand(Basic, 100, 100, 0.05);
@benchmark SparseArrays.dot($a, $b)
@benchmark idot(Basic(0), $a, $b)
@benchmark Grad(idot)(Val(1), Basic(0), $a, $b)


text:
"""
Gallery of `SeisInvExpt`.
# Arguments
* `attrib::Symbol` : name
  * `=:pizza` all around sources and receivers
  * `=:downhole` sources and receivers along a drill-bit
"""
function SeisInvExpt(attrib::Symbol; born_flag=false, rfields=[:Vx,:Vz,:P], α=0.0, parameterization=[:χvp, :χρ, :null])

	if(attrib==:pizza)
		# starting model
		model0 = Gallery.Seismic(:acou_homo2);
		model = deepcopy(model0)
		# add some noise to starting model
		Models.Seismic_addon!(model0, randn_perc=0.1, fields=[:χvp,:χρ])

		print(model)

		# add perturbations
		for ellip_loc in [[500.,0.], [0.,500.], [-500.,0.], [0.,-500.]]
			Models.Seismic_addon!(model, ellip_rad=50., ellip_loc=ellip_loc,
				ellip_pert=0.1, fields=[:χvp,:χρ])
		end
		Models.Seismic_addon!(model, randn_perc=0.1, fields=[:χvp,:χρ])

		# sources, receivers
		acqgeom=Acquisition.Geom_circ(nss=5,nr=20,rad=[900.,900.])

		acqsrc=Acquisition.Src_fixed_mod(acqgeom.nss,1,[:P],mod=model,nλ=3)
		tgrid=acqsrc.tgrid
		igrid=broadcast(x->range(x[1],stop=x[end],step=50.),model.mgrid)
		igrid_interp_scheme=:B2
	elseif(attrib==:downhole)
		gri = [range(-120.,stop=30.,step=1.0), range(-30.,stop=30.,step=1.0)]
		model = Models.Seismic_zeros(gri)
		Models.adjust_bounds!(model, [2500., 3500.], [2500., 3500.], [2500.,3500.])
		model0=deepcopy(model)
		print(model)
		Models.Seismic_addon!(model, randn_perc=0.1, fields=[:χvp,:χρ])
		Models.Seismic_addon!(model0, randn_perc=0.1, fields=[:χvp,:χρ])
		Models.Seismic_addon!(model, ellip_rad=[2000., 5.], ellip_loc=[0.,0.],ellip_pert=0.1, ellip_α=α, fields=[:χvp],)

		acqgeom=Acquisition.Geom_fixed(-75., -50., 0.0, -100., -50., 0.0, 2, 20, :vertical, :vertical)
		acqsrc=Acquisition.Src_fixed_mod(acqgeom.nss,1,[:P],mod=model,nλ=8, tmaxfrac=0.8)
		tgrid=acqsrc.tgrid
		igrid=[range(-40.,stop=30.,length=100), range(-30.,length=1,step=60.)]
		igrid_interp_scheme=:B2
	else
		error("invalid attrib")
	end

	if(born_flag)
		pa = FWI.Param(acqsrc, acqgeom, tgrid, FWI.ModFdtdBorn(),
		       rfields=rfields,
		       model0, modm0=model0, modm_obs=model,
		       igrid_interp_scheme=igrid_interp_scheme,
		       igrid=igrid, parameterization=parameterization, verbose=false);

	else
		pa = FWI.Param(acqsrc, acqgeom, tgrid, FWI.ModFdtd(),
		       model0, modm_obs=model, modm0=model0,
		       rfields=rfields,
		       igrid_interp_scheme=igrid_interp_scheme,
		       igrid=igrid, parameterization=parameterization, verbose=false);
	end

	return pa
end


text:
export dualize

# MOI dualize
function dualize(
    primal_model::MOI.ModelLike;
    dual_names::DualNames = EMPTY_DUAL_NAMES,
    variable_parameters::Vector{VI} = VI[],
    ignore_objective::Bool = false,
)
    # Creates an empty dual problem
    dual_problem = DualProblem{Float64}()
    return dualize(
        primal_model,
        dual_problem,
        dual_names,
        variable_parameters,
        ignore_objective,
    )
end

function dualize(
    primal_model::MOI.ModelLike,
    dual_problem::DualProblem{T};
    dual_names::DualNames = EMPTY_DUAL_NAMES,
    variable_parameters::Vector{VI} = VI[],
    ignore_objective::Bool = false,
) where {T}
    # Dualize with the optimizer already attached
    return dualize(
        primal_model,
        dual_problem,
        dual_names,
        variable_parameters,
        ignore_objective,
    )
end

function dualize(
    primal_model::MOI.ModelLike,
    dual_problem::DualProblem{T},
    dual_names::DualNames,
    variable_parameters::Vector{VI},
    ignore_objective::Bool,
) where {T}
    # Throws an error if objective function cannot be dualized
    supported_objective(primal_model)

    # Query all constraint types of the model
    con_types = MOI.get(primal_model, MOI.ListOfConstraints())
    supported_constraints(con_types) # Throws an error if constraint cannot be dualized

    # Set the dual model objective sense
    set_dual_model_sense(dual_problem.dual_model, primal_model)

    # Get Primal Objective Coefficients
    primal_objective =
        get_primal_objective(primal_model, variable_parameters, T)

    # Add variables to the dual model and their dual cone constraint.
    # Return a dictionary from dual variables to primal constraints constants (obj coef of dual var)
    dual_obj_affine_terms = add_dual_vars_in_dual_cones(
        dual_problem.dual_model,
        primal_model,
        dual_problem.primal_dual_map,
        dual_names,
        con_types,
    )

    add_primal_parameter_vars(
        dual_problem.dual_model,
        primal_model,
        dual_problem.primal_dual_map,
        dual_names,
        variable_parameters,
        primal_objective,
        ignore_objective,
    )
    add_quadratic_slack_vars(
        dual_problem.dual_model,
        primal_model,
        dual_problem.primal_dual_map,
        dual_names,
        variable_parameters,
        primal_objective,
        ignore_objective,
    )

    # Add dual equality constraint and get the link dictionary
    scalar_affine_terms = add_dual_equality_constraints(
        dual_problem.dual_model,
        primal_model,
        dual_problem.primal_dual_map,
        dual_names,
        primal_objective,
        con_types,
        variable_parameters,
    )

    if ignore_objective
        # do not add objective
    else
        # Fill Dual Objective Coefficients Struct
        dual_objective = get_dual_objective(
            dual_problem,
            dual_obj_affine_terms,
            primal_objective,
            con_types,
            scalar_affine_terms,
            variable_parameters,
        )
        # Add dual objective to the model
        set_dual_objective(dual_problem.dual_model, dual_objective)
    end
    return dual_problem
end

# JuMP dualize
function dualize(model::JuMP.Model; dual_names::DualNames = EMPTY_DUAL_NAMES)
    # Create an empty JuMP model
    JuMP_model = JuMP.Model()

    if JuMP.mode(model) != JuMP.AUTOMATIC # Only works in AUTOMATIC mode
        error(
            "Dualization does not support solvers in $(model.moi_backend.mode) mode",
        )
    end
    # Dualize and attach to the model
    dualize(
        backend(model),
        DualProblem(backend(JuMP_model));
        dual_names = dual_names,
    )
    fill_obj_dict_with_variables!(JuMP_model)
    return JuMP_model
end

function dualize(
    model::JuMP.Model,
    optimizer_constructor;
    dual_names::DualNames = EMPTY_DUAL_NAMES,
)
    # Dualize the JuMP model
    dual_JuMP_model = dualize(model; dual_names = dual_names)
    # Set the optimizer
    JuMP.set_optimizer(dual_JuMP_model, optimizer_constructor)
    return dual_JuMP_model
end

# dualize docs
"""
    dualize(args...; kwargs...)

The `dualize` function works in three different ways. The user can provide:

* A `MathOptInterface.ModelLike`

The function will return a `DualProblem` struct that has the dualized model
and `PrimalDualMap{Float64}` for users to identify the links between primal
and dual model. The `PrimalDualMap{Float64}` maps variables and constraints
from the original primal model into the respective objects of the dual model.

* A `MathOptInterface.ModelLike` and a `DualProblem{T}`

* A `JuMP.Model`

The function will return a JuMP model with the dual representation of the problem.

* A `JuMP.Model` and an optimizer constructor

The function will return a JuMP model with the dual representation of the problem with
the optimizer constructor attached.

On each of these methods, the user can provide the following keyword arguments:

* `dual_names`: of type `DualNames` struct. It allows users to set more intuitive names
for the dual variables and dual constraints created.

* `variable_parameters`: A vector of MOI.VariableIndex containing the variables that
should not be considered model variables during dualization. These variables will behave
like constants during dualization. This is specially useful for the case of bi-level
modelling, where the second level depends on some decisions from the upper level.

* `ignore_objective`: a boolean indicating if the objective function should be
added to the dual model. This is also useful for bi-level modelling, where the second
level model is represented as a KKT in the upper level model.

"""
function dualize end

function fill_obj_dict_with_variables!(model::JuMP.Model)
    all_indices = MOI.get(
        model,
        JuMP.MOI.ListOfVariableIndices(),
    )::Vector{MOI.VariableIndex}
    for vi in all_indices
        model.obj_dict[Symbol(
            MOI.get(backend(model), MOI.VariableName(), vi),
        )] = VariableRef(model, vi)
    end
    return model
end


text:
using Statistics
using LinearAlgebra
using SparseArrays

include("../src/fdfd_maxwell_solver.jl")
include("../src/geometry_code.jl")
include("../src/simulation_unit_cell.jl")
include("../src/misc.jl")


text:
using jInv.Mesh;
using ForwardHelmholtz
#import EikonalInv.expandModelNearest
#import EikonalInv.addAbsorbingLayer
plotting = false;
if plotting
	using PyPlot;
	close("all")
end

include("My_sub2ind.jl");

# m = readdlm("SEGmodel2Dsalt.dat"); m = m';
m = ones(256,128);
m = m*1e-3;
Minv = getRegularMesh([0.0,13.5,0.0,4.2],collect(size(m))-1);

pad = 25;

m = ones(size(m));
m = 1./m.^2

if plotting
	figure();
	imshow(1./sqrt(m'))
end
f = 1.5;
w = 2*pi*f

maxOmega = getMaximalFrequency(m,Minv);
ABLamp = maxOmega;

println("omega*h:");
println(w*Minv.h*sqrt(maximum(m)));
pad = pad*ones(Int64,Minv.dim);
H = GetHelmholtzOperator(Minv,m,w,ones(size(m))*0.01,true,pad,ABLamp,true)[1];
SH = H + GetHelmholtzShiftOP(m, w,0.1);
n = Minv.n+1; n_tup = tuple(n...);
src = div.(n,2);
src[end] = 1;
q = zeros(Complex,n_tup)
q[My_sub2ind(n,src)] = 1/(Minv.h[1]^2);

s = H\q[:];
s = real(reshape(s,n_tup));
s[My_sub2ind(n,src)] = 0.0;

if plotting
	figure();
	imshow(s');title("Helmholtz Solution");
end

sh = SH\q[:];
sh = real(reshape(sh,n_tup));
sh[My_sub2ind(n,src)]= 0.0;

if plotting
	figure();
	imshow(sh');title("Shifted Helmholtz Solution");
end


text:
"""
    as_dataframe(clicks::Vector{T}, starttime::DateTime, wavfile) where {T<:AbstractClick}
    as_dataframe(clicks::Vector{T}, starttime::DateTime) where {T<:AbstractClick}
    as_dataframe(clicks::Vector{T}) where {T<:AbstractClick}

Converts a vector of `ClickPointer` or `Click` structures to tabular form as a
`DataFrame`. If optional argument `starttime` is provided, the returned
`DataFrame` will also include `start_datetime` and `stop_datetime` columns,
with the clicks' start and end times relative to the start of the sound file
converted to `DateTime`s relative to `starttime` (i.e., the time of the first
sample in the analyzed file).  The source file for the clicks can be included
via th optional argument `wavfile` (either a `String` or `Vector{String}`).
"""
function as_dataframe(clicks::Vector{T}, starttime::DateTime) where {T<:AbstractClick}
    df = DataFrame(
        (start=left(c),
         stop=right(c),
         start_datetime=starttime + Dates.Millisecond(round(left(c) * 1e3)),
         stop_datetime=starttime + Dates.Millisecond(round(right(c) * 1e3)),
         samplerate=samplerate(c),
         waveform=samples(c)) for c in clicks)
    return df
end

function as_dataframe(clicks::Vector{T}) where {T<:AbstractClick}
    df = DataFrame((start=left(c),
                    stop=right(c),
                    samplerate=samplerate(c),
                    waveform=samples(c)) for c in clicks)
    return df
end

function as_dataframe(clicks::Vector{T}, starttime::DateTime, wavfile) where {T<:AbstractClick}
    df = as_dataframe(clicks, starttime)
    df[!, :wavfile] .= wavfile
    return df
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: iot_1click_devices_service
using AWS.Compat
using AWS.UUIDs

"""
    claim_devices_by_claim_code(claim_code)
    claim_devices_by_claim_code(claim_code, params::Dict{String,<:Any})

Adds device(s) to your account (i.e., claim one or more devices) if and only if you
received a claim code with the device(s).

# Arguments
- `claim_code`: The claim code, starting with \"C-\", as provided by the device
  manufacturer.

"""
claim_devices_by_claim_code(claimCode; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("PUT", "/claims/$(claimCode)"; aws_config=aws_config)
claim_devices_by_claim_code(claimCode, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("PUT", "/claims/$(claimCode)", params; aws_config=aws_config)

"""
    describe_device(device_id)
    describe_device(device_id, params::Dict{String,<:Any})

Given a device ID, returns a DescribeDeviceResponse object describing the
 details of the
device.

# Arguments
- `device_id`: The unique identifier of the device.

"""
describe_device(deviceId; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("GET", "/devices/$(deviceId)"; aws_config=aws_config)
describe_device(deviceId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("GET", "/devices/$(deviceId)", params; aws_config=aws_config)

"""
    finalize_device_claim(device_id)
    finalize_device_claim(device_id, params::Dict{String,<:Any})

Given a device ID, finalizes the claim request for the associated device.
 Claiming a
device consists of initiating a claim, then publishing a device event,
 and finalizing the
claim. For a device of type button, a device event can
 be published by simply clicking the
device.

# Arguments
- `device_id`: The unique identifier of the device.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"tags"`: A collection of key/value pairs defining the resource tags. For example, {
  \"tags\": {\"key1\": \"value1\", \"key2\": \"value2\"} }. For more information, see AWS
  Tagging Strategies.
"""
finalize_device_claim(deviceId; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("PUT", "/devices/$(deviceId)/finalize-claim"; aws_config=aws_config)
finalize_device_claim(deviceId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("PUT", "/devices/$(deviceId)/finalize-claim", params; aws_config=aws_config)

"""
    get_device_methods(device_id)
    get_device_methods(device_id, params::Dict{String,<:Any})

Given a device ID, returns the invokable methods associated with the device.

# Arguments
- `device_id`: The unique identifier of the device.

"""
get_device_methods(deviceId; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("GET", "/devices/$(deviceId)/methods"; aws_config=aws_config)
get_device_methods(deviceId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("GET", "/devices/$(deviceId)/methods", params; aws_config=aws_config)

"""
    initiate_device_claim(device_id)
    initiate_device_claim(device_id, params::Dict{String,<:Any})

Given a device ID, initiates a claim request for the associated device.
 Claiming a device
consists of initiating a claim, then publishing a device event,
 and finalizing the claim.
For a device of type button, a device event can
 be published by simply clicking the
device.

# Arguments
- `device_id`: The unique identifier of the device.

"""
initiate_device_claim(deviceId; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("PUT", "/devices/$(deviceId)/initiate-claim"; aws_config=aws_config)
initiate_device_claim(deviceId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("PUT", "/devices/$(deviceId)/initiate-claim", params; aws_config=aws_config)

"""
    invoke_device_method(device_id)
    invoke_device_method(device_id, params::Dict{String,<:Any})

Given a device ID, issues a request to invoke a named device method (with possible
parameters). See the \"Example POST\" code snippet below.

# Arguments
- `device_id`: The unique identifier of the device.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"deviceMethod"`: The device method to invoke.
- `"deviceMethodParameters"`: A JSON encoded string containing the device method request
  parameters.
"""
invoke_device_method(deviceId; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("POST", "/devices/$(deviceId)/methods"; aws_config=aws_config)
invoke_device_method(deviceId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("POST", "/devices/$(deviceId)/methods", params; aws_config=aws_config)

"""
    list_device_events(device_id, from_time_stamp, to_time_stamp)
    list_device_events(device_id, from_time_stamp, to_time_stamp, params::Dict{String,<:Any})

Using a device ID, returns a DeviceEventsResponse object containing an
 array of events for
the device.

# Arguments
- `device_id`: The unique identifier of the device.
- `from_time_stamp`: The start date for the device event query, in ISO8061 format. For
  example,
 2018-03-28T15:45:12.880Z
- `to_time_stamp`: The end date for the device event query, in ISO8061 format. For
  example,
 2018-03-28T15:45:12.880Z

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to return per request. If not set, a
  default value of
 100 is used.
- `"nextToken"`: The token to retrieve the next set of results.
"""
list_device_events(deviceId, fromTimeStamp, toTimeStamp; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("GET", "/devices/$(deviceId)/events", Dict{String, Any}("fromTimeStamp"=>fromTimeStamp, "toTimeStamp"=>toTimeStamp); aws_config=aws_config)
list_device_events(deviceId, fromTimeStamp, toTimeStamp, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("GET", "/devices/$(deviceId)/events", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("fromTimeStamp"=>fromTimeStamp, "toTimeStamp"=>toTimeStamp), params)); aws_config=aws_config)

"""
    list_devices()
    list_devices(params::Dict{String,<:Any})

Lists the 1-Click compatible devices associated with your AWS account.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"deviceType"`: The type of the device, such as \"button\".
- `"maxResults"`: The maximum number of results to return per request. If not set, a
  default value of
 100 is used.
- `"nextToken"`: The token to retrieve the next set of results.
"""
list_devices(; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("GET", "/devices"; aws_config=aws_config)
list_devices(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("GET", "/devices", params; aws_config=aws_config)

"""
    list_tags_for_resource(resource-arn)
    list_tags_for_resource(resource-arn, params::Dict{String,<:Any})

Lists the tags associated with the specified resource ARN.

# Arguments
- `resource-arn`: The ARN of the resource.

"""
list_tags_for_resource(resource_arn; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("GET", "/tags/$(resource-arn)"; aws_config=aws_config)
list_tags_for_resource(resource_arn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("GET", "/tags/$(resource-arn)", params; aws_config=aws_config)

"""
    tag_resource(resource-arn, tags)
    tag_resource(resource-arn, tags, params::Dict{String,<:Any})

Adds or updates the tags associated with the resource ARN. See AWS IoT 1-Click Service
Limits for the maximum number of tags allowed per
 resource.

# Arguments
- `resource-arn`: The ARN of the resource.
- `tags`: A collection of key/value pairs defining the resource tags. For example, {
  \"tags\": {\"key1\": \"value1\", \"key2\": \"value2\"} }. For more information, see AWS
  Tagging Strategies.

"""
tag_resource(resource_arn, tags; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("POST", "/tags/$(resource-arn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resource_arn, tags, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("POST", "/tags/$(resource-arn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), params)); aws_config=aws_config)

"""
    unclaim_device(device_id)
    unclaim_device(device_id, params::Dict{String,<:Any})

Disassociates a device from your AWS account using its device ID.

# Arguments
- `device_id`: The unique identifier of the device.

"""
unclaim_device(deviceId; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("PUT", "/devices/$(deviceId)/unclaim"; aws_config=aws_config)
unclaim_device(deviceId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("PUT", "/devices/$(deviceId)/unclaim", params; aws_config=aws_config)

"""
    untag_resource(resource-arn, tag_keys)
    untag_resource(resource-arn, tag_keys, params::Dict{String,<:Any})

Using tag keys, deletes the tags (key/value pairs) associated with the specified
 resource
ARN.

# Arguments
- `resource-arn`: The ARN of the resource.
- `tag_keys`: A collections of tag keys. For example, {\"key1\",\"key2\"}

"""
untag_resource(resource_arn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("DELETE", "/tags/$(resource-arn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resource_arn, tagKeys, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("DELETE", "/tags/$(resource-arn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), params)); aws_config=aws_config)

"""
    update_device_state(device_id)
    update_device_state(device_id, params::Dict{String,<:Any})

Using a Boolean value (true or false), this operation
 enables or disables the device given
a device ID.

# Arguments
- `device_id`: The unique identifier of the device.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"enabled"`: If true, the device is enabled. If false, the device is
 disabled.
"""
update_device_state(deviceId; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("PUT", "/devices/$(deviceId)/state"; aws_config=aws_config)
update_device_state(deviceId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_1click_devices_service("PUT", "/devices/$(deviceId)/state", params; aws_config=aws_config)


text:
module IntrinsicMetrics

using NearestNeighbors
using StaticArrays
using Distances
using LightGraphs, SimpleWeightedGraphs
using RecipesBase

include("intrinsic_metric.jl")
include("plotting.jl")

export IntrinsicMetric,
    ambientmetric, adjgraph, points, npoints, ambientdim, points3d,
    evaluate, pairwise!, pairwise, colwise!, colwise,
    Euclidean, Minkowski, Cityblock, Chebyshev

end


text:
# entire file copied from https://www.cs.purdue.edu/homes/dgleich/codes/hyperkron/hyperkron.jl
"""
authors:
David F. Gleich
Arjun S. Ramani
Nicole Eikmeier
Charun Upara (update to Julia 1.4)
Joshua Turner (update to Julia 1.4)
"""

using SparseArrays
using LinearAlgebra
using Random
"""
`kronecker_graph`
=================

Create an instance of a Kronecker graph using the grass-hopping algorithm.
"""
## Regions

# Increment to the next region of a Kronecker graph
@inline function _next_region!(cur::AbstractArray{T,1}, m::Integer) where T <: Integer
  lastk = length(cur)
  for k=length(cur):-1:1
    lastk = k              # save the last state
    cur[k] += 1            # increment the last element
    if cur[k] == m+1       # if there is spill
      if k > 1             # if there is an array
        continue # recur on prefix
      else
        cur[k] = 0
      end
    end
    break
  end
  for k=lastk+1:length(cur)
    cur[k] = cur[k-1]
  end
  cur
end

#=
cur = [1,1,1,1]
@show _next_region!(cur, 3)
@show _next_region!(cur, 3)
@show _next_region!(cur, 3)
@show _next_region!(cur, 3)
@show _next_region!(cur, 3)
@time _next_region!(cur, 3)

cur = [1,1,1]
@show _next_region!(cur, 4)
while cur[1] != 0
  @show _next_region!(cur, 4)
end
=#

function _first_region(a::AbstractArray,k::Integer)
  return ones(UInt32,k)
end

function _all_regions(a::AbstractArray,k::Integer)
  # CU: no need to call this function here, can just call ones() directly?
  cur = _first_region(a,k)
  # CU: the type seems to be hard-coded, so I'm setting it to UInt32 to make things easier
  rval = Array{UInt32}[]
  while cur[1] != 0
    push!(rval, copy(cur))
    _next_region!(cur, 4)
  end
  return rval
end

@show R =_all_regions([0.99,0.5,0.5,0.2],3)
@show typeof(R)
@time _all_regions([0.99,0.5,0.5,0.2],3)
@time _all_regions([0.99,0.5,0.5,0.2],4)
@time _all_regions([0.99,0.5,0.5,0.2],5)

# Unranking

# unrank([1,2,2,4], 1) returns [1,2,2,4]
# unrank([1,2,2,4], 2) returns [1,2,4,2]
# unrank([0,1,1,3], 3) returns [1,4,2,2]

# This assumes ms is sorted

function num_multiset_permutations(ms::AbstractArray{T,1}, start_pos::Integer = 1) where T<:Integer
  nperm = 1
  count = 1
  countfac = 1
  @inbounds for i in start_pos+1:length(ms)
    nperm *= (i-start_pos+1)
    if ms[i] == ms[i-1]
      count += 1
      countfac *= count
    else
      nperm = div(nperm,countfac)
      # new element
      count = 1
      countfac = 1
    end
  end
  nperm = div(nperm,countfac)
end

# Find the next distinct character after a[k]
# Returns the index where it occurs or length(a)+1
@inline function _next_distinct_character(a, k::Integer, m::Integer)
  # check for easy out
  if k>=m
    return m+1
  else
    ak = a[k]
    if ak != a[k+1] # many cases
      return k+1
    elseif ak == a[m] # some cases
      return m+1
    else
      # we have to search, so there are probably many dups
      # since we assume sorted, we can bisect
      lo = 1
      hi = m-k # we need to search [k+lo:k+hi]
      while hi - lo > 2
        mid = (hi + lo) >> 1
        if ak == a[k+mid]
          lo = mid
        else
          hi = mid
        end
      end
      ak != a[k+lo+1] && return k+lo+1
      return k+lo+2
    end
  end
end

#=
@show _next_distinct_character([5,6,6,6,8],2,4)
@show _next_distinct_character([5,6,6,6,6,8],1,6)
@show _next_distinct_character([5,6,6,6,6,8],6,6)-6
=#
# This version doesn't use recursion, and unranks the sorted set ms in place
# It also gets the multiplicity of each element to update nperm without
# recomputing it from scratch, it also can take in the total number
# of permutations and just update as we go.
#
# This version doesn't use recursion, and unranks the sorted set ms in place
# It also gets the multiplicity of each element to update nperm without
# recomputing it from scratch, it also can take in the total number
# of permutations and just update as we go.
# This version tries to avoid integer division by storing nperm as a float
@inline function unrank5!(ms::AbstractArray{T,1}, n::Integer, nperm, m) where T<:Integer
  multcur = 0
  imultcur = 1.0
  multnext = 0
  npermf = float(nperm)
  @inbounds for k=1:m # for each location in the output
    if n == 1
      break
    end

    if multcur == 0 # update
      multcur = _next_distinct_character(ms, k, m) - k # multiplicity of the current element
      imultcur = 1.0/multcur
    end

    # update nperm for removing multcur
    npermf *= multcur
    npermf /= (m-k+1) # shrink the length
    nperm = round(typeof(nperm),npermf)

    #=
    if nperm != num_multiset_permutations(ms, k+1)
      @show 0, ms, k, num_multiset_permutations(ms, k+1), nperm
    end
    @assert nperm == num_multiset_permutations(ms, k+1)
    =#

    if n <= nperm # this means we start with ms[k]
      multcur -= 1
      if multcur != 0
        imultcur = 1.0/multcur
      end
      continue
    else # we need to search
      n -= nperm

      if k+multcur <= m
        if ms[k] == ms[k+multcur]
          #@show ms, k, _next_distinct_character(ms,k,m), multcur
          @assert ms[k] != ms[k+multcur]
        end
      end

      pos = k+multcur
      while pos <= m
        # we are going to swap ms[k] and ms[pos]
        # but to update, we need to know the multiplicity of ms[pos]

        multnext = _next_distinct_character(ms,pos,m) - pos

        # so the element ms[pos] has multiplicity multnext
        # and element ms[k] has multiplicity multcur-1

        ms[pos], ms[k] = ms[k], ms[pos] # swap leading positions
        npermf *= multnext
        npermf *= imultcur
        nperm = round(typeof(nperm),npermf)

        #=
        if nperm != num_multiset_permutations(ms, k+1)
          @show 1, ms, k, num_multiset_permutations(ms, k+1), nperm
        end
        @assert nperm == num_multiset_permutations(ms, k+1)
        =#

        if n <= nperm
          break
        else
          n -= nperm
          pos += multnext
          multcur = multnext
          imultcur = 1.0/multcur
        end
      end

      multcur = 0 # reset the search
    end
  end
  return ms
end
unrank5!(ms, n) = unrank5!(ms, n, num_multiset_permutations(ms), length(ms))

@inline function mydivmod(val,n)
  if n==2
    d2 = val >> 1
    d1 = val - d2*n
  elseif n==4
    d2 = val >> 2
    d1 = val - d2*n
  else
    d1 = mod(val,n)
    d2 = div(val,n)
  end
  return d2, d1
end
@show mydivmod(5,2), divrem(5,2)

@inline function _morton_decode3(mind::AbstractArray{T,1}, n, sz) where T<:Integer
  row = 0
  col = 0
  tub = 0
  base = 1
  for val in mind
    # CU: ind2sub() was deprecated in Julia 0.7, changed to a new way instead
    s1,s2,s3 = Base._ind2sub(sz, val) #basically backward map where base keeps track of the block size n^
    row += base*(s1-1)
    col += base*(s2-1)
    tub += base*(s3-1)
    base *= n
  end
  return (row+1, col+1, tub+1)
end
@show _morton_decode3([8,8,8,8],2,(2,2,2))
#@time _morton_decode([4,4,4,4],2)

#_randgeo(p) = ceil(Int,-randexp() / log1p(-p)) # jsut a copy of Distributions.jl
#_randgeo(p) = ceil(Int,-randexp() / log1p(-p)) # jsut a copy of Distributions.jl
#_randgeo_true(p) = ceil(Int,-randexp() / log1p(-p)) # jsut a copy of Distributions.jl
#_randgeo(p) = _randgeo_true(max(0.01*p, p + 0.0*p*randn()))
_randgeo(p) = ceil(Int,-randexp() / log1p(-p))
function _generate_region!(
    ei::Array{Ta}, ej::Array{Ta}, ek::Array{Ta},
    r::AbstractArray{Tb,1},
    e::AbstractArray{Tb,1},
    p, small_n::Integer, sz,
    onlytris::Bool, onlysym::Bool) where {Ta, Tb <:Integer}

    N = num_multiset_permutations(r) # total area
    m = length(r)
    i = 0

    # ugh, do we need this? should make these tolerances.
    if N*p < 1.0e-12 || p < 1.0e-15
      return
    end

    gap = _randgeo(p)

    while i+gap <= N
      i += gap
      copy!(e, r)
      #@show e, i, N, m
      mult_edge = unrank5!(e, i, N, m) # generate the ith permutation
      #src, dst = _morton_decode(mult_edge, small_n)
      src,dst,tub = _morton_decode3(mult_edge, small_n, sz)
      # only add the edge if it's in the lower region of the matrix.
      # This means we are assuming a symmetric generation.
      if (onlysym && src <= dst <= tub) || !onlysym
        if (onlytris && (src != dst != tub != src)) || !onlytris
          push!(ei, src)
          push!(ej, dst)
          push!(ek, tub)
        end
      end
      gap = _randgeo(p)
    end
end

# K = initial tensor
# k = r in paper
# The rest is optional
function fast_hyperkron_edges(K,k;
  noise::Float64=0.0,onlytris::Bool=false,onlysym::Bool=true)
  n = size(K,1)
  sz = size(K)
  v = vec(K) # vectorized

  vn = v/sum(v)
  V = collect(v for _ in 1:k) # create a set of parameters for each level
  #V = collect(max.(0.0,min.(2.0,sum(v).*normalize!(vec(kron_params(vn[1],vn[2]+noise*(2*rand()-1.0),vn[4]+noise*(2*rand()-1.0),vn[8])),1))) for _ in 1:k )
  #@show V

  # Dense = adjacency matrix, every spot is filled in, even if there's a 0. Very large graph = expensive to store
  # Sparse = Store a list of non-zeroes edges. For i and j, there will be two one-dimension array

  # hold edges
  ei = zeros(Int,0)
  ej = similar(ei)
  ek = similar(ei)

  nregions = 0
  # v is the vectoriezd tensor
  region = _first_region(v, k)

  edge_in_region = similar(region) # allocated once to avoid lots of allocs
  while region[1] != 0
    nregions += 1
    # compute the probability of erdos-reyes regions
    p = 1.0
    level = 1
    for j in region
      p *= V[level][j]
      level += 1
    end
    _generate_region!(ei, ej, ek, region, edge_in_region, p, n, sz, onlytris, onlysym)
    _next_region!(region, length(v))
  end
  return ei, ej, ek
end

function kron_params(a,b,c,d)
  K = zeros(2,2,2)
  K[1,1,1] = a
  K[1,1,2] = b; K[1,2,1] = b;   K[2,1,1] = b
  K[1,2,2] = c; K[2,1,2] = c;   K[2,2,1] = c
  K[2,2,2] = d
  return K
end

##
function hyperedges_to_graph_with_triangles(n::Int, hedges)
  ei = zeros(Int64,0)
  ej = zeros(Int64,0)
  for i=1:length(hedges[1])
    push!(ei, hedges[1][i]); push!(ej, hedges[2][i]) # (i -> j)
    push!(ei, hedges[1][i]); push!(ej, hedges[3][i] ) # (i -> k)
    push!(ei, hedges[2][i]); push!(ej, hedges[3][i] ) # (j -> k)
  end
  A = sparse(ei,ej,1.0,n,n)
  A = max.(A,A') # symmetrize
  A = A - spdiagm(0 => diag(A)) #remove diagonal entries
  fill!(A.nzval, 1.0) # set all non-zeros
  return A
end

##
""" Generate a HyperKron Graph.

    hyperkron_graph(K,k)
    hyperkron_graph(K,k; motif=:triangle, onlytris::Bool=false)

The code calls "hypergraph_to_edges_with_" 'motif' in order
to generate the graph to return. Consequently, you can implement
your own motifs via this routine.
"""
function hyperkron_graph(K,k;
    motif_to_graph::Function = hyperedges_to_graph_with_triangles, kwargs...)

  hedges = fast_hyperkron_edges(K,k; kwargs...)
  return motif_to_graph(size(K,1)^k, hedges), hedges
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    ShiftedPath(path, offset)

Traverse an object with the `path` shifted by an `offset`
that can be positive or negative.
"""
struct ShiftedPath{P<:Path} <: Path
  path::P
  offset::Int
end

function traverse(object, path::ShiftedPath)
  p = traverse(object, path.path)
  n = length(p)
  o = path.offset
  s = o ≥ 0 ? o : abs((n + o) % n)
  i1 = Iterators.cycle(p)
  i2 = Iterators.drop(i1, s % n)
  Iterators.take(i2, n)
end


text:
function _convert(
        node   :: treeclassifier.NodeMeta{S},
        list   :: AbstractVector{T},
        labels :: AbstractVector{T}) where {S, T}

    if node.is_leaf
        return Leaf{T}(list[node.label], labels[node.region])
    else
        left = _convert(node.l, list, labels)
        right = _convert(node.r, list, labels)
        return Node{S, T}(node.feature, node.threshold, left, right)
    end
end

function _convertOOB(
        node   :: treeclassifier.NodeMeta{S},
        list   :: AbstractVector{T},
        labels :: AbstractVector{T},
        inds   :: AbstractVector{Int},
        oob    :: AbstractMatrix{S}) where {S, T}

    tree = _convert(node, list, labels)
    return TreeOOB{S, T}(tree, inds, oob)
end

function build_tree(
        inds                :: AbstractVector{Int},
        labels              :: AbstractVector{T},
        features            :: AbstractMatrix{S},
        n_subfeatures        = 0,
        max_depth            = -1,
        min_samples_leaf     = 1,
        min_samples_split    = 2,
        min_purity_increase  = 0.0;
        loss                 = util.entropy :: Function,
        rng                  = Random.GLOBAL_RNG) where {S, T}

    if max_depth == -1
        max_depth = typemax(Int)
    end
    if n_subfeatures == 0
        n_subfeatures = size(features, 2)
    end

    rng = mk_rng(rng)::Random.AbstractRNG
    t = treeclassifier.fit(
        X                   = features[inds, :],
        Y                   = labels[inds],
        W                   = nothing,
        loss                = loss,
        max_features        = Int(n_subfeatures),
        max_depth           = Int(max_depth),
        min_samples_leaf    = Int(min_samples_leaf),
        min_samples_split   = Int(min_samples_split),
        min_purity_increase = Float64(min_purity_increase),
        rng                 = rng)

    return _convertOOB(t.root, t.list, labels[t.labels],
        setdiff(collect(1:length(labels)), inds),
        features[setdiff(collect(1:length(labels)), inds),:])
end

function apply_tree_oob(tree::TreeOOB) where {S, T}
    return apply_tree(tree.tree, tree.oobfeat)
end

"""
Build a random forest which keep track of the oob sample
"""
function build_forest_oob(
        labels              :: AbstractVector{T},
        features            :: AbstractMatrix{S},
        n_subfeatures       = -1,
        n_trees             = 10,
        partial_sampling    = 0.7,
        max_depth           = -1,
        min_samples_leaf    = 1,
        min_samples_split   = 2,
        min_purity_increase = 0.0;
        rng                 = Random.GLOBAL_RNG) where {S, T}

    if n_trees < 1
        throw("the number of trees must be >= 1")
    end
    if !(0.0 < partial_sampling <= 1.0)
        throw("partial_sampling must be in the range (0,1]")
    end

    if n_subfeatures == -1
        n_features = size(features, 2)
        n_subfeatures = round(Int, sqrt(n_features))
    end

    t_samples = length(labels)
    n_samples = floor(Int, partial_sampling * t_samples)

    forest = Vector{TreeOOB{S, T}}(undef, n_trees)

    entropy_terms = util.compute_entropy_terms(n_samples)
    loss = (ns, n) -> util.entropy(ns, n, entropy_terms)

    if rng isa Random.AbstractRNG
        Threads.@threads for i in 1:n_trees
            inds = rand(rng, 1:t_samples, n_samples)
            forest[i] = build_tree(
                inds,
                labels,
                features,
                n_subfeatures,
                max_depth,
                min_samples_leaf,
                min_samples_split,
                min_purity_increase,
                loss = loss,
                rng = rng)
        end
    elseif rng isa Integer # each thread gets its own seeded rng
        Threads.@threads for i in 1:n_trees
            Random.seed!(rng + i)
            inds = rand(1:t_samples, n_samples)
            forest[i] = build_tree(
                inds,
                labels,
                features,
                n_subfeatures,
                max_depth,
                min_samples_leaf,
                min_samples_split,
                min_purity_increase,
                loss = loss)
        end
    else
        throw("rng must of be type Integer or Random.AbstractRNG")
    end
    return EnsembleOOB{S, T}(length(labels), features, forest)
end

"""
Prediction on out of bag sample for a random forest
"""
function apply_forest_oob(forest::EnsembleOOB{S, T}) where {S, T}
    n_samp = forest.nsamp
    n_trees = length(forest.trees)
    predictions = Array{T}(undef, n_samp)
    for i in 1:n_samp
        votes = Array{T}(undef, n_trees)
        len = 0
        obs_feat = forest.features[i, :]
        for j in 1:n_trees
            if i in forest.trees[j].inds
                len += 1
                votes[len] = apply_tree(forest.trees[j].tree, obs_feat)
            end
        end
        trunc_votes = votes[1:len]
        if T <: Float64
            predictions[i] = mean(trunc_votes)
        else
            predictions[i] = majority_vote(trunc_votes)
        end
    end
    return predictions
end


text:
using VariationalInequality
using JuMP

# https://cdr.lib.unc.edu/indexablecontent/uuid:778ca632-74ca-4858-8c3c-6dcfc7e6e703
# Example 3.8. This example is used for testing the RPM in (Fukushima, 1986).
m = VIPModel()

@variable(m, x1)
@variable(m, x2)
@variable(m, x3)

@NLconstraint(m, const1, x1^2 + 0.4x2^2 + 0.6x3^2 <= 1)

@NLexpression(m, F1, 2x1 + 0.2x1^3 - 0.5x2 + 0.1x3 - 4)
@NLexpression(m, F2, -0.5x1 + x2 + 0.1x2^3 + 0.5)
@NLexpression(m, F3, 0.5x1 - 0.2x2 + 2x3 - 0.5)

correspond(m, [F1, F2, F3], [x1, x2, x3])

sol1, Fval1, gap1 = solveVIP(m, algorithm=:fixed_point, max_iter=1000, step_size=0.1)
@assert 0<= gap1 < 1e-6

# The above `solveVIP` sets the value of variables at the solution
clearValues(m)

sol2, Fval2, gap2 = solveVIP(m, algorithm=:extra_gradient, max_iter=1000, step_size=0.1)
@assert 0<= gap2 < 1e-6

# sol2, Fval2, gap2 = solveVIP(m, algorithm=:hyperplane, max_iter=1000, step_size=0.1)

@show sol1
@show sol2
@show gap1
@show gap2

println("x1 = ", sol1[x1] )
println("x2 = ", sol1[x2] )
println("x3 = ", sol1[x3] )

@test isapprox(sol1[x1], 1; atol=1e-5)
@test isapprox(sol1[x2], 0; atol=1e-5)
@test isapprox(sol1[x3], 0; atol=1e-5)


text:
"""
    PolyaGammaPSWSampler(b::Int, z::Real)

PSW sampler ([1]) for a Polya-Gamma distribution with parameters `b` and `z`,
and Laplace transform

```math
\\mathcal{L}(t) = \\cosh^b(z) \\cosh^{-b}(\\sqrt{2t + z^2})
```

References

* [1] <https://doi.org/10.1080/01621459.2013.829001>
"""
struct PolyaGammaPSWSampler{T <: Real} <: Sampleable{Univariate, Continuous}
    b::Int
    z::T
end

function Base.rand(rng::AbstractRNG, s::PolyaGammaPSWSampler)
    out = 0.0
    for _ in 1:s.b
        out += rpg_devroye_1(rng, s.z)
    end
    return out
end

function rpg_devroye_1(rng::AbstractRNG, z::Real)
    t = 0.64
    z = abs(z) * 0.5
    yz = π^2 / 8 + z^2 / 2
    while true
        if rand(rng) < mass_texpon(z, t)
            x = t + randexp(rng) / yz
        else
            x = rtigauss(rng, z, t)
        end
        s = acoef(0, x, t)
        y = s * rand(rng)
        n = 0
        while true
            n += 1
            if n % 2 == 1
                s -= acoef(n, x, t)
                y <= s && return x / 4
            else
                s += acoef(n, x, t)
                y > s && break
            end
        end
    end
end

function mass_texpon(z::Real, t::Real)
    yz = π^2 / 8 + z^2 / 2
    b = inv(sqrt(t)) * (t * z - 1)
    a = -inv(sqrt(t)) * (t * z + 1)
    t0 = log(yz) + yz * t
    tb = t0 - z + normlogcdf(b)
    ta = t0 + z + normlogcdf(a)
    qdivp = fourinvπ * (exp(tb) + exp(ta))
    return 1.0 / (1.0 + qdivp)
end

function rtigauss(rng::AbstractRNG, z::Real, t::Real)
    z = abs(z)
    μ = inv(z)
    x = t + 1.0
    if μ > t
        α = 0.0
        while rand(rng) > α
            e1 = randexp(rng)
            e2 = randexp(rng)
            while e1^2 > (2 * e2 / t)
                e1 = randexp(rng)
                e2 = randexp(rng)
            end
            x = t / (1 + t * e1)^2
            α = exp(-0.5 * z^2 * x)
        end
    else
        while x > t
            y = randn(rng)^2
            x = μ + μ^2 * y / 2 - μ * √(4 * μ * y + (μ * y)^2) / 2
            if rand(rng) > μ / (μ + x)
                x = μ^2 / x
            end
        end
    end
    return x
end

function acoef(n::Int, x::Real, t::Real)
    if x > t
        π * (n + 0.5) * exp(-(n + 0.5)^2 * π^2 * x / 2)
    else
        2 / x / √x / sqrthalfπ * (n + 0.5) * exp(-2 * (n + 0.5)^2 / x)
    end
end

function Distributions.mean(s::PolyaGammaPSWSampler)
    s.b * inv(2.0 * s.z) * tanh(s.z / 2.0)
end

function Distributions.var(s::PolyaGammaPSWSampler)
    s.b * inv(4 * s.z^3) * (sinh(s.z) - s.z) * (sech(s.z / 2)^2)
end


text:
abstract type InitialConditionChronology end

"""
    InterProblemChronology()

    Type struct to select an information sharing model between stages that uses results from the most recent stage executed to calculate the initial conditions. This model takes into account solutions from stages defined finer resolutions
"""

struct InterProblemChronology <: InitialConditionChronology end

"""
    InterProblemChronology()

    Type struct to select an information sharing model between stages that uses results from the same recent stage to calculate the initial conditions. This model ignores solutions from stages defined finer resolutions.
"""
struct IntraProblemChronology <: InitialConditionChronology end


text:
using Base.Test
using CUDArt
using CUDNN
@show CUDNN_VERSION

import CUDNN.cudnnConvolutionForward
import CUDNN.cudnnConvolutionBackwardFilter
import CUDNN.cudnnConvolutionBackwardData
import CUDNN.cudnnGetConvolutionNdForwardOutputDim
import CUDNN.cudnnGetPoolingNdForwardOutputDim
import CUDNN.cudnnPoolingForward
import CUDNN.cudnnPoolingBackward
include(Knet.dir("src","util","conv_pool_cpu.jl"))

padding=0
stride=1
x = reshape(Float32[1:16;], 4, 4, 1, 1); tx = CudaArray(x); @show x
w = reshape(Float32[1:4;], 2, 2, 1, 1); tw = CudaArray(w); @show w
y = zeros(Float32,3,3,1,1); ty = CudaArray(y); @show y
cudnnConvolutionForward(tx,tw,ty; padding=padding, stride=stride); y = to_host(ty); @show y
y2 = zeros(Float32,3,3,1,1); @show y2
cudnnConvolutionForward(x,w,y2; padding=padding, stride=stride); @show y2
@test_approx_eq y y2

dy = rand(Float32, size(y)); tdy = CudaArray(dy); @show dy
dw = zeros(Float32, size(w)); tdw = CudaArray(dw); @show dw
cudnnConvolutionBackwardFilter(tx,tdy,tdw); dw = to_host(tdw); @show dw
dw2 = zeros(Float32, size(w));  @show dw2
cudnnConvolutionBackwardFilter(x,dy,dw2); @show dw2
@test_approx_eq dw dw2

dx = zeros(Float32, size(x)); tdx = CudaArray(dx); @show dx
cudnnConvolutionBackwardData(tw, tdy, tdx); dx = to_host(tdx); @show dx
dx2 = zeros(Float32, size(x));  @show dx2
cudnnConvolutionBackwardData(w, dy, dx2); @show dx2
@test_approx_eq dx dx2

using CUDNN: PD, CUDNN_POOLING_MAX, cudnnGetPoolingNdForwardOutputDim
# x = rand(Float32,18,18,3,100); tx = CudaArray(x); @show x
x = reshape(Float32[1:25;], 5, 5, 1, 1); tx = CudaArray(x); @show x
psize, padding, stride = 3, 0, 3
pd1 = PD(2, psize, padding, stride, CUDNN_POOLING_MAX)
@assert cudnnGetPoolingNdForwardOutputDim(pd1, tx) == cudnnGetPoolingNdForwardOutputDim(x, window=psize, padding=padding, stride=stride, mode=0)
@show ydims = cudnnGetPoolingNdForwardOutputDim(x, window=psize, padding=padding, stride=stride, mode=0)
y = zeros(Float32, ydims); ty = CudaArray(y);
y2 = zeros(y)
cudnnPoolingForward(tx, ty; window=psize, padding=padding, stride=stride, mode=0); y = to_host(ty); @show y
cudnnPoolingForward(x, y2; window=psize, padding=padding, stride=stride, mode=0); @show y2
@test_approx_eq y y2

dx = zeros(Float32, size(x)); tdx = CudaArray(dx);
dx2 = zeros(Float32, size(x));
dy = rand(Float32, size(y)); tdy = CudaArray(dy); @show dy
cudnnPoolingBackward(ty, tdy, tx, tdx; window=psize, padding=padding, stride=stride, mode=0); dx = to_host(tdx); @show dx
cudnnPoolingBackward(y, dy, x, dx2; window=psize, padding=padding, stride=stride, mode=0); @show dx2
@test_approx_eq dx dx2


text:
using AFLplusplus_jll
using Clang_jll

Clang_jll.clang() do clang
	src = joinpath(@__DIR__, "afl-shim.c")
	out = joinpath(@__DIR__, "afl-shim")

	run(`$clang $src -o $out`)

	AFLplusplus_jll.afl_clang_fast() do afl_clang
		afl_path = joinpath(dirname(afl_clang), "../lib/afl")

		src = joinpath(@__DIR__, "afl-test.c")
		out = joinpath(@__DIR__, "afl-test")

		withenv("AFL_PATH" => afl_path, "AFL_CC" => clang) do
			run(`$afl_clang $src -o $out`)
		end
	end
end