text:
#Basic math
using DSP, Statistics, LinearAlgebra, Random

a = rand(Float32, 10, 10)
ad = AFArray(a)

@testset "Math" begin
    @test sum(abs2, Array(ad + 2) - (a .+ 2)) < 1e-6
    @test sum(abs2, Array(ad .+ 2) - (a .+ 2)) < 1e-6
    @test sum(abs2, Array(2 + ad) - (2 .+ a)) < 1e-6
    @test sum(abs2, Array(2 .+ ad) - (2 .+ a)) < 1e-6
    @test sum(abs2, Array(ad - 2) - (a .- 2)) < 1e-6
    @test sum(abs2, Array(ad .- 2) - (a .- 2)) < 1e-6
    @test sum(abs2, Array(2 - ad) - (2 .- a)) < 1e-6
    @test sum(abs2, Array(2 .- ad) - (2 .- a)) < 1e-6
    @test sum(abs2, Array(ad * 2) - (a * 2)) < 1e-6
    @test sum(abs2, Array(ad .* 2) - (a .* 2)) < 1e-6
    @test sum(abs2, Array(2 * ad) - (2 * a)) < 1e-6
    @test sum(abs2, Array(2 .* ad) - (2 .* a)) < 1e-6
    @test sum(abs2, Array(ad / 2) - (a / 2)) < 1e-6
    @test sum(abs2, Array(ad ./ 2) - (a ./ 2)) < 1e-6
    @test sum(abs2, Array(2 ./ ad) - (2 ./ a)) < 1e-6
    @test sum(abs2, Array(ad .^ 2.0) - (a .^ 2)) < 1e-6
    @test sum(abs2, Array(2 .^ ad) - (2 .^ a)) < 1e-6

    #Trig functions
    @test sum(abs2, Array(sin.(ad)) - sin.(a)) < 1e-6
    @test sum(abs2, Array(cos.(ad)) - cos.(a)) < 1e-6
    @test sum(abs2, Array(tan.(ad)) - tan.(a)) < 1e-6
    @test sum(abs2, Array(sinh.(ad)) - sinh.(a)) < 1e-6
    @test sum(abs2, Array(cosh.(ad)) - cosh.(a)) < 1e-6
    @test sum(abs2, Array(tanh.(ad)) - tanh.(a)) < 1e-6

    @test sum(fill(AFArray, 1, (1, 2))) == 2
    @test sum(fill(AFArray{Float32}, 1, (1, 2))) == 2
    @test sum(fill(AFArray{Float32,2}, 1, (1, 2))) == 2
    @test sum(ones(AFArray{Float32}, (1, 2))) == 2f0
    @test sum(ones(AFArray{Float32,2}, (1, 2))) == 2f0
    @test eltype(zeros(AFArray{Float32}, (1, 2))) == Float32
    @test eltype(zeros(AFArray{Float32,2}, (1, 2))) == Float32
    @test typeof(zeros(ad)) == AFArray{Float32, 2}
    @test typeof(ones(ad)) == AFArray{Float32, 2}
end

@testset "Arrays" begin
    @test err_to_string(Cuint(0)) == "Success"
    @test_throws ErrorException set_device(-5)
    arr1 = @inferred AFArray{Int,1}([1, 2])
    @test eltype(arr1) == Int
    @test ndims(arr1) == 1
    @test (@inferred size(arr1)) == (2,)
    arr2 = @inferred copy(arr1)
    arr3 = @inferred deepcopy(arr2)
    @test typeof(arr3) == typeof(arr2)
    arr4 = @inferred Array{Int,1}(arr1)
    @test arr4 == [1, 2]
    arr5 = AFArray{Int,2}([1 2 3; 4 5 6])
    @test (@inferred size(arr5)) == (2,3)
    @test size(cat(2, arr1, arr5)) == (2, 4)
    @test size(hcat(arr1, arr5)) == (2, 4)
    @test size(vcat(arr1', arr5')) == (4, 2)
    arr6 = @inferred AFArray([1., 2.])
    arr7 = @inferred Array(arr6)
    @test @inferred eltype(arr6) == eltype(arr7)
    @test @inferred ndims(arr6) == ndims(arr7)
    @test @inferred size(arr6) == size(arr7)
    @test @inferred any(@inferred isnan(arr6)) == false
    @test @inferred any(@inferred isinf(arr6)) == false
    @test @inferred any(@inferred any(isnan(arr6), 1)) == false
    @test @inferred any(@inferred any(isinf(arr6), 1)) == false
    @test @inferred all(@inferred isnan(arr6)) == false
    @test @inferred all(@inferred isinf(arr6)) == false
    @test @inferred all(@inferred all(isnan(arr6), 1)) == false
    @test @inferred all(@inferred all(isinf(arr6), 1)) == false
    @test sum(arr5) == 21
    @test sum(arr6) == 3.
    @test sum(complex(arr6)) == 3 + 0im
    @test sum(complex(arr6, arr6)) == 3 + 3im
    @test sum(abs2(arr5)) == 91
    @test sum(abs2(arr6)) == 5.0
    arr9 = AFArray{Complex{Float64},2}([1.0+3im 2. 3.; 4 5 6])
    @test sum(arr9) == 21+3im
    @test sum(abs2(arr9)) == 100.0
    @test sum(conj(arr9)) == 21-3im
    @test sum(real(arr9)) == 21
    @test sum(imag(arr9)) == 3
    @test typeof(@inferred AFArray{Complex{Float32},2}(arr9)) == AFArray{Complex{Float32},2}
    @test typeof(@inferred AFArray{UInt32}(arr9)) == AFArray{UInt32,2}
    b = @inferred(1 + arr1)
    @test sum(b) == 5
end

a1 = rand(2, 3) + rand(2, 3)im
a2 = rand(2, 3) + rand(2, 3)im
af1 = @inferred AFArray(a1)
af2 = @inferred AFArray(a2)

for op in [:+, :-, :.+, :.-, :.*, :./, :.^]
    @test @eval sum($op(a1, a2)) ≈ sum(@inferred $op(af1, af2))
end

a3 = randn(Float32, 2, 3)
af3 = AFArray(a3)
a4 = randn(Float32, 2, 3)
af4 = AFArray(a4)

for op in [:.>, :.>=, :.!=, :.<, :.<=]
    @test @eval all($op(a3, a4) .== Array(@inferred $op(af3, af4)))
    @test @eval all($op(0, a4) .== Array(@inferred $op(0, af4)))
    @test @eval all($op(a3, 0) .== Array(@inferred $op(af3, 0)))
end

for op in [:.==]
    @test @eval all($op(a3, a4) .== Array($op(af3, af4)))
    @test @eval all($op(0, a4) .== Array($op(0, af4)))
    @test @eval all($op(a3, 0) .== Array($op(af3, 0)))
end

c1 = rand()

for op in [:*, :.+, :.-, :.*, :./, :.^]
    @test @eval sum($op(c1, a2)) ≈ sum(@inferred $op(c1, af2))
end

for op in [:*, :/, :.^, :.+, :.-, :.*, :./]
    @test @eval sum($op(a1, c1)) ≈ sum(@inferred $op(af1, c1))
end

as = Int32[-2 -1 0 1 2]
@test all(AFArray(sign.(as)) .== sign(AFArray(as)))
@test !any(!(AFArray(sign.(as)) .== sign(AFArray(as))))

@test (@inferred size(AFArray(rand(1,2,3)))) == (1,2,3)
@test (@inferred size(AFArray(rand(1,2,3,4)))) == (1,2,3,4)

am = rand(Float32, 3, 3)
amf = AFArray(am)

@test sum(amf * amf) ≈ sum(am * am)
@test sum(amf * amf') ≈ sum(am * am')
#@test sum(A_mul_Bt(amf, amf)) ≈ sum(am * am')
@test sum(amf' * amf) ≈ sum(am' * am)
#@test sum(At_mul_B(amf, amf)) ≈ sum(am' * am)
@test sum(amf' * amf') ≈ sum(am' * am')
#@test sum(At_mul_Bt(amf, amf)) ≈ sum(am' * am')

@test Array(sum(amf, 1)) ≈ sum(am, dims=1)
@test Array(sum(amf, 2)) ≈ sum(am, dims=2)
# @test Array(sum(amf, [1,2])) ≈ sum(am, dims=[1,2])
# @test Array(sum(amf, [2,1])) ≈ sum(am, dims=[2,1])
# @test Array(sum(amf, 1:2)) ≈ sum(am, 1:2)
# @test Array(sum(amf, (1,2))) ≈ sum(am, (1,2))

@test Array(prod(amf, 1)) ≈ prod(am, dims=1)
@test Array(prod(amf, 2)) ≈ prod(am, dims=2)

@test all(Array(minimum(amf, 1)) .== minimum(am, dims=1))
@test all(Array(minimum(amf, 2)) .== minimum(am, dims=2))

@test all(Array(maximum(amf, 1)) .== maximum(am, dims=1))
@test all(Array(maximum(amf, 2)) .== maximum(am, dims=2))
@test maximum(amf) == maximum(am)
@test minimum(amf) == minimum(am)
@test mean(amf) ≈ mean(am)
@test std(amf) ≈ std(am)
@test var(amf) ≈ var(am)
@test median(amf) == median(am)

@test Array(vec(amf)) == vec(am)
@test typeof(norm(amf)) == Float32
@test typeof(norm(amf, 1)) == Float32
@test typeof(norm(amf, Inf)) == Float32
@test @inferred(norm(amf)) ≈ norm(am)
@test @inferred(norm(amf, 1)) ≈ norm(am, 1)
@test @inferred(norm(amf, Inf)) ≈ norm(am, Inf)
u,s,v = @inferred(svd(amf))

signal = rand(Float32, 100)
signalf = AFArray(signal)
filt = rand(Float32, 3)
filtf = AFArray(filt)
@test sum(conv(signal, filt)) ≈ sum(@inferred conv(signalf, filtf))

srand(AFArray, rand(Int))

@testset "Sizes" begin
    @test size(rand(AFArray, 1)) == (1,)
    @test size(rand(AFArray, 1, 2)) == (1,2)
    @test size(rand(AFArray{Float64}, 1)) == (1,)
    @test size(rand(AFArray{Float64}, 1, 2)) == (1,2)
    @test size(rand(AFArray{Float64}, 1)) == (1,)
    @test size(rand(AFArray{Float64}, 1, 2)) == (1,2)

    @test size(randn(AFArray, 1)) == (1,)
    @test size(randn(AFArray, 1, 2)) == (1,2)
    @test size(randn(AFArray{Float64}, 1)) == (1,)
    @test size(randn(AFArray{Float64}, 1, 2)) == (1,2)
    @test size(randn(AFArray{Float64}, 1)) == (1,)
    @test size(randn(AFArray{Float64}, 1, 2)) == (1,2)
end

a3 = rand(UInt8, 2, 3)
af3 = AFArray(a3)
a4 = rand(UInt8, 2, 3)
af4 = AFArray(a4)
c1 = rand(UInt8)

for op in [:&, :|, :xor, :min, :max]
    @test @eval all($op.(a3, a4) .== Array(@inferred $op(af3, af4)))
    @test @eval all($op.(c1, a4) .== Array(@inferred $op(c1, af4)))
    @test @eval all($op.(a3, c1) .== Array(@inferred $op(af3, c1)))
end

a3 = rand(UInt64, 2, 3)
af3 = AFArray(a3)
a4 = rand(UInt64, 2, 3)
af4 = AFArray(a4)
c1 = rand(UInt64)

for op in [:&, :|, :xor, :min, :max]
    @test @eval all($op.(a3, a4) .== Array(@inferred $op(af3, af4)))
    @test @eval all($op.(c1, a4) .== Array(@inferred $op(c1, af4)))
    @test @eval all($op.(a3, c1) .== Array(@inferred $op(af3, c1)))
end

a3 = abs.(rand(Int64, 2, 3))
af3 = AFArray(a3)
a4 = mod.(rand(Int64, 2, 3), 10)
af4 = AFArray(a4)
c1 = mod(rand(Int64), 10)

for op in [:<<, :>>, :&, :|, :xor, :min, :max]
    @test @eval all($op.(a3, a4) .== Array($op.(af3, af4)))
    @test @eval all($op.(a3, a4) .== Array(@inferred $op(af3, af4)))
    @test @eval all($op.(c1, a4) .== Array(@inferred $op(c1, af4)))
    @test @eval all($op.(a3, c1) .== Array(@inferred $op(af3, c1)))
end

a3 = rand(Bool, 2, 3)
af3 = AFArray(a3)
a4 = rand(Bool, 2, 3)
af4 = AFArray(a4)
c1 = rand(Bool)

for op in [:&, :|, :xor, :min, :max]
    @test @eval all($op.(a3, a4) .== Array(@inferred $op(af3, af4)))
    @test @eval all($op.(c1, a4) .== Array(@inferred $op(c1, af4)))
    @test @eval all($op.(a3, c1) .== Array(@inferred $op(af3, c1)))
end

a5 = rand(Float32, 5)
af5 = AFArray(a5)

for op in [:abs, :acos, :asin, :atan, :cbrt, :ceil, :cos, :cosh,
           :exp, :expm1, :floor, :erfc,
           :imag, :isinf, :isnan, :lgamma, :log, :log10,
           :log1p, :log2, :real, :round, :sign, :signbit, :sin, :sinh, :sqrt, :tan,
           :tanh, :trunc]
    @testset "$op" begin
        @test @eval $op.(a5) ≈ Array(@inferred $op(af5))
    end
end

a6 = rand(Float32, 3, 3)
af6 = AFArray(a6)

for op in [:identity, :vec, :transpose]
    @testset "$op" begin
        @test @eval $op(a6) ≈ Array(@inferred $op(af6))
    end
end

for op in [:maximum, :minimum, :sum, :prod, :sizeof, :var, :std, :mean, :median]
    @testset "$op" begin
        @test @eval $op(a6) ≈ @inferred $op(af6)
    end
end

aa = [true false true]
bb = [1 2 3]
cc = [4 5 6; 7 8 9]
a = AFArray(aa)
b = AFArray(bb)
c = AFArray(cc)

@test Array(select(a, b, c)) == ifelse.(aa, bb, cc)
@test Array(select(a, 0, b)) == ifelse.(aa, 0, bb)
@test Array(select(a, b, 0)) == ifelse.(aa, bb, 0)

@test size(reshape(c, 6)) == (6, )
@test size(reshape(c, (3,2))) == (3, 2)

#plot(rand(AFArray, 10), rand(AFArray, 10))

s = rand(AFArray{Float32}, 10)
sh = Array(s)
@testset "Sort" begin
    @test Array(@inferred sort(s)) == sort(sh)
    @test Array(@inferred sortperm(s)) == sortperm(sh)
end

@testset "IR" begin
    x = AFArray([1f0; 2f0; 3f0])
    b = AFArray([1.; 2.])
    @test Array(@inferred fir(b, x)) == [1f0; 4f0; 7f0]
    a = AFArray([1f0])
    b = AFArray{Float32}(b)
    @test Array(@inferred iir(b, a, x)) == [1f0; 4f0; 7f0]
end

ss = similar(s)
@test typeof(ss) == typeof(s)
ss[1] = 10
@test s[1] != 10f0 && ss[1] == 10f0
ssa = similar(s, Float64)
@test typeof(ssa) == AFArray{Float64,1}

aa = [NaN Inf 3]
a = AFArray(aa)

@test Array(@inferred isfinite(a)) == isfinite.(aa)


text:
# GMM estimation for a sample from Chi^2(theta)
# compare to two method of moments estimators (see chi2mm.m)
using Distributions, Plots
function main()
n = 30
theta = 3.0
reps = 1000
results = zeros(reps)
W = eye(2)
for i = 1:1000
    y = rand(Chisq(theta), n)
    obj = theta -> ([theta.-mean(y); theta.-0.5*var(y)]'W*[theta.-mean(y); theta.-0.5*var(y)])
    thetahat, junk, junk = fminunc(obj, [3.0])
    results[i] = sqrt(n)*(thetahat[1]-theta)
end
histogram(results,nbins=50)
gui()
return
end
main()

ext: jl

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
using Test
include("sudoku_local.jl")
# add here the file with your unit tests


text:
"""
    @pywith context symbol block

This is a python context-manager block.  The above is equivalent to
```python
with context as symbol:
    block
```
in Python.

___TODO:___ Currently this doesn't handle exceptions the way it's supposed to due to Python
and Julia having different scoping rules for their `try`, `catch` blocks.  Making this work
correctly would require a somewhat complicated macro that looked for all assignments in the
block and repeated them after the `try` construct.
"""
macro pywith(context, as::Symbol, block::Expr)
    # is it possible to escape only portions of the output, or do I have to do this?
    sys = gensym()
    mgr = gensym()
    exit = gensym()
    value = gensym()
    esc(quote
        @pyimport sys as $sys
        $mgr = $context
        $exit = pytypeof($mgr)[:__exit__]
        $value = pytypeof($mgr)[:__enter__]($mgr)
        let $as = $value
            $block
        end
        $exit($mgr, $sys.exc_info()...)
    end)
end

macro pywith(context, block::Expr)
    sys = gensym()
    mgr = gensym()
    exit = gensym()
    esc(quote
        @pyimport sys
        $mgr = $context
        $exit = pytypeof($mgr)[:__exit__]
        pytypeof($mgr)[:__enter__]($mgr)
        $block
        $exit($mgr, sys.exc_info()...)
    end)
end
export @pywith


text:
@testset "Vector spaces" begin
    @testset "VectorSpace" begin
        𝒫 = ParameterSpace()
        𝒯 = Taylor(0)
        # different types of vector spaces cannot be equal
        @test !(𝒫 == 𝒯)
        # different types of vector spaces cannot be subset of one another
        @test !(𝒫 ⊆ 𝒯)
        # intersection and union between different types of vector spaces is not allowed
        @test_throws MethodError intersect(𝒫, 𝒯)
        @test_throws MethodError union(𝒫, 𝒯)
    end

    @testset "ParameterSpace" begin
        𝒫 = ParameterSpace()
        @test 𝒫 ⊆ 𝒫
        @test 𝒫 ∩ 𝒫 == 𝒫 ∪ 𝒫 == 𝒫
        @test dimension(𝒫) == 1
        @test RadiiPolynomial._firstindex(𝒫) == 1
        @test RadiiPolynomial._lastindex(𝒫) == 1
        @test indices(𝒫) == Base.OneTo(1)
        @test RadiiPolynomial._findposition(1, 𝒫) == 1
    end

    @testset "Taylor" begin
        𝒯₂, 𝒯₃ = Taylor(2), Taylor(3)
        @test !(𝒯₂ == 𝒯₃)
        @test 𝒯₂ ⊆ 𝒯₃
        @test 𝒯₂ ∩ 𝒯₃ == 𝒯₂
        @test 𝒯₂ ∪ 𝒯₃ == 𝒯₃
        @test order(𝒯₂) == 2
        @test dimension(𝒯₂) == 3
        @test RadiiPolynomial._firstindex(𝒯₂) == 0
        @test RadiiPolynomial._lastindex(𝒯₂) == 2
        @test indices(𝒯₂) == 0:2
        @test RadiiPolynomial._findindex_constant(𝒯₂) == 0
        @test RadiiPolynomial._findposition(1, 𝒯₂) == 2
        @test RadiiPolynomial._findposition(0:2, 𝒯₂) == 1:3
        @test RadiiPolynomial._findposition([0, 2], 𝒯₂) == [1, 3]
        @test RadiiPolynomial._findposition(:, 𝒯₂) == Colon()
    end

    @testset "Fourier" begin
        ℱ₂, ℱ₃ = Fourier(2, 1.0), Fourier(3, 1.0)
        @test !(ℱ₂ == ℱ₃)
        @test ℱ₂ ⊆ ℱ₃
        @test ℱ₂ ∩ ℱ₃ == ℱ₂
        @test ℱ₂ ∪ ℱ₃ == ℱ₃
        @test order(ℱ₂) == 2
        @test frequency(ℱ₂) == 1.0
        @test dimension(ℱ₂) == 5
        @test RadiiPolynomial._firstindex(ℱ₂) == -2
        @test RadiiPolynomial._lastindex(ℱ₂) == 2
        @test indices(ℱ₂) == -2:2
        @test RadiiPolynomial._findindex_constant(ℱ₂) == 0
        @test RadiiPolynomial._findposition(-1, ℱ₂) == 2
        @test RadiiPolynomial._findposition(0:1, ℱ₂) == 3:4
        @test RadiiPolynomial._findposition([-2, 2], ℱ₂) == [1, 5]
        @test RadiiPolynomial._findposition(:, ℱ₂) == Colon()
        @test convert(Fourier{Float64}, ℱ₂) == ℱ₂
        @test convert(Fourier{Int}, ℱ₂) == Fourier(2, 1)
        @test promote_type(Fourier{Float64}, Fourier{Float64}) == Fourier{Float64}
        @test promote_type(Fourier{Float64}, Fourier{Int}) == Fourier{Float64}
    end

    @testset "Chebyshev" begin
        𝒞₂, 𝒞₃ = Chebyshev(2), Chebyshev(3)
        @test !(𝒞₂ == 𝒞₃)
        @test 𝒞₂ ⊆ 𝒞₃
        @test 𝒞₂ ∩ 𝒞₃ == 𝒞₂
        @test 𝒞₂ ∪ 𝒞₃ == 𝒞₃
        @test order(𝒞₂) == 2
        @test dimension(𝒞₂) == 3
        @test RadiiPolynomial._firstindex(𝒞₂) == 0
        @test RadiiPolynomial._lastindex(𝒞₂) == 2
        @test indices(𝒞₂) == 0:2
        @test RadiiPolynomial._findindex_constant(𝒞₂) == 0
        @test RadiiPolynomial._findposition(1, 𝒞₂) == 2
        @test RadiiPolynomial._findposition(0:2, 𝒞₂) == 1:3
        @test RadiiPolynomial._findposition([0, 2], 𝒞₂) == [1, 3]
        @test RadiiPolynomial._findposition(:, 𝒞₂) == Colon()
    end

    @testset "TensorSpace" begin
        𝒯₂, ℱ₂, 𝒞₂ = Taylor(2), Fourier(2, 1.0), Chebyshev(2)
        𝑇 = TensorSpace((𝒯₂, ℱ₂, 𝒞₂))
        @test spaces(𝑇) == (𝒯₂, ℱ₂, 𝒞₂)
        @test 𝒯₂ ⊗ ℱ₂ ⊗ 𝒞₂ == (𝒯₂ ⊗ ℱ₂) ⊗ 𝒞₂ == 𝒯₂ ⊗ (ℱ₂ ⊗ 𝒞₂) == 𝑇
        @test 𝑇 ⊗ 𝑇 == TensorSpace((𝒯₂, ℱ₂, 𝒞₂, 𝒯₂, ℱ₂, 𝒞₂))
        @test 𝑇[1] == 𝒯₂
        @test 𝑇[1:3] == 𝑇[[1, 2, 3]] == 𝑇[:] == 𝑇
        @test Base.tail(𝑇) == TensorSpace((ℱ₂, 𝒞₂))
        @test Base.front(𝑇) == TensorSpace((𝒯₂, ℱ₂))
        @test 𝑇 ⊆ Taylor(3) ⊗ ℱ₂ ⊗ 𝒞₂
        @test 𝑇 ∩ (Taylor(3) ⊗ Fourier(1, 1.0) ⊗ Chebyshev(3)) == 𝒯₂ ⊗ Fourier(1, 1.0) ⊗ 𝒞₂
        @test 𝑇 ∪ (Taylor(3) ⊗ Fourier(1, 1.0) ⊗ Chebyshev(3)) == Taylor(3) ⊗ ℱ₂ ⊗ Chebyshev(3)
        @test dimension(𝑇) == 45
        @test dimension(𝑇, 1) == 3
        @test dimensions(𝑇) == (3, 5, 3)
        @test RadiiPolynomial._firstindex(𝑇) == (0, -2, 0)
        @test RadiiPolynomial._lastindex(𝑇) == (2, 2, 2)
        @test indices(𝑇) == Base.Iterators.product(0:2, -2:2, 0:2)
        @test RadiiPolynomial._findindex_constant(𝑇) == (0, 0, 0)
        @test RadiiPolynomial._findposition((1, 2, 0), 𝑇) == 14
        @test RadiiPolynomial._findposition((:, -2:2, :), 𝑇) == collect(1:45)
        @test RadiiPolynomial._findposition([(1, 2, 0), (2, 2, 2)], 𝑇) == [14, 45]
        @test RadiiPolynomial._findposition((:, :, :), 𝑇) == Colon()
        @test RadiiPolynomial._findposition(:, 𝑇) == Colon()
        @test order(𝑇) == (2, 2, 2)
        @test order(𝑇, 1) == 2
        @test_throws MethodError frequency(𝑇)
        @test frequency(𝑇, 2) == 1.0
        @test convert(TensorSpace{Tuple{Taylor,Fourier{Float64},Chebyshev}}, 𝑇) == 𝑇
        @test convert(TensorSpace{Tuple{Taylor,Fourier{Int},Chebyshev}}, 𝑇) == Taylor(2) ⊗ Fourier(2, 1) ⊗ Chebyshev(2)
    end

    @testset "CartesianSpace" begin
        𝒯₁¹ = CartesianPower(Taylor(1), 1)
        𝒞₂² = CartesianPower(Chebyshev(2), 2)
        ℱ₃³ = CartesianPower(Fourier(3, 1.0), 3)
        𝒯₁¹𝒞₂²ℱ₃³ = CartesianProduct((𝒯₁¹, 𝒞₂², ℱ₃³))
        𝒯₂¹𝒞₃²ℱ₃³ = CartesianProduct((CartesianPower(Taylor(2), 1), CartesianPower(Chebyshev(3), 2), ℱ₃³))
        @test space(𝒯₁¹) == Taylor(1)
        @test spaces(𝒯₁¹) == [Taylor(1)]
        @test spaces(𝒯₁¹𝒞₂²ℱ₃³) == (𝒯₁¹, 𝒞₂², ℱ₃³)
        @test nb_cartesian_product(ℱ₃³) == nb_cartesian_product(𝒯₁¹𝒞₂²ℱ₃³) == 3
        @test 𝒯₁¹^1 == CartesianPower(𝒯₁¹, 1)
        @test 𝒯₁¹ × 𝒞₂² × ℱ₃³ == 𝒯₁¹ × (𝒞₂² × ℱ₃³) == (𝒯₁¹ × 𝒞₂²) × ℱ₃³ == (𝒯₁¹ × 𝒞₂²) × CartesianProduct(tuple(ℱ₃³)) == 𝒯₁¹𝒞₂²ℱ₃³
        @test Base.front(𝒯₁¹𝒞₂²ℱ₃³) == 𝒯₁¹ × 𝒞₂²
        @test Base.tail(𝒯₁¹𝒞₂²ℱ₃³) == 𝒞₂² × ℱ₃³
        @test 𝒞₂²[1] == Chebyshev(2)
        @test 𝒞₂²[1:2] == 𝒞₂²[[1, 2]] == 𝒞₂²[:] == 𝒞₂²
        @test 𝒯₁¹𝒞₂²ℱ₃³[1] == 𝒯₁¹
        @test 𝒯₁¹𝒞₂²ℱ₃³[1:3] == 𝒯₁¹𝒞₂²ℱ₃³[[1, 2, 3]] == 𝒯₁¹𝒞₂²ℱ₃³[:] == 𝒯₁¹𝒞₂²ℱ₃³
        @test 𝒯₁¹𝒞₂²ℱ₃³ ⊆ 𝒯₂¹𝒞₃²ℱ₃³
        @test 𝒯₁¹𝒞₂²ℱ₃³ ∩ 𝒯₂¹𝒞₃²ℱ₃³ == 𝒯₁¹𝒞₂²ℱ₃³
        @test 𝒯₁¹𝒞₂²ℱ₃³ ∪ 𝒯₂¹𝒞₃²ℱ₃³ == 𝒯₂¹𝒞₃²ℱ₃³
        @test dimension(𝒯₁¹𝒞₂²ℱ₃³) == 29
        @test dimension(𝒯₁¹𝒞₂²ℱ₃³, 1) == dimension(𝒯₁¹, 1) == 2
        @test dimensions(𝒯₁¹) == [2]
        @test dimensions(𝒯₁¹𝒞₂²ℱ₃³) == (2, 6, 21)
        @test RadiiPolynomial._firstindex(𝒯₁¹𝒞₂²ℱ₃³) == RadiiPolynomial._firstindex(𝒞₂²) == 1
        @test RadiiPolynomial._lastindex(𝒞₂²) == 6
        @test RadiiPolynomial._lastindex(𝒯₁¹𝒞₂²ℱ₃³) == 29
        @test indices(𝒞₂²) == Base.OneTo(6)
        @test indices(𝒯₁¹𝒞₂²ℱ₃³) == Base.OneTo(29)
        @test order(𝒯₁¹𝒞₂²ℱ₃³) == ([1], [2, 2], [3, 3, 3])
        @test order(𝒯₁¹𝒞₂²ℱ₃³, 1) == [1]
        @test order(𝒯₁¹, 1) == 1
        @test_throws MethodError frequency(𝒯₁¹𝒞₂²ℱ₃³)
        @test frequency(𝒯₁¹𝒞₂²ℱ₃³, 3) == frequency(ℱ₃³) == [1.0, 1.0, 1.0]
        @test frequency(ℱ₃³, 1) == 1.0
        @test RadiiPolynomial._component_findposition(1, 𝒯₁¹) == RadiiPolynomial._component_findposition(1:1, 𝒯₁¹) == 1:2
        @test RadiiPolynomial._component_findposition(2, 𝒞₂²) == RadiiPolynomial._component_findposition(2:2, 𝒞₂²) == 4:6
        @test RadiiPolynomial._component_findposition(1, 𝒯₁¹𝒞₂²ℱ₃³) == RadiiPolynomial._component_findposition(1:1, 𝒯₁¹𝒞₂²ℱ₃³) == 1:2
        @test RadiiPolynomial._component_findposition(2, 𝒯₁¹𝒞₂²ℱ₃³) == RadiiPolynomial._component_findposition(2:2, 𝒯₁¹𝒞₂²ℱ₃³) == 3:8
        @test convert(CartesianProduct{Tuple{CartesianPower{Taylor}, CartesianPower{Chebyshev}, CartesianPower{Fourier{Float64}}}}, 𝒯₁¹𝒞₂²ℱ₃³) == 𝒯₁¹𝒞₂²ℱ₃³
        @test convert(CartesianProduct{Tuple{CartesianPower{Taylor}, CartesianPower{Chebyshev}, CartesianPower{Fourier{Int}}}}, 𝒯₁¹𝒞₂²ℱ₃³) == 𝒯₁¹ × 𝒞₂² × Fourier(3, 1)^3
    end
end


text:
# Part 1

mutable struct Position
    horizontal::Int64
    depth::Int64
end
Position() = Position(0, 0)

function move!(position::Position, direction::Symbol, amount::Integer)
    if direction == :forward
        position.horizontal += amount
    elseif direction == :down
        position.depth += amount
    elseif direction == :up
        position.depth -= amount
    end
end

function parse_line(line::AbstractString)
    direction, amount = split(line, " ")
    return (Symbol(direction), parse(Int64, amount))
end

function get_final_position(commands::AbstractString, position=Position())
    for line in split(commands, "\n")
        isempty(line) && continue
        direction, amount = parse_line(line)
        move!(position, direction, amount)
    end
    return position
end

"""Multiply two components to get final answer."""
get_ans(position::Position) = position.horizontal * position.depth
function get_ans(commands::AbstractString, initial=Position())
    return get_ans(get_final_position(commands, initial))
end

example = """
forward 5
down 5
forward 8
up 3
down 8
forward 2
"""

@assert get_ans(example) == 150

function _load_input(name::AbstractString)
    path = joinpath("inputs", name)
    return read(path, String)
end

commands = _load_input("day_02.txt")
get_ans(commands)

# Part 2

mutable struct PositionAim
    horizontal::Int64
    depth::Int64
    aim::Int64
end
PositionAim() = PositionAim(0, 0, 0)

get_ans(position::PositionAim) = position.horizontal * position.depth

function move!(state::PositionAim, direction::Symbol, amount::Integer)
    if direction == :forward
        state.horizontal += amount
        state.depth += state.aim * amount
    elseif direction == :down
        state.aim += amount
    elseif direction == :up
        state.aim -= amount
    end
end

@assert get_ans(example, PositionAim()) == 900

get_ans(commands, PositionAim())


text:
export
    CBSEnv

# CBS submodule
module CBSEnv

using ..CRCBS
using Parameters, LightGraphs, DataStructures
using GraphUtils: get_graph, get_vtx

@with_kw struct State <: AbstractGraphState
    vtx::Int        = -1 # vertex of graph
    t::Int          = -1
end
Base.convert(::Type{State},s::GraphState) = State(vtx=get_vtx(s),t=get_t(s))
# CRCBS.get_vtx(s::State) = s.vtx
CRCBS.get_t(s::State) = s.t
@with_kw struct Action <: AbstractGraphAction
    e::Edge{Int}    = Edge(-1,-1)
    dt::Int         = 1
end
Base.convert(::Type{Action},a::GraphAction) = Action(e=get_e(a),dt=get_dt(a))
CRCBS.get_e(a::Action) = a.e
CRCBS.get_dt(a::Action) = a.dt
@with_kw struct LowLevelEnv{C<:AbstractCostModel,H<:AbstractCostModel,G<:AbstractGraph,T} <: GraphEnv{State,Action,C}
    graph::G                    = Graph()
    goal::State                 = State()
    agent_idx::Int              = -1
    constraints::T              = discrete_constraint_table(nv(graph),nv(graph)^2,agent_idx)
    cost_model::C               = SumOfTravelTime()
    heuristic::H                = NullHeuristic() #PerfectHeuristic(graph,Vector{Int}(),Vector{Int}())
end
# CRCBS.get_graph(env::LowLevelEnv)            = env.graph
# GraphUtils.get_graph(env::LowLevelEnv)            = env.graph
# CRCBS.get_agent_id(env::LowLevelEnv)         = env.agent_idx
# CRCBS.get_constraints(env::LowLevelEnv)      = env.constraints
# CRCBS.get_goal(env::LowLevelEnv)             = env.goal
# CRCBS.get_cost_model(env::LowLevelEnv)       = env.cost_model
# CRCBS.get_heuristic_model(env::LowLevelEnv)  = env.heuristic
CRCBS.base_env_type(env::LowLevelEnv)        = LowLevelEnv

CRCBS.get_possible_actions(env::LowLevelEnv,s)  = map(v->Action(e=Edge(get_vtx(s),v)),outneighbors(get_graph(env),get_vtx(s)))
CRCBS.get_next_state(s::State,a::Action)    = State(get_e(a).dst,get_t(s)+get_dt(a))
CRCBS.get_next_state(env::LowLevelEnv,s,a)  = get_next_state(s,a)
CRCBS.wait(s::State)                        = Action(e=Edge(get_vtx(s),get_vtx(s)))
CRCBS.wait(env::LowLevelEnv,s)              = Action(e=Edge(get_vtx(s),get_vtx(s)))

function CRCBS.build_env(mapf::MAPF{E,S,G}, node::ConstraintTreeNode, idx::Int) where {S,G,E<:LowLevelEnv}
    if idx > 0
        t_goal = -1
        n = PathNode{State,Action}(sp=mapf.goals[idx])
        s_constraints, _ = search_constraints(mapf.env,get_constraints(node,idx),n)
        for c in s_constraints
            t_goal = max(t_goal,get_time_of(c)+1)
        end
        return typeof(mapf.env)(
            graph = get_graph(mapf.env),
            constraints = get_constraints(node,idx),
            goal = State(mapf.goals[idx],t=t_goal),
            agent_idx = idx,
            cost_model = get_cost_model(mapf.env),
            heuristic = get_heuristic_model(mapf.env),
            )
    else
        return typeof(mapf.env)(
            graph = get_graph(mapf.env),
            goal = state_type(mapf)(),
            agent_idx = idx,
            cost_model = get_cost_model(mapf.env),
            heuristic = get_heuristic_model(mapf.env),
            )
    end
end

end


text:
using Test
using EmonCMS
using JuliaDB
using Dates
using Unitful

function maketestfeed(N, interval, t0)
  peak1 = 200.0
  peak2 = 300.0
  result = []
  Ni, Nj, Nk, Nl, Nm = sort(rand(1:N,5))
  t0 -= interval
  t0 *= 1000
  for (t,P) in zip(1:Ni, range(0,peak1;length=Ni))
    push!(result, Union{Int64,Float64}[t0+t*interval*1000,P])
  end
  for t in Nj:Nk
    push!(result, Union{Int64,Float64}[t0+t*interval*1000,peak2])
  end
  for t in Nl:Nm
    push!(result, Union{Int64,Float64}[t0+t*interval*1000,peak2])
  end

  Nmissing = Nj-Ni-1 + Nl-Nk-1

  return result, (Ni-1)*interval*peak1/2 + ((Nk-Nj) + (Nm-Nl))*interval*peak2, Nmissing
end

function getblock(::Any, feed, starttime, endtime, interval)
  starttime *= 1000
  endtime *= 1000
  tarray = getindex.(feed,1)
  startidx = searchsortedfirst(tarray,starttime)
  endidx = searchsortedlast(tarray,endtime)
  if endidx > length(tarray)
    return feed[startidx:end]
  end
  return feed[startidx:endidx]
end

# Set up a dummy database first
tdir = mktempdir()
EmonCMS.writeconnfile(EmonCMS.Connection("https://localhost", "dummykey"), EmonCMS.connectionfile(tdir))
emonds = EmonDataSet(tdir)

@testset "Feed setup" begin
  feedstuple = (id=[2,3], unit=["W", "W"], name=["TestFeed", "TotalPower"], starttime=[5000,5000], interval=[10,10])
  feedstable = table(feedstuple, pkey=[:id, :name])
  save(feedstable, EmonCMS.feedsfile(emonds.path))
end

@testset "Feed reading" begin
  N = 50000
  feeds = feedlist(emonds)
  interval = feeds[1][:interval]
  starttime = feeds[1][:starttime]
  feed, expectedenergy, expectednmissing = maketestfeed(N, interval, starttime)
  @test (feed[1][1] ÷ 1000) == starttime
  endtime = feed[end][1] ÷ 1000

  feedfile = joinpath(tdir, "TestFeed.juliadb")
  save(EmonCMS.updatefeedtable(getblock, nothing, EmonCMS.createfeedtable(), feed, starttime, endtime, interval, "TestFeed"), feedfile)
  feedtable = load(feedfile)
  save(feedtable, joinpath(tdir, "TotalPower.juliadb"))
  push!(feed, [(endtime+interval)*1000,10])
  feedtable2 = EmonCMS.updatefeedtable(getblock, nothing, feedtable, feed, starttime, endtime+interval, interval, "TestFeed")
  @test all(select(feedtable,:time) .== select(feedtable2,:time)[1:end-1])
  @test length(feedtable2) == length(feedtable) + 1
  testfeed = getfeed(emonds, "TestFeed")
  @test length(testfeed) == length(feed) + expectednmissing-1
  @test count(ismissing, select(testfeed,:value)) == expectednmissing
  energytable = energyperperiod(testfeed, Second(interval))
  energy = select(energytable, :energy)
  @test sum(skipmissing(energy)) ≈ (expectedenergy * u"J" |> u"kW*hr")
end

@testset "Postprocessing" begin
  names, energies, counts = energysummary(emonds; years=[1970], period=Day(1), totalpowerfeeds=["TotalPower", "TotalPower"])
  @test length(names) == 2
  @test names[2] == "TestFeed"
  @test names[1] == "Unknown"
  @test all(skipmissing(energies[:,1]) .== skipmissing(energies[:,2]))
end

@testset "Export" begin
  csvdir = mktempdir()
  exportcsv(emonds, csvdir)
  for tname in select(feedlist(emonds), :name)
    ref = EmonCMS.loadfeed(emonds,tname)
    imported = EmonCMS.importfeed(joinpath(csvdir, "$tname.csv"))
    @test all(select(ref,:time) .== select(imported,:time))
    @test all(skipmissing(select(ref,:value) .≈ select(imported,:value)))
    @test count(ismissing, select(ref,:value)) == count(ismissing, select(imported,:value))
  end
  @test feedlist(emonds) == loadtable(joinpath(csvdir,"feedlist.csv"))
end


text:
#### Landmarks specification
import Bridge: _b, _b!, B!, σ!, b!, σ, b, auxiliary

"""
    MarslandShardlow{T} <: ContinuousTimeProcess{State{PointF}}

## Arguments
- `a`: Hamiltonian kernel parameter
- `c`: kernel multiplicate parameter
- `γ`:  noise level
- `λ`:  mean reversion parameter (heath-bath parameter in Marsland-Shardlow (2017))
-  `n`: number of landmarks
"""
struct MarslandShardlow{T} <: ContinuousTimeProcess{State{PointF}}
    a::T # kernel std parameter
    c::T # kernel multiplicate parameter
    γ::T # noise level
    λ::T # mean reversion
    n::Int
end

"""
    MarslandShardlowAux{S,T} <: ContinuousTimeProcess{State{PointF}}

## Arguments
- `a`: Hamiltonian kernel parameter
- `c`: kernel multiplicate parameter
- `γ`:  noise level
- `λ`:  mean reversion parameter (heath-bath parameter in Marsland-Shardlow (2017))
- `xT`: State at time T used for constructing the auxiliary process
-  `n`: number of landmarks
-  `gramT`: Gram-matrix at xT
"""
struct MarslandShardlowAux{S,T} <: ContinuousTimeProcess{State{PointF}}
    a::T        # kernel std parameter
    c::T        # kernel multiplicate parameter
    γ::T        # noise level
    λ::T        # mean reversion
    xT::State{Point{S}}  # use x = State(P.v, zero(P.v)) where v is conditioning vector
    n::Int
    gramT::Matrix{T} # Gram-matrix at xT
end

struct Noisefield{T}
    δ::Point{T}         # locations of noise field
    γ::Point{T}         # scaling at noise field (used to be called lambda)
    τ::T                # std of Gaussian kernel noise field
end

"""
    Landmarks{S,T} <: ContinuousTimeProcess{State{PointF}}

## Arguments
- `a`: Hamiltonian kernel parameter
- `c`: kernel multiplicate parameter
-  `n`:: Int64 number of landmarks
- `db`::Float64 square domain bound used for construction of noise fields
- `nfstd`:  standard deviation of noisefields (assumed to be the same for all noisefields)
- `nfs`:  vector of noisefields
"""
struct  Landmarks{S,T} <: ContinuousTimeProcess{State{PointF}}
    a::T            # kernel std
    c::T            # kernel multiplicate parameter
    n::Int64        # numer of landmarks
    db::Array{Float64,1} # domainbound
    nfstd::Float64
    nfs::Vector{Noisefield{S}}  # vector containing pars of noisefields
end

"""
    LandmarksAux{S,T} <: ContinuousTimeProcess{State{PointF}}

## Arguments
- `a`: Hamiltonian kernel parameter
- `c`: kernel multiplicate parameter
- `xT`: State at time T used for constructing the auxiliary process
-  `n`: number of landmarks
- `nfs`:  vector of noisefields
"""
struct LandmarksAux{S,T} <: ContinuousTimeProcess{State{PointF}}
    a::T            # kernel std
    c::T            # kernel multiplicate parameter
    xT::State{Point{S}}     # use x = State(P.v, zero(P.v)) where v is conditioning vector
    n::Int64        # numer of landmarks
    nfs::Vector{Noisefield{S}}  # vector containing pars of noisefields
    gramT::Matrix{T} # Gram-matrix at xT
end

"""
    MarslandShardlowAux(P::MarslandShardlow, xT, gramT) = MarslandShardlowAux(P.a,P.c, P.γ, P.λ, xT, P.n, gramT)
"""
MarslandShardlowAux(P::MarslandShardlow, xT, gramT) = MarslandShardlowAux(P.a,P.c, P.γ, P.λ, xT, P.n, gramT)

"""
    LandmarksAux(P::Landmarks, xT, gramT) = LandmarksAux(P.a,P.c, xT, P.n, P.nfs, gramT)
"""
LandmarksAux(P::Landmarks, xT, gramT) = LandmarksAux(P.a,P.c, xT, P.n, P.nfs, gramT)

"""
    auxiliary(P::Union{MarslandShardlow, Landmarks}, xT, gramT)

Construct auxiliary process corresponding to `P` and `xT`
"""
function auxiliary(P::Union{MarslandShardlow, Landmarks}, xT, gramT)
    if isa(P,MarslandShardlow)
        return MarslandShardlowAux(P, xT, gramT)
    elseif isa(P,Landmarks)
        return LandmarksAux(P, xT, gramT)
    end
end

const LandmarkModel = Union{Landmarks, LandmarksAux, MarslandShardlow, MarslandShardlowAux}

Bridge.constdiff(::Union{MarslandShardlow, MarslandShardlowAux,LandmarksAux}) = true
Bridge.constdiff(::Landmarks) = false

"""
    kernel(q, P::LandmarkModel)

kernel in Hamiltonian: P.c * exp(-Bridge.inner(q)/(2*P.a^2))
"""
function kernel(q, P::LandmarkModel)
   P.c * exp(-Bridge.inner(q)/(2*P.a^2))
end

"""
    ∇kernel(q, P::LandmarkModel)

gradient of kernel in hamiltonian
"""
function ∇kernel(q, P::LandmarkModel)
    -P.c * P.a^(-2) * kernel(q, P) * q
end

"""
Needed for b! in case P is auxiliary process
"""
function ∇kernel(q, qT, P::LandmarkModel)
     -P.a^(-2) *P.c  * kernel(qT, P) * q
end

"""
    hamiltonian(x, P)

Hamiltonian for deterministic part of landmarks model
"""
function hamiltonian(x, P)
    s = 0.0
    for i in axes(x, 2), j in axes(x, 2)
        s += 1/2*dot(x.p[i], x.p[j])*kernel(x.q[i] - x.q[j], P)
    end
    s
end

Bridge.b(t::Float64, x, P::Union{Landmarks,MarslandShardlow})= Bridge.b!(t, x, copy(x), P)

Bridge.σ(t, x, dm, P) =  Bridge.σ!(t, x, dm , 0*x, P)

################ MS model #########################################################################################

"""
    Bridge.b!(t, x, out, P::MarslandShardlow)

Evaluate drift of landmarks in (t,x) and save to out
x is a state and out as well
"""
function Bridge.b!(t, x, out, P::MarslandShardlow)
    zero!(out)
    pp = out.p
    qq = out.q
    @inbounds  for i in 1:P.n
        for j in 1:P.n
            qq[i] += p(x,j) * kernel(q(x,i) - q(x,j), P)
            pp[i] += -P.λ * qq[i] - dot(p(x,i), p(x,j)) * ∇kernel(q(x,i) - q(x,j), P)
        end
    end
    out
end

"""
    Bridge.b!(t, x, out, Paux::MarslandShardlowAux)

Evaluate drift of landmarks auxiliary process in (t,x) and save to out
x is a state and out as well
"""
function Bridge.b!(t, x, out, Paux::MarslandShardlowAux, xT)
    zero!(out)
    pp = out.p
    qq = out.q
    @inbounds  for i in 1:Paux.n
        for j in 1:Paux.n
            qq[i] += p(x,j) * Paux.gramT[i,j]
        end
    end
    pp .= - Paux.λ * qq
    out
end

"""
    Bridge.B(t, Paux::MarslandShardlowAux)

Compute tildeB(t) for landmarks auxiliary process
"""
function Bridge.B(t, Paux::MarslandShardlowAux, xT) # not AD safe
    X = zeros(UncF, 2Paux.n, 2Paux.n)
    @inbounds for i in 1:Paux.n
        for j in 1:Paux.n
            X[2i-1,2j] =  Paux.gramT[i,j] * one(UncF)
            X[2i,2j] = -Paux.λ * Paux.gramT[i,j] * one(UncF)
        end
    end
    X
end

"""
    Bridge.B!(t,X,out, Paux::MarslandShardlowAux)

Compute B̃(t) * X (B̃ from auxiliary process) and write to out
Both B̃(t) and X are of type UncMat
"""
function Bridge.B!(t,X,out, Paux::MarslandShardlowAux, xT)
    out .= 0.0 * out
    @inbounds  for i in 1:Paux.n  # separately loop over even and odd indices
        for k in 1:2Paux.n # loop over all columns
            for j in 1:Paux.n
                 out[2i-1,k] += Paux.gramT[i,j] * X[p(j), k]
                 out[2i,k] += -Paux.λ * Paux.gramT[i,j] * X[p(j), k]
            end
        end
    end
    out
end

function Bridge.β(t, Paux::MarslandShardlowAux, xT) # Not AD save
    State(zeros(PointF,Paux.n), zeros(PointF,Paux.n))
end

"""
    Bridge.σ!(t, x, dm, out, P::MarslandShardlow)

Compute σ(t,x) * dm and write to out
"""
function Bridge.σ!(t, x, dm, out, P::MarslandShardlow)
    zero!(out)
    out.p .= dm*P.γ
    out
end

"""
    Bridge.σ!(t, x, dm, out, P::MarslandShardlowAux, xT)

Compute σ(t,x) * dm and write to out
"""
function Bridge.σ!(t, x, dm, out, P::MarslandShardlowAux, xT)
    zero!(out)
    out.p .= dm*P.γ
    out
end

"""
    σt!(t, x_, y::State{Pnt}, out, P::MarslandShardlow) where Pnt
compute σ(t,x)' y, where y::State
the result is a vector of points that is written to out
"""
function σt!(t, x_, y::State{Pnt}, out, P::MarslandShardlow) where Pnt
    zero!(out)
    out .= P.γ * y.p
    out
end

"""
    Bridge.a(t,  P::Union{MarslandShardlow, MarslandShardlowAux})

Returns matrix a(t) for Marsland-Shardlow model
"""
function Bridge.a(t,  P::MarslandShardlow)
    I = Int[]
    X = UncF[]
    γ2 = P.γ^2
    for i in 1:P.n
        push!(I, 2i)
        push!(X, γ2*one(UncF))
    end
    sparse(I, I, X, 2P.n, 2P.n)
end

function Bridge.a(t,  P::MarslandShardlowAux, xT)
    I = Int[]
    X = UncF[]
    γ2 = P.γ^2
    for i in 1:P.n
        push!(I, 2i)
        push!(X, γ2*one(UncF))
    end
    sparse(I, I, X, 2P.n, 2P.n)
end

Bridge.a(t, x, P::MarslandShardlow) = Bridge.a(t, P)
Bridge.a(t, x, P::MarslandShardlowAux, xT) = Bridge.a(t, P, xT)

function Bridge.a!(t, x, A, P::MarslandShardlowAux, xT)
    zero!(A)
    G = γ2*one(UncF)
    for i ∈ (P.n+1):2P.n
        A[i,i] = G
    end
end

"""
    σ̃(t,  P::Union{MarslandShardlow, MarslandShardlowAux})

Return sparse matrix  matrix σ̃(t)
"""
function σ̃(t,  P::MarslandShardlow)
    Iind = Int[]
    Jind = Int[]
    X = UncF[]
    γ = P.γ
    @inbounds for i in 1:P.n
        push!(Iind, 2i)
        push!(Jind,i)
        push!(X, γ*one(UncF))
    end
    sparse(Iind, Jind, X, 2P.n, P.n)
end

function σ̃(t,  P::MarslandShardlowAux, xT)
    Iind = Int[]
    Jind = Int[]
    X = UncF[]
    γ = P.γ
    @inbounds for i in 1:P.n
        push!(Iind, 2i)
        push!(Jind,i)
        push!(X, γ*one(UncF))
    end
    sparse(Iind, Jind, X, 2P.n, P.n)
end

"""
    amul(t, x::State, xin::State, P::Union{MarslandShardlow, MarslandShardlowAux})

Multiply a(t,x) times xin (which is of type state)
Returns variable of type State
"""
function amul(t, x::State, xin::State, P::MarslandShardlow)
    out = copy(xin)
    zero!(out.q)
    out.p .= P.γ^2 .* xin.p
    out
end
function amul(t, x::State, xin::State, P::MarslandShardlowAux, xT)
    out = copy(xin)
    zero!(out.q)
    out.p .= P.γ^2 .* xin.p
    out
end

function amul(t, x::State, xin::Vector{<:Point}, P::MarslandShardlow)
    out = copy(x)
    zero!(out.q)
    out.p .= P.γ^2 .* vecofpoints2state(xin).p
    out
end

function amul(t, x::State, xin::Vector{<:Point}, P::MarslandShardlowAux, xT)
    out = copy(x)
    zero!(out.q)
    out.p .= P.γ^2 .* vecofpoints2state(xin).p
    out
end

################ AHS model #########################################################################################

"""
    K̄(q,τ)

K̄(q,τ) = exp(-Bridge.inner(q)/(2*τ^2))
Kernel for noisefields of AHS-model
"""
function K̄(q::T,τ) where T
     exp(-Bridge.inner(q)/(2.0*τ*τ))
end

"""
    ∇K̄(q,τ)

Gradient of kernel for noisefields
"""
function ∇K̄(q::T,τ) where T
     -1.0/(τ*τ) * K̄(q,τ) * q
end

"""
    ∇K̄(q, qT, τ)

Needed for b! in case P is auxiliary process
"""
function ∇K̄(q, qT, τ)
    -1.0/(τ*τ) * K̄(qT,τ) * q
end

"""
    Define z(q) = < ∇K̄(q - δ,τ), λ >
    Required for Stratonovich -> Ito correction in AHS-model
"""
z(q,τ,δ,λ) =  Bridge.inner(∇K̄(q - δ,τ),λ)

"""
    Define ∇z(q) = ∇ < ∇K̄(q - δ,τ), λ >
    Required for Stratonovich -> Ito correction in AHS-model
"""
∇z(q::T,τ,δ,λ) where T =  ForwardDiff.gradient(x -> z(x,τ,δ,λ),q)::T

# function for specification of diffusivity of landmarks
"""
    σq(q, nf::Noisefield) = Diagonal(nf.γ * K̄(q - nf.δ,nf.τ))

Suppose one noise field nf
Returns diagonal matrix with noisefield for position at point location q (can be vector or Point)
"""
σq(q, nf::Noisefield) = Diagonal(nf.γ * K̄(q - nf.δ,nf.τ))

"""
    σp(q, p, nf::Noisefield) = -Diagonal(p .* nf.γ .* ∇K̄(q - nf.δ,nf.τ))

Suppose one noise field nf
Returns diagonal matrix with noisefield for momentum at point location q (can be vector or Point)
"""
σp(q, p, nf::Noisefield) = -Diagonal(p .* nf.γ .* ∇K̄(q - nf.δ,nf.τ))

"""
    For AHS model compute total noise field on position experienced at a point x.
    Useful for plotting purposes.

    Example usage:
        σq(Point(0.0, 0.0), nfs)
        σq([0.0; 0.0], nfs)
"""
function σq(x, nfs::Array{<:Noisefield,1})
    out = σq(x, nfs[1])
    @inbounds for j in 2:length(nfs)
            out += σq(x, nfs[j])
    end
    out
end

σq(nfs) = (x) -> σq(x,nfs) # inefficient

"""
    construct_nfs(db::Array{Float64,1}, nfstd, γ)

Construct sequence of Noisefields for AHS model
db: domainbound. Vector [db[1], db[2]] (in case d=2), where sources are places on square grid specified by
        (-db[1]:2nfstd:db[1]) x -db[2]:2nfstd:db[2]
    nfstd: standard deviation of noise fields (the smaller: the more noise fields we use)
    γ: if set to one, then the value of the  noise field on the positions is approximately 1 at all locations in the domain
"""
function construct_nfs(db::Array{Float64,1}, nfstd, γ)
    r1 = -db[1]:2nfstd:db[1]
    if d==1
        nfloc = Point.(collect(r1))[:]
        nfscales = [2/pi*γ*Point(1.0) for x in nfloc]  # intensity
    elseif d==2
        r2 = -db[2]:2nfstd:db[2]
        nfloc = Point.(collect(product(r1, r2)))[:]
        nfscales = [2/pi*γ*Point(1.0, 1.0) for x in nfloc]  # intensity
    elseif d==3
        error("first test carefully whether construct_nfs is ok for d=3")
        r2 = -db[2]:2nfstd:db[2]
        r3 = -db[3]:2nfstd:db[3]
        nfloc = Point.(collect(product(r1, r2,r3)))[:]
        nfscales = [2/pi*γ*Point(1.0, 1.0, 1.0) for x in nfloc]  # intensity
    end
    #return Noisefield{Float64}[Noisefield(δ, λ, nfstd) for (δ, λ) in zip(nfloc, nfscales)]
    [Noisefield(δ, λ, nfstd) for (δ, λ) in zip(nfloc, nfscales)]
end

"""
Evaluate drift of landmarks in (t,x) and save to out
x is a state and out as well
"""
function Bridge.b!(t, x, out, P::Landmarks)
    zero!(out)
    pp = out.p
    qq = out.q
    @inbounds for i in 1:P.n
        for j in 1:P.n
            qq[i] += p(x,j)*kernel(q(x,i) - q(x,j), P)
            pp[i] +=  -dot(p(x,i), p(x,j)) * ∇kernel(q(x,i) - q(x,j), P)
        end
         if itostrat
             for k in 1:length(P.nfs)
                 nf = P.nfs[k]
                 K̄qi = K̄(q(x,i)-nf.δ,nf.τ)
                 zqi = z(q(x,i),nf.τ,nf.δ,nf.γ)
                 qq[i] += (0.5 * zqi * K̄qi) * nf.γ
                 pp[i] += 0.5 * dot(p(x,i),nf.γ) * ( zqi * ∇K̄(q(x,i)-nf.δ,nf.τ) -K̄qi * ∇z(q(x,i),nf.τ,nf.δ,nf.γ) )
             end
         end
    end
    out
end

"""
    Bridge.b!(t, x, out, Paux::LandmarksAux)

Evaluate drift of landmarks auxiliary process in (t,x) and save to out
x is a state and out as well
"""
function Bridge.b!(t, x, out, Paux::LandmarksAux, xT)
    zero!(out)
    pp = out.p
    qq = out.q
    @inbounds for i in 1:Paux.n
        for j in 1:Paux.n
            qq[i] += p(x,j) * Paux.gramT[i,j]
        end
        if itostrat
            qT = q(xT,i)
            p_i = p(x,i)
            for k in 1:length(Paux.nfs)
                # approximate q by qT
                nf = Paux.nfs[k]
                qq[i] += 0.5 * z(qT,nf.τ,nf.δ,nf.γ) * K̄(qT-nf.δ,nf.τ) * nf.γ
                pp[i] += 0.5 * dot(p_i,nf.γ) * ( z(qT,nf.τ,nf.δ,nf.γ) * ∇K̄(qT - nf.δ,nf.τ) -K̄(qT - nf.δ,nf.τ) * ∇z(qT,nf.τ,nf.δ,nf.γ) )
            end
        end
    end
    out
end

"""
    Bridge.B(t, Paux::LandmarksAux, xT)

Compute tildeB(t) for landmarks auxiliary process
"""
function Bridge.B(t, Paux::LandmarksAux, xT)
    X = zeros(UncF, 2Paux.n, 2Paux.n)
    @inbounds  for i in 1:Paux.n
        for j in 1:Paux.n
            X[2i-1,2j] =  Paux.gramT[i,j] * one(UncF)
        end
        if itostrat
            @inbounds for k in 1:length(Paux.nfs)
                nf = Paux.nfs[k]
                qT = q(xT,i)
                X[2i,2i] += 0.5 * ( z(qT,nf.τ,nf.δ,nf.γ) * ∇K̄(qT - nf.δ,nf.τ) -K̄(qT - nf.δ,nf.τ) * ∇z(qT,nf.τ,nf.δ,nf.γ) )  * nf.γ'
            end
        end
    end
    X
end

"""
    Bridge.B!(t,X,out, Paux::LandmarksAux, xT)

Compute B̃(t) * X (B̃ from auxiliary process) and write to out
Both B̃(t) and X are of type UncMat
"""
function Bridge.B!(t, X, out, Paux::LandmarksAux, xT)
    out .= 0.0 * out
    u = zero(UncF)
    @inbounds  for i in 1:Paux.n  # separately loop over even and odd indices
        for k in 1:2Paux.n # loop over all columns
            for j in 1:Paux.n
                 out[2i-1,k] += Paux.gramT[i,j]*X[p(j), k]
            end
            if itostrat
                u = 0.0*u
                for k in 1:length(Paux.nfs)
                    nf = P.nfs[k]
                    qT = q(xT,i)
                    u += 0.5 * ( z(qT,nf.τ,nf.δ,nf.γ) * ∇K̄(qT - nf.δ,nf.τ) -K̄(qT - nf.δ,nf.τ) * ∇z(qT,nf.τ,nf.δ,nf.γ) )  * nf.γ'
                end
                out[2i,k] = u * X[2i,k]
            end
        end
    end
    out
end

function Bridge.β(t, Paux::LandmarksAux, xT)
    out = zeros(PointF,Paux.n)
    if itostrat
        @inbounds  for i in 1:Paux.n
            for k in 1:length(Paux.nfs)
                nf = Paux.nfs[k]
                qT = q(xT,i)
                out[i] += 0.5 * z(qT,nf.τ,nf.δ,nf.γ) * K̄(qT-nf.δ,nf.τ) * nf.γ # simply take q at endpoint
            end
        end
        return(State(out,zeros(PointF,Paux.n)))
    else
        return (State(zeros(PointF,Paux.n), zeros(PointF,Paux.n)) )
    end
end

"""
    Bridge.σ!(t, x_, dm, out, P::Union{Landmarks,LandmarksAux})

Compute sigma(t,x) * dm where dm is a vector and sigma is the diffusion coefficient of landmarks
write to out which is of type State
"""
function Bridge.σ!(t, x_, dm, out, P::Landmarks)
    x = x_
    zero!(out)
    qq = out.q
    pp = out.p
    @inbounds for i in 1:P.n
        for j in 1:length(P.nfs)
            qq[i] += σq(q(x, i), P.nfs[j]) * dm[j]
            pp[i] += σp(q(x, i), p(x, i), P.nfs[j]) * dm[j]
        end
    end
    out
end

function Bridge.σ!(t, x_, dm, out, P::LandmarksAux, xT)
    x = xT
    zero!(out)
    qq = out.q
    pp = out.p
    @inbounds for i in 1:P.n
        for j in 1:length(P.nfs)
            qq[i] += σq(q(x, i), P.nfs[j]) * dm[j]
            pp[i] += σp(q(x, i), p(x, i), P.nfs[j]) * dm[j]
        end
    end
    out
end

"""
Compute sigma(t,x)' * y where y is a state and sigma is the diffusion coefficient of landmarks
returns a vector of points of length P.nfs
"""
function σtmul(t, x_, y::State{Pnt}, P::Landmarks) where Pnt
    x = x_
    out = zeros(Pnt, length(P.nfs))
    @inbounds for j in 1:length(P.nfs)
        for i in 1:P.n
            out[j] += σq(q(x, i), P.nfs[j])' * y.q[i] +
                        σp(q(x, i), p(x, i), P.nfs[j])' * y.p[i]
        end
    end
    out
end

function σtmul(t, x_, y::State{Pnt}, P::LandmarksAux, xT) where Pnt
    x = xT
    out = zeros(Pnt, length(P.nfs))
    @inbounds for j in 1:length(P.nfs)
        for i in 1:P.n
            out[j] += σq(q(x, i), P.nfs[j])' * y.q[i] +
                        σp(q(x, i), p(x, i), P.nfs[j])' * y.p[i]
        end
    end
    out
end

"""
    σ̃(t,  Paux::LandmarksAux, xT)

Return matrix σ̃(t) for LandmarksAux
"""
function σ̃(t,  Paux::LandmarksAux, xT)
    x = xT
    out = zeros(UncF, 2Paux.n, length(Paux.nfs))
    @inbounds  for i in 1:Paux.n
        for j in 1:length(Paux.nfs)
            out[2i-1,j] = σq(q(x, i), Paux.nfs[j])
            out[2i,j] = σp(q(x, i), p(x, i), Paux.nfs[j])
        end
    end
    out
end

"""
    compute σ(t,x)' y, where y::State
    the result is a vector of points that is written to out
"""
function σt!(t, x_, y::State{Pnt}, out, P::Landmarks) where Pnt
    zero!(out)
    x = x_
    @inbounds for j in 1:length(P.nfs)
         for i in 1:P.n
            out[j] += σq(q(x, i), P.nfs[j])' * q(y, i) +
                    σp(q(x, i), p(x, i), P.nfs[j])' * p(y, i)

        end
    end
    out
end

function σt!(t, x_, y::State{Pnt}, out, P::LandmarksAux, xT) where Pnt
    zero!(out)
    x = xT
    @inbounds for j in 1:length(P.nfs)
         for i in 1:P.n
            out[j] += σq(q(x, i), P.nfs[j])' * q(y, i) +
                    σp(q(x, i), p(x, i), P.nfs[j])' * p(y, i)

        end
    end
    out
end

function Bridge.a(t, x_, P::Landmarks)
    x = x_
    out = zeros(Unc{deepeltype(x)}, 2P.n,2P.n)
    @inbounds for i in 1:P.n, k in i:P.n,  j in 1:length(P.nfs)
                a11 = σq(q(x,i),P.nfs[j])
                a21 = σp(q(x,i),p(x,i),P.nfs[j])
                a12 = σq(q(x,k),P.nfs[j])
                a22 = σp(q(x,k),p(x,k),P.nfs[j])
                out[2i-1,2k-1] += a11 * a12'
                out[2i-1,2k] += a11 * a22'
                out[2i,2k-1] += a21 * a12'
                out[2i,2k] += a21 * a22'
    end
    @inbounds for i in 2:2P.n,  k in 1:i-1
            out[i,k] = out[k,i]
    end
    out
end

function Bridge.a(t, x_, P::LandmarksAux, xT)
    x = xT
    out = zeros(Unc{deepeltype(x)}, 2P.n,2P.n)
    @inbounds for i in 1:P.n, k in i:P.n,  j in 1:length(P.nfs)
                a11 = σq(q(x,i),P.nfs[j])
                a21 = σp(q(x,i),p(x,i),P.nfs[j])
                a12 = σq(q(x,k),P.nfs[j])
                a22 = σp(q(x,k),p(x,k),P.nfs[j])
                out[2i-1,2k-1] += a11 * a12'
                out[2i-1,2k] += a11 * a22'
                out[2i,2k-1] += a21 * a12'
                out[2i,2k] += a21 * a22'
    end
    @inbounds for i in 2:2P.n,  k in 1:i-1
            out[i,k] = out[k,i]
    end
    out
end

Bridge.a(t, P::LandmarksAux, xT) =  Bridge.a(t, 0, P, xT)

"""
    amul(t, x::State, xin::Vector{<:Point}, P::Union{Landmarks,LandmarksAux})

Multiply a(t,x) times a vector of points
Returns a State
(first multiply with sigma', via function σtmul, next left-multiply this vector with σ)
"""
function amul(t, x::State, xin::Vector{<:Point}, P::Landmarks)
    #vecofpoints2state(Bridge.a(t, x, P)*xin)
    out = copy(x)
    zero!(out)
    Bridge.σ!(t, x, σtmul(t, x, vecofpoints2state(xin), P),out,P)
end

function amul(t, x::State, xin::Vector{<:Point}, P::LandmarksAux, xT)
    #vecofpoints2state(Bridge.a(t, x, P)*xin)
    out = copy(x)
    zero!(out)
    Bridge.σ!(t, x, σtmul(t, x, vecofpoints2state(xin), P),out,P, xT)
end

"""
Multiply a(t,x) times a state
Returns a state
(first multiply with sigma', via function σtmul, next left-multiply this vector with σ)
"""
function amul(t, x::State, xin::State, P::Landmarks)
    out = copy(x)
    zero!(out)
    Bridge.σ!(t, x, σtmul(t, x, xin, P),out,P)
end
function amul(t, x::State, xin::State, P::LandmarksAux, xT)
    out = copy(x)
    zero!(out)
    Bridge.σ!(t, x, σtmul(t, x, xin, P),out,P, xT)
end

"""
    Bridge.a!(t, x_, out, P::Union{Landmarks,LandmarksAux})
"""
function Bridge.a!(t, x_, out, P::Landmarks)
    zero!(out)
    x = x_
    @inbounds for i in 1:P.n,  k in i:P.n, j in 1:length(P.nfs)
        a11 = σq(q(x,i),P.nfs[j])
        a21 = σp(q(x,i),p(x,i),P.nfs[j])
        a12 = σq(q(x,k),P.nfs[j])
        a22 = σp(q(x,k),p(x,k),P.nfs[j])
        out[2i-1,2k-1] += a11 * a12'
        out[2i-1,2k] += a11 * a22'
        out[2i,2k-1] += a21 * a12'
        out[2i,2k] += a21 * a22'
    end
    @inbounds for i in 2:2P.n, k in 1:i-1
            out[i,k] = out[k,i]
    end
    out
end

function Bridge.a!(t, x_, out, P::LandmarksAux, xT)
    zero!(out)
    x = xT
    @inbounds for i in 1:P.n,  k in i:P.n, j in 1:length(P.nfs)
        a11 = σq(q(x,i),P.nfs[j])
        a21 = σp(q(x,i),p(x,i),P.nfs[j])
        a12 = σq(q(x,k),P.nfs[j])
        a22 = σp(q(x,k),p(x,k),P.nfs[j])
        out[2i-1,2k-1] += a11 * a12'
        out[2i-1,2k] += a11 * a22'
        out[2i,2k-1] += a21 * a12'
        out[2i,2k] += a21 * a22'
    end
    @inbounds for i in 2:2P.n, k in 1:i-1
            out[i,k] = out[k,i]
    end
    out
end

function getγ(P)
    if isa(P,MarslandShardlow)
        out = P.γ
    elseif isa(P,Landmarks)
        out = P.nfs[1].γ[1]*pi/2
    end
    out
end

function dimwiener(P)
    if P isa MarslandShardlow
        out = P.n
    elseif P isa Landmarks
        out = length(P.nfs)
    end
    out
end

"""
    gramkernel(q, P; ϵ = 10^(-12))

Gram matrix for kernel with vector of landmarks given by q::Vector(PointF)
ϵ*I is added to avoid numerical problems that destroy PSDness of the Gram matrix
"""
function gramkernel(q, P; ϵ = 10^(-12))
    K =  [(kernel(q[i]- q[j],P) + (i==j)*ϵ) * one(UncF)   for i  in eachindex(q), j in eachindex(q)]
    PDMat(deepmat(K))
end

"""
    gram_matrix!(out, q, P)

Write into out the Gram matrix for kernel with vector of landmarks given by q::Vector(PointF)
"""
function gram_matrix!(out, q, P)
    for j ∈ 1:P.n, i ∈ j:P.n
        out[i,j] = out[j,i] = kernel(q[i]- q[j],P)
    end
    out
end

"""
    gram_matrix(q, P)

Gram matrix for kernel with vector of landmarks given by q::Vector(PointF)
"""
function gram_matrix(q, P)
    out = zeros(P.n, P.n)
    for j ∈ 1:P.n, i ∈ j:P.n
        out[i,j] = out[j,i] = kernel(q[i]- q[j],P)
    end
    out
end

"""
    hamiltonian(x::NState, P::MarslandShardlow)

Returns Hamiltonian at `x` for deterministic system (no Wiener noise)
"""
function hamiltonian(x::NState, P::MarslandShardlow)
    s = 0.0
    for i in 1:P.n, j in 1:P.n
        s += dot(x.p[i], x.p[j])*kernel(x.q[i] - x.q[j], P)
    end
    0.5 * s
end


text:
using Test, Random, CUDA, Knet, AutoGrad, Statistics

using CUDA.CUDNN:
    cudnnNormalizationForward,
    cudnnNormalizationForward!,
    cudnnNormalizationForwardInference,
    cudnnNormalizationForwardTraining,
    cudnnNormalizationBackward,
    cudnnActivationDescriptor,
    cudnnNormMode_t,
        CUDNN_NORM_PER_ACTIVATION, # 0, bnScale, bnBias tensor dims are 1xCxHxWx.. (one value per CHW...-slice, normalized over N slice)
        CUDNN_NORM_PER_CHANNEL,    # 1, bnScale, bnBias tensor dims are 1xCx1x1 (one value per C-dim normalized over Nx1xHxW subtensors)
    cudnnNormOps_t,
        CUDNN_NORM_OPS_NORM,                # 0, /* do normalization only */
        CUDNN_NORM_OPS_NORM_ACTIVATION,     # 1, /* do Norm, then activation */
        CUDNN_NORM_OPS_NORM_ADD_ACTIVATION, # 2, /* do Norm, then elemWiseAdd, then activation */
    cudnnNormAlgo_t,
        CUDNN_NORM_ALGO_STANDARD, # 0
        CUDNN_NORM_ALGO_PERSIST,  # 1
    cudnnActivationMode_t,
        CUDNN_ACTIVATION_SIGMOID,      # 0
        CUDNN_ACTIVATION_RELU,         # 1
        CUDNN_ACTIVATION_TANH,         # 2
        CUDNN_ACTIVATION_CLIPPED_RELU, # 3
        CUDNN_ACTIVATION_ELU,          # 4
        CUDNN_ACTIVATION_IDENTITY,     # 5
    cudnnNanPropagation_t,
        CUDNN_NOT_PROPAGATE_NAN, # 0
        CUDNN_PROPAGATE_NAN,     # 1
    cudnnTensorFormat_t,
        CUDNN_TENSOR_NCHW,        # 0, /* row major (wStride = 1, hStride = w) */
        CUDNN_TENSOR_NHWC,        # 1, /* feature maps interleaved ( cStride = 1 )*/
        CUDNN_TENSOR_NCHW_VECT_C, # 2, /* each image point is vector of element of C, vector length in data type */
    handle

@testset "cudnn/normalization" begin

    x, s, b, s2, b2, s3, b3 = (Param(CUDA.randn(Float64,x...)) for x in ((5,4,3,2),(1,1,3,1),(1,1,3,1),(5,4,3,1),(5,4,3,1),(5,1,1,1),(5,1,1,1)))

    # @gcheck needs Inference and Training to compute the same function
    # For this they need to use the same mean/var
    # @gcheck sums the output of non-scalar functions by default
    # This has the effect of multiplying scale with a 0 mean array
    # So we'll use getindex(y,1) instead
    function normtest(x,s,b; dims=(1,2,4), o...)
        m = mean(value(x);dims)
        v = var(value(x);mean=m,corrected=false,dims)
        y = cudnnNormalizationForward(x,m,v,b,s; training=Knet.training(), o...)
        y[1]
    end

    @test normtest(x,s,b) ≈ value(@diff(normtest(x,s,b)))
    @test @gcheck normtest(x,s,b)
    @test @gcheck normtest(x,s2,b2; mode = CUDNN_NORM_PER_ACTIVATION, dims=4)
    @test @gcheck normtest(x,s,b; algo = CUDNN_NORM_ALGO_PERSIST)
    @test @gcheck normtest(x,s3,b3; algo = CUDNN_NORM_ALGO_PERSIST, format = CUDNN_TENSOR_NHWC, dims=(2,3,4))
    @test @gcheck normtest(x,s,b; alpha = 2)
    @test @gcheck normtest(x,s,b; epsilon = 0)
    @test @gcheck normtest(x,s3,b3; format = CUDNN_TENSOR_NHWC, dims=(2,3,4))
    @test @gcheck normtest(x,s,b; exponentialAverageFactor = 0.01)
    @test @gcheck normtest(x,s,b; savedMean = similar(s), savedInvVariance = similar(s))
    # cudnn-8.0.5: Currently, CUDNN_NORM_OPS_NORM_ACTIVATION and CUDNN_NORM_OPS_NORM_ADD_ACTIVATION are not supported.
    #@test_skip @gcheck normtest(x,s,b; normOps = CUDNN_NORM_OPS_NORM_ACTIVATION, activationMode = CUDNN_ACTIVATION_RELU)
    #@test_skip @gcheck normtest(x,s,b; normOps = CUDNN_NORM_OPS_NORM_ADD_ACTIVATION, activationMode = CUDNN_ACTIVATION_RELU, z)
    #@test_skip @gcheck normtest(x,s,b; groupCnt = 2) # Currently only groupCnt=1 is supported
    #@test_skip @gcheck normtest(x,s,b; beta = 2) # we don't use beta in functional code
end


text:
cd(@__DIR__)
using Pkg
CI = get(ENV, "CI", nothing) == "true" || get(ENV, "GITHUB_TOKEN", nothing) !== nothing
CI && Pkg.activate(@__DIR__)
CI && Pkg.instantiate()
using DrWatson
using Documenter, DataFrames, Parameters, Dates, JLD2, UnPack
using DocumenterTools: Themes

# %%
# download the themes
for file in ("juliadynamics-lightdefs.scss", "juliadynamics-darkdefs.scss", "juliadynamics-style.scss")
    download("https://raw.githubusercontent.com/JuliaDynamics/doctheme/master/$file", joinpath(@__DIR__, file))
end
# create the themes
for w in ("light", "dark")
    header = read(joinpath(@__DIR__, "juliadynamics-style.scss"), String)
    theme = read(joinpath(@__DIR__, "juliadynamics-$(w)defs.scss"), String)
    write(joinpath(@__DIR__, "juliadynamics-$(w).scss"), header*"\n"*theme)
end
# compile the themes
Themes.compile(joinpath(@__DIR__, "juliadynamics-light.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-light.css"))
Themes.compile(joinpath(@__DIR__, "juliadynamics-dark.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-dark.css"))

isdir(datadir()) && rm(datadir(); force = true, recursive = true)

using Literate
Literate.markdown("src/workflow.jl", "src"; credit = false)

makedocs(modules = [DrWatson, UnPack],
    sitename= "DrWatson",
    authors = "George Datseris and contributors.",
    doctest = false,
    format = Documenter.HTML(
        prettyurls = CI,
        assets = [
            "assets/logo.ico",
            asset("https://fonts.googleapis.com/css?family=Quicksand|Montserrat|Source+Code+Pro|Lora&display=swap", class=:css),
            ],
        ),
    pages = [
        "Introduction" => "index.md",
        "DrWatson Workflow Tutorial" => "workflow.md",
        "Project Setup" => "project.md",
        "Naming Simulations" => "name.md",
        "Saving Tools" => "save.md",
        "Running & Listing Simulations" => "run&list.md",
        "Real World Examples" => "real_world.md"
        ],
)

if CI
    deploydocs(
        repo = "github.com/JuliaDynamics/DrWatson.jl.git",
        target = "build",
        push_preview = true
    )
end


text:
using Sockets
using DirectSockets
using Test

@testset "DirectSockets" begin
  sock1 = DirectSockets.UDPSocket()
  sock2 = DirectSockets.UDPSocket()
  @test sock1 isa DirectSockets.UDPSocket
  @test sock2 isa DirectSockets.UDPSocket
  @test bind(sock1, ip"0.0.0.0", 9000)
  @test send(sock2, ip"127.0.0.1", 9000, UInt8[1,2,3,4])
  @test recv(sock1) == UInt8[1,2,3,4]
  @test send(sock2, ip"127.0.0.1", 9000, UInt8[5,6,7,8])
  buf = zeros(UInt8, DirectSockets.BUFSIZE)
  n = recv!(sock1, buf)
  @test n == 4
  @test buf[1:4] == UInt8[5,6,7,8]
  @test send(sock2, ip"127.0.0.1", 9000, UInt8[9,10,11,12])
  host_port, data = recvfrom(sock1)
  @test host_port isa Sockets.InetAddr{IPv4}
  @test data == UInt8[9,10,11,12]
  @test send(sock1, host_port, UInt8[13,14,15,16])
  @test recv(sock2) == UInt8[13,14,15,16]
  @test close(sock1) === nothing
  @test close(sock2) === nothing
  sock3 = DirectSockets.UDPSocket()
  @test bind(sock3, Sockets.InetAddr(ip"0.0.0.0", 9000))
  @test close(sock3) === nothing
  @test_throws ErrorException recv!(sock3, buf)
  @test_throws ErrorException recv(sock3)
end


text:
# Interface -------

function contributors(owner, repo, attempts = 3; auth = AnonymousAuth(), options...)
    contributors(auth, owner, repo, attempts; options...)
end

function contributors(auth::Authorization, owner, repo, attempts = 3; headers = Dict(), options...)
    r = attempt_stats_request(auth, owner, repo, "contributors", attempts; headers = headers, options...)

    data = JSON.parse(r.data)
    map!(data) do udata
        udata["author"] = User(udata["author"])
        udata
    end
end

function commit_activity(owner, repo, attempts = 3; auth = AnonymousAuth(), options...)
    commit_activity(auth, owner, repo, attempts; options...)
end

function commit_activity(auth::Authorization, owner, repo, attempts = 3; headers = Dict(), options...)
    r = attempt_stats_request(auth, owner, repo, "commit_activity", attempts; headers = headers, options...)
    JSON.parse(r.data)
end

function code_frequency(owner, repo, attempts = 3; auth = AnonymousAuth(), options...)
    code_frequency(auth, owner, repo, attempts; options...)
end

function code_frequency(auth::Authorization, owner, repo, attempts = 3; headers = Dict(), options...)
    r = attempt_stats_request(auth, owner, repo, "code_frequency", attempts; headers = headers, options...)
    JSON.parse(r.data)
end

function participation(owner, repo, attempts = 3; auth = AnonymousAuth(), options...)
    participation(auth, owner, repo, attempts; options...)
end

function participation(auth::Authorization, owner, repo, attempts = 3; headers = Dict(), options...)
    r = attempt_stats_request(auth, owner, repo, "participation", attempts; headers = headers, options...)
    JSON.parse(r.data)
end

function punch_card(owner, repo, attempts = 3; auth = AnonymousAuth(), options...)
    punch_card(auth, owner, repo, attempts; options...)
end

function punch_card(auth::Authorization, owner, repo, attempts = 3; headers = Dict(), options...)
    r = attempt_stats_request(auth, owner, repo, "punch_card", attempts; headers = headers, options...)
    JSON.parse(r.data)
end

# Utility -------

function attempt_stats_request(auth, owner, repo, stat, attempts; headers = Dict(), options...)
    authenticate_headers(headers, auth)

    for a in attempts:-1:0
        r = get(URI(API_ENDPOINT; path = "/repos/$owner/$repo/stats/$stat");
                headers = headers,
                options...)

        handle_error(r)

        if r.status == 200
            return r
        end

        sleep(2.0)
    end

    throw(StatsError("Unsuccessfully attempted to retrieve $stat $attempts times."))
end


text:
## Seeded PageRank animation
# We are going to make an animation where we move from vertex to vertex
# following the nearest vertex in a seeded PageRank simulation

using MatrixNetworks
using Plots

function seeded_pagerank_animation(A,xy,seed,α,nsteps)
  curseed = seed
  started = Set{Int}() # this is the list of points we've started.
  @assert(nsteps <= size(A,1))
  @assert(is_connected((A)))
  anim = @animate for i=1:nsteps
    push!(started, curseed)
    f2 = seeded_pagerank(A, α, curseed)
    p = sortperm(f2) # this gets things to draw in the right z-order
    scatter(xy[p,1],xy[p,2], marker_z=-abs.(log10.(f2[p])).^0.5,colorbar=false,alpha=0.5,
      color=:magma,
      size=(1200,1200),
      markerstrokecolor=nothing,markerstrokewidth=0,label="",framestyle=:none,background=:black)
    curseed = -1
    for v in sortperm(f2,rev=true)
      if v in started
        continue
      else
        curseed = v
        break
      end
    end
  end
  return anim
end
A = load_matrix_network_all("geoclusters-example").A
xy = load_matrix_network_all("geoclusters-example").xy
anim = seeded_pagerank_animation(A,xy,1,0.5,100)
gif(anim,"geoclusters-animation.gif",fps=10)

## Now use the actual artistsim graph
names = readlines("../2-unit-1-demos/artistsim.names")
data = readdlm("../2-unit-1-demos/artistsim.smat")
xy = readdlm("artistsim.xy")
A = sparse(Int.(data[2:end,1]).+1,
           Int.(data[2:end,2]).+1,
           Int.(data[2:end,3]),
           Int(data[1,1]),Int(data[1,1]))
Asym = max.(A,A') # make the graph undirected
Acc,ccfilt = largest_component(Asym)
xycc = xy[ccfilt,:]

##
anim = seeded_pagerank_animation(Acc,xycc,200,0.5,200)
gif(anim,"artistsim-seeded-pagerank-animation.gif",fps=10)


text:
function singular(cond)
  quote let
    local expected = $(string(cond))
    $(esc(cond)) || error("$expected, turns out to be false")
  end end
end

function exact(cond,prediction)
  quote let
    local got = $(esc(cond))
    local expected = $(esc(prediction))
    local isPass = typeof(expected) <: Function ? expected(got): expected == got
    isPass || error("got $(repr(got)), expected to be $(repr(expected))")
  end end
end

function singleConditiontest(cond,leftMost,prediction)
  quote let
    local condition = $(esc(cond))
    local expected = $(esc(prediction))
    local fn = $(esc(leftMost))
    fn(condition, expected) || error("got $(repr(condition)), expected to be $fn $(repr(expected))")
  end end
end

function multipleConditionTest(cond,prediction,rightMostFn,fns)
  quote let
    local condition = $(esc(cond))
    local expected= $(esc(prediction))
    local rightMostFn = $(esc(rightMostFn))
    local fns = $(Expr(:vcat,reverse(fns)...))
    local allTrue = rightMostFn(condition,expected)
    local fn
    for fn in fns
      allTrue = fn(allTrue)
    end
    allTrue || error("got $(repr(condition)), expected to be $(join(fns,' ')) $rightMostFn $(repr(expected))")
  end end
end

macro is(args...)
  len = length(args)
  test, = args
  cond, leftMost = test.args
  if len == 1
    if test.head == :(=>)
      exact(test.args...)
    else
      singular(test)
    end
  elseif len > 2
    multipleConditionTest(cond,args[end], args[end - 1], [leftMost, [args[i] for i in 2:(len- 2)]])
  else
    singleConditiontest(cond,leftMost,args[end])
  end
end

not(x) = x == false
not(x,y) = x != y

truthy(x) = x != false
falsy(x) = !truthy(x)
falsey = falsy

above(x,y) = x > y
above{T<:Array}(array::T,x) = length(array) > x
isAbove = isabove = above

least(x,y) = x >= y
least{T<:Array}(array::T,x) = length(array) >= x
atleast = atLeast = least

below(x,y) = x < y
below{T<:Array}(array::T,x) = below(array) < x
isbelow = isBelow = below

most(x,y) = x <= y
most{T<:Array}(array::T,x) = below(array) <= x
ismost = isMost = most

empty(x) = length(x) == 0


text:
export SendDesc, WaitDesc, send, wait, signal

const Cond = Int32

struct SendDesc
    dbrec::ISem32
    db::ISem64
end

struct WaitDesc
    cond::GDS_WAIT_COND
    sema::ISem32
    flag::ISem32
end

@inline function send(info::SendDesc)
    release(ind.dbrec)
    threadfence_system()
    release(info.db)
end

@inline function wait(info::WaitDesc)
    wait(info.sema, info.cond)
end

@inline function signal(info::WaitDesc)
    release(info.flag)
end


text:
module PandorasBoxColoredWaterfall

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/coloredWaterfall" ColoredWaterfall(x::Integer, y::Integer, color="LightSkyBlue")

segmentCache = Dict{String, Any}()

const colors = sort(collect(keys(Ahorn.XNAColors.colors)))

const waterSegmentMatrix = [
    1 0 0 0 0 0 0 1 1 1;
    1 0 0 0 0 0 0 1 1 1;
    1 0 0 0 0 0 0 1 1 1;
    1 0 0 0 0 0 0 1 1 1;
    1 0 0 0 0 0 0 1 1 1;
    1 0 0 0 0 0 0 1 1 1;
    1 0 0 0 0 0 0 1 1 1;
    1 0 0 0 0 0 0 1 1 1;
    1 1 0 0 0 0 0 0 1 1;
    1 1 0 0 0 0 0 0 1 1;
    1 1 1 0 0 0 0 0 0 1;
    1 1 1 0 0 0 0 0 0 1;
    1 1 1 0 0 0 0 0 0 1;
    1 1 1 0 0 0 0 0 0 1;
    1 1 1 0 0 0 0 0 0 1;
    1 1 1 0 0 0 0 0 0 1;
    1 1 1 0 0 0 0 0 0 1;
    1 1 1 0 0 0 0 0 0 1;
    1 1 0 0 0 0 0 0 1 1;
    1 1 0 0 0 0 0 0 1 1;
    1 1 0 0 0 0 0 0 1 1;
]

function getColor(color)
    if haskey(Ahorn.XNAColors.colors, color)
        return Ahorn.XNAColors.colors[color]

    else
        try
            return ((Ahorn.argb32ToRGBATuple(parse(Int, replace(color, "#" => ""), base=16))[1:3] ./ 255)..., 1.0)

        catch

        end
    end

    return (1.0, 1.0, 1.0, 1.0)
end

function getSegment(color)
    if haskey(segmentCache, color)
        return segmentCache[color]

    else
        baseColor = getColor(color)
        fillColor = baseColor .* 0.3
        surfaceColor = baseColor .* 0.8

        waterSegment = Ahorn.matrixToSurface(
            waterSegmentMatrix,
            [
                fillColor,
                surfaceColor
            ]
        )

        segmentCache[color] = waterSegment

        return waterSegment
    end
end

function getHeight(entity::ColoredWaterfall, room::Maple.Room)
    waterEntities = filter(e -> e.name == "water", room.entities)
    waterRects = Ahorn.Rectangle[
        Ahorn.Rectangle(
            Int(get(e.data, "x", 0)),
            Int(get(e.data, "y", 0)),
            Int(get(e.data, "width", 8)),
            Int(get(e.data, "height", 8))
        ) for e in waterEntities
    ]

    width, height = room.size
    x, y = Int(get(entity.data, "x", 0)), Int(get(entity.data, "y", 0))
    tx, ty = floor(Int, x / 8) + 1, floor(Int, y / 8) + 1

    wantedHeight = 8 - y % 8
    while wantedHeight < height - y
        rect = Ahorn.Rectangle(x, y + wantedHeight, 8, 8)

        if any(Ahorn.checkCollision.(waterRects, Ref(rect)))
            break
        end

        if get(room.fgTiles.data, (ty + 1, tx), '0') != '0'
            break
        end

        wantedHeight += 8
        ty += 1
    end

    return wantedHeight
end

const placements = Ahorn.PlacementDict(
    "Waterfall (Pandora's Box)" => Ahorn.EntityPlacement(
        ColoredWaterfall
    )
)

Ahorn.editingOptions(entity::ColoredWaterfall) = Dict{String, Any}(
    "color" => colors
)

function Ahorn.selection(entity::ColoredWaterfall, room::Maple.Room)
    x, y = Ahorn.position(entity)
    height = getHeight(entity, room)

    return Ahorn.Rectangle(x, y, size(waterSegmentMatrix, 2), height)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::ColoredWaterfall, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))
    color = get(entity.data, "color", "LightSkyBlue")

    height = getHeight(entity, room)
    segment = getSegment(color)
    segmentHeight, segmentWidth = size(waterSegmentMatrix)

    Ahorn.Cairo.save(ctx)

    Ahorn.rectangle(ctx, 0, 0, segmentWidth, height)
    Ahorn.clip(ctx)

    for i in 0:segmentHeight:ceil(Int, height / segmentHeight) * segmentHeight
        Ahorn.drawImage(ctx, segment, 0, i)
    end

    Ahorn.restore(ctx)
end

end

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
# This file demonstrates a relatively simple way to implement the AbstractTrees
# interface for a classic binary tree. See "binarytree_infer.jl" for a more performant
# (but laborious) approach.

if !isdefined(@__MODULE__, :BinaryNode)
    include("binarytree_core.jl")
end

## Things we need to define
function AbstractTrees.children(node::BinaryNode)
    if isdefined(node, :left)
        if isdefined(node, :right)
            return (node.left, node.right)
        end
        return (node.left,)
    end
    isdefined(node, :right) && return (node.right,)
    return ()
end

## Things that make printing prettier
AbstractTrees.printnode(io::IO, node::BinaryNode) = print(io, node.data)

## Optional enhancements
# These next two definitions allow inference of the item type in iteration.
# (They are not sufficient to solve all internal inference issues, however.)
Base.eltype(::Type{<:TreeIterator{BinaryNode{T}}}) where T = BinaryNode{T}
Base.IteratorEltype(::Type{<:TreeIterator{BinaryNode{T}}}) where T = Base.HasEltype()

## Let's test it. First build a tree.
root = BinaryNode(0)
l = leftchild(1, root)
r = rightchild(2, root)
lr = rightchild(3, l)

print_tree(root)
collect(PostOrderDFS(root))
@static if isdefined(@__MODULE__, :Test)
    @testset "binarytree_easy.jl" begin
        @test [node.data for node in PostOrderDFS(root)] == [3, 1, 2, 0]
        @test [node.data for node in PreOrderDFS(root)] == [0, 1, 3, 2]
        @test [node.data for node in Leaves(root)] == [3, 2]
    end
end


text:
const maxdegree = 2
const CONST_TYPE = Float32


text:
#---------------------------------------------------------
# # [Ideal age](@id ideal-age)
#---------------------------------------------------------

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/1_ideal_age.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/1_ideal_age.ipynb)

#md # !!! note
#md #     All the AIBECS tutorials and how-to guides are available as Jupyter notebooks.
#md #     You can view them with [nbviewer](https://nbviewer.jupyter.org/)or execute them online with [binder](https://mybinder.org/) by clicking on the badges above!
#md #     (Note that binder can be slow to launch and its memory caps can be a problem when running.)

# The tracer equation for the ideal age is
#
# $$\left(\partial_t + \mathbf{T}\right) \boldsymbol{a} = 1 - \frac{\boldsymbol{a}}{τ} \, (\boldsymbol{z} \le z_0),$$
#
# where the sink term on the right clamps the age to $0$ at the surface (where $\boldsymbol{z} \le z_0$).
# The smaller the timescale $\tau$, the quicker $\boldsymbol{a}$ is restored to $0$ at the surface.

# AIBECS can interpret tracer equations as long as you arrange them under the generic form:
#
# $$\big(\partial_t + \mathbf{T}(\boldsymbol{p}) \big) \boldsymbol{x} = \boldsymbol{G}(\boldsymbol{x}, \boldsymbol{p}),$$
#
# where $\mathbf{T}(\boldsymbol{p})$ is the transport, $\boldsymbol{G}(\boldsymbol{x}, \boldsymbol{p})$ is the net local sources and sinks, and $\boldsymbol{p}$ is the vector of model parameters.
# We will then use the AIBECS to simulate the ideal age by finding the steady-state of the system, i.e., the solution of
#
# $$\partial_t \boldsymbol{x} = \boldsymbol{F}(\boldsymbol{x}, \boldsymbol{p}) = \boldsymbol{G}(\boldsymbol{x}, \boldsymbol{p}) - \mathbf{T}(\boldsymbol{p}) \, \boldsymbol{x} = 0.$$

# In this tutorial, we will simulate the ideal age by
# 1. defining functions for `T(p)` and `G(x,p)`,
# 1. defining the parameters `p`,
# 1. generating the state function `F(x,p)` and solving the associated steady-state problem,
# 1. and finally making a plot of our simulated ideal age.

# We start by telling Julia that we want to use the AIBECS package and the OCIM2 circulation
#nb # (the Ocean Circulation Inverse Model, see [1](https://doi.org/10.1029/2018JC014716) for details).
#md # (the Ocean Circulation Inverse Model[^1]).

#md # [^1]:
#md #     DeVries, T., & Holzer, M. (2019). Radiocarbon and helium isotope constraints on deep ocean ventilation and mantle‐³He sources. Journal of Geophysical Research: Oceans, 124, 3036–3057. doi:[10.1029/2018JC014716](https://doi.org/10.1029/2018JC014716)

using AIBECS
grd, TOCIM2 = OCIM2.load()

#md # !!! note
#md #     If it's your first time, Julia will ask you to download the OCIM2, in which case you should accept (i.e., type `y` and "return").
#md #     Once downloaded, AIBECS will remember where it downloaded the file and it will only load it from your laptop.
#nb # > **Note**
#nb # > If it's your first time, Julia will ask you to download the OCIM2, in which case you should accept (i.e., type `y` and "return").
#nb # > Once downloaded, AIBECS will remember where it downloaded the file and it will only load it from your laptop.

# `grd` is an `OceanGrid` object containing information about the 3D grid of the OCIM2 circulation and `TOCIM2` is the transport matrix representing advection and diffusion.

# The local sources and sinks for the age take the form

function G(x,p)
    @unpack τ, z₀ = p
    return @. 1 - x / τ * (z ≤ z₀)
end

# as per the tracer equation.
# The `@unpack` line unpacks the parameters `τ` and `z₀`.
# The `return` line returns the net sources and sinks.
# (The `@.` "macro" tells Julia that the operations apply to every element.)

# We can define the vector `z` of depths with `depthvec`.

z = depthvec(grd)

# Now we must construct a type for `p` the parameters.
# This type must contain our parameters `τ` and `z₀`.

struct IdealAgeParameters{U} <: AbstractParameters{U}
    τ::U
    z₀::U
end

# The type is now ready for us to generate an instance of the parameter `p`.
# Let's use `τ = 1.0` (s) and `z₀` the minimum depth of the model.

p = IdealAgeParameters(1.0, 30.0)

# We now use the AIBECS to generate the state function $\boldsymbol{F}$ (and its Jacobian) via

F = AIBECSFunction(TOCIM2, G)

# Now that `F` and `p` are defined, we are going to solve for the steady-state.
# But first, we must create a `SteadyStateProblem` object that contains `F`, `p`, and an initial guess `x_init` for the age.
# (`SteadyStateProblem` comes from [SciMLBase](https://github.com/JuliaDiffEq/SciMLBase.jl).)

# Let's make a vector of 0's for our initial guess.

nb = sum(iswet(grd))  # number of wet boxes
x_init = zeros(nb)    # Start with age = 0 everywhere

# Now we can create our `SteadyStateProblem` instance

prob = SteadyStateProblem(F, x_init, p)

# And finally, we can `solve` this problem, using the AIBECS `CTKAlg()` algorithm,

age = solve(prob, CTKAlg())

# This should take a few seconds.

# To conclude this tutorial, let's have a look at the age using AIBECS' plotting recipes and [Plots.jl](https://github.com/JuliaPlots/Plots.jl).

using Plots

# We first convert the age in years
# (because the default SI unit we used, i.e., seconds, is a bit small relative to global ocean timescales).

age_in_yrs = age * u"s" .|> u"yr"

# And we take a horizontal slice at about 2000m.

plothorizontalslice(age_in_yrs, grd, depth=2000u"m", color=:magma)

# Or look at the horiontal mean

plothorizontalmean(age_in_yrs, grd)

# That's it for this tutorial...
# Good job!


text:
# cannot go over the speed limit
# cannot go backward
function speed_limit_actions(mdp::GenerativeMergingMDP, s::AugScene)
    valid_acts = Int64[]
    sizehint!(valid_acts, 7)
    ego = get_by_id(s.scene, EGO_ID)
    acc = s.ego_info.acc
    for a in actions(mdp)
        next_v = ego.state.v + action_map(mdp, acc, a).a*mdp.dt
        if 0 <= next_v <= mdp.env.main_lane_vmax
            push!(valid_acts, a)
        end
    end
    return valid_acts
end

# jerk constraints
function acceleration_limit_actions(mdp::GenerativeMergingMDP, s::AugScene)
    acc = s.ego_info.acc
    valid_acts = Int64[]
    # if we were hard braking we must release the brake
    # if isapprox(acc, mdp.max_deceleration)
    #     for a in actions(mdp)
    #         if action_map(mdp, 0.0, a).a >= 0.0
    #             push!(valid_acts, a)
    #         end
    #     end
    #     return valid_acts
    # else
    for a in actions(mdp)
        next_acc = action_map(mdp, acc, a).a
        if abs(next_acc) <= mdp.comfortable_acceleration
            push!(valid_acts, a)
        end
    end
    return valid_acts
    # end
end

function collision_lookahead(mdp::GenerativeMergingMDP, s::AugScene)
    valid_acts = Int64[]
    action_list = LatLonAccel[LatLonAccel(0.,0.0) for veh in s.scene if veh.id != EGO_ID]
    for a in actions(mdp)
        sp = deepcopy(s.scene)
        acc = action_map(mdp, s.ego_info.acc, a)
        tick!(sp, mdp.env.roadway, [acc, action_list...], mdp.dt)
        if !icollision_checker(sp, EGO_ID)
            push!(valid_acts, a)
        end
    end
    return valid_acts
end

function idm_limit_actions(mdp::GenerativeMergingMDP, s::AugScene)
    scene = s.scene
    ego_ind = findfirst(EGO_ID, scene)
    ego = scene[ego_ind]
    acc = s.ego_info.acc
    fore_res = get_front_neighbor(mdp.env, s.scene, EGO_ID)
    if fore_res.ind == nothing
        return actions(mdp)
    end
    v_ego = scene[ego_ind].state.v
    v_oth = scene[fore_res.ind].state.v
    headway = fore_res.Δs
    track_longitudinal!(mdp.ego_idm, v_ego, v_oth, headway)
    des_acc = mdp.ego_idm.a
    valid_acts = Int64[]
    for a in actions(mdp)
        next_acc = action_map(mdp, acc, a)
        if next_acc.a <= des_acc
            push!(valid_acts, a)
        end
    end
    return valid_acts
end

function safe_actions(mdp::GenerativeMergingMDP, s::AugScene)
    speed = speed_limit_actions(mdp, s)
    # jerk = acceleration_limit_actions(mdp, s)
    collision_free = collision_lookahead(mdp, s)
    available_acts = intersect(speed, collision_free)
    # available_acts = intersect(speed, jerk, collision_free)
    # available_acts = intersect(jerk, speed, collision_free)
    # ego = get_by_id(s.scene, EGO_ID)
    # if get_lane(mdp.env.roadway, ego) == main_lane(mdp.env)
    #     idm = idm_limit_actions(mdp, s)
    #     available_acts = intersect(speed, jerk, idm, collision_free)
    # end
    # push!(available_acts, HARD_BRAKE)
    # if isempty(available_acts)
        # return Int64[HARD_BRAKE]
    # end
    sort!(available_acts)
    return available_acts
end

function safe_actions(mdp::MDP{S,A}, o::AbstractArray) where {S,A}
    s = POMDPs.convert_s(S, o, mdp)
    return safe_actions(mdp, s)
end

function best_action(acts::Vector{A}, val::AbstractArray{T}, problem::M) where {A, T <: Real, M <: Union{POMDP, MDP}}
    best_a = acts[1]
    best_ai = actionindex(problem, best_a)
    best_val = val[best_ai]
    for a in acts
        ai = actionindex(problem, a)
        if val[ai] > best_val
            best_val = val[ai]
            best_ai = ai
            best_a = a
        end
    end
    return best_a::A
end

function masked_linear_epsilon_greedy(max_steps::Int64, eps_fraction::Float64, eps_end::Float64)
    # define function that will be called to select an action in DQN
    # only supports MDP environments
    function action_masked_epsilon_greedy(policy::AbstractNNPolicy, env::MDPEnvironment, obs, global_step::Int64, rng::AbstractRNG)
        eps = DeepQLearning.update_epsilon(global_step, eps_fraction, eps_end, max_steps)
        acts = safe_actions(policy.problem, obs) #XXX using pomdp global variable replace by safe_actions(mask, obs)
        val = actionvalues(policy, obs) #change this
        if rand(rng) < eps
            return (rand(rng, acts), eps)
        else
            return (best_action(acts, val, env.problem), eps)
        end
    end
    return action_masked_epsilon_greedy
end

function masked_evaluation()
    function masked_evaluation_policy(policy::AbstractNNPolicy, env::MDPEnvironment, n_eval::Int64, max_episode_length::Int64, verbose::Bool)
        avg_r = 0
        violations = 0
        avg_steps = 0
        for i=1:n_eval
            done = false
            r_tot = 0.0
            step = 0
            obs = reset!(env)
            DeepQLearning.resetstate!(policy)
            while !done && step <= max_episode_length
                acts = safe_actions(env.problem, obs)
                val = actionvalues(policy, obs)
                act = best_action(acts, val, env.problem)
                obs, rew, done, info = step!(env, act)
                r_tot += rew
                step += 1
            end
            # println("Episode reward: $r_tot")
            r_tot < 0. ? violations += 1 : nothing
            avg_r += r_tot
            avg_steps += step
        end
        # log_value(LOGGER, "eval_reward", avg_r/n_eval)
        # log_value(LOGGER, "collisions", violations/n_eval*100)
        # log_value(LOGGER, "avg_steps", avg_steps/n_eval)
        if verbose
            @printf("Evaluation ... Avg Reward %1.2f | Violations (%%) %2.2f | Avg Steps %3.2f", avg_r/n_eval, violations/n_eval*100, avg_steps/n_eval)
        end
        avg_r /= n_eval
        return  avg_r
    end
    return masked_evaluation_policy
end

struct MaskedPolicy{P<:AbstractNNPolicy} <: AbstractNNPolicy
    policy::P
end

function POMDPs.action(p::MaskedPolicy, s)
    acts = safe_actions(p.policy.problem, s)
    val = actionvalues(p.policy, s)
    act = best_action(acts, val, p.policy.problem)
    @assert in(act, acts)
    return act
end

struct MaskedRandomPolicy <: Policy
    problem::GenerativeMergingMDP
    rng::MersenneTwister
end

function POMDPs.action(p::MaskedRandomPolicy, s)
    acts = safe_actions(p.problem, s)
    return rand(p.rng, acts)
end


text:
export tap_on_unsubscribe, TapBeforeUnsubscription, TapAfterUnsubscription

import Base: show

"""
    TapBeforeUnsubscription

One of the strategies for `tap_on_unsubscribe` operator. With `TapBeforeUnubscription` tap callback will be called before actual unsubscription.

See also: [`tap_on_unsubscribe`](@ref), [`TapAfterUnsubscription`](@ref)
"""
struct TapBeforeUnsubscription end

"""
    TapAfterUnsubscription

One of the strategies for `tap_on_unsubscribe` operator. With `TapBeforeUnsubscription` tap callback will be called after actual unsubscription.

See also: [`tap_on_unsubscribe`](@ref), [`TapBeforeUnsubscription`](@ref)
"""
struct TapAfterUnsubscription end

"""
    tap_on_unsubscribe(tapFn::F, strategy::S = TapBeforeUnsubscription()) where { F <: Function }

Creates a tap operator, which performs a side effect on the unsubscription on the source Observable, but return an Observable that is identical to the source. Tap callback triggers only once.

# Arguments
- `tapFn::Function`: side-effect tap function with `() -> Nothing` signature
- `strategy`: (optional), specifies the order of a side-effect and an actual unsubscription, uses `TapBeforeUnsubscription` by default

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

# Examples
```jldoctest
using Rocket

source = from([ 1, 2, 3 ])
subscription = subscribe!(source |> tap_on_unsubscribe(() -> println("Someone unsubscribed")), logger())
unsubscribe!(subscription)
;

# output

[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Data: 3
[LogActor] Completed
Someone unsubscribed
```

See also: [`TapBeforeUnsubscription`](@ref), [`TapAfterUnsubscription`](@ref), [`tap`](@ref), [`tap_on_subscribe`](@ref), [`tap_on_complete`](@ref), [`logger`](@ref)
"""
tap_on_unsubscribe(tapFn::F, strategy::S = TapBeforeUnsubscription()) where { F <: Function, S } = TapOnUnsubscribeOperator{F, S}(tapFn, strategy)

struct TapOnUnsubscribeOperator{F, S} <: InferableOperator
    tapFn    :: F
    strategy :: S
end

operator_right(::TapOnUnsubscribeOperator, ::Type{L}) where L = L

function on_call!(::Type{L}, ::Type{L}, operator::TapOnUnsubscribeOperator{F, S}, source) where { L, F, S }
    return proxy(L, source, TapOnUnsubscribeProxy{F, S}(operator.tapFn, operator.strategy))
end

struct TapOnUnsubscribeProxy{F, S} <: SourceProxy
    tapFn    :: F
    strategy :: S
end

source_proxy!(::Type{L}, proxy::TapOnUnsubscribeProxy{F, T}, source::S) where { L, S, F, T } = TapOnUnsubscribeSource{L, S, F, T}(proxy.tapFn, proxy.strategy, source)

@subscribable struct TapOnUnsubscribeSource{L, S, F, T} <: Subscribable{L}
    tapFn    :: F
    strategy :: T
    source   :: S
end

mutable struct TapOnUnsubscribeSubscription{S, F, T} <: Teardown
    is_unsubscribed :: Bool
    tapFn           :: F
    strategy        :: T
    subscription    :: S
end

as_teardown(::Type{ <: TapOnUnsubscribeSubscription }) = UnsubscribableTeardownLogic()

function on_subscribe!(source::TapOnUnsubscribeSource{L, S, F, T}, actor) where { L, S, F, T }
    return TapOnUnsubscribeSubscription(false, source.tapFn, source.strategy, subscribe!(source.source, actor))
end

function on_unsubscribe!(subscription::TapOnUnsubscribeSubscription)
    result = __on_unsubscribe_with_tap(subscription.strategy, subscription)
    subscription.is_unsubscribed = true
    return result
end

function __on_unsubscribe_with_tap(::TapBeforeUnsubscription, subscription::TapOnUnsubscribeSubscription)
    if !subscription.is_unsubscribed
        subscription.tapFn()
    end
    return unsubscribe!(subscription.subscription)
end

function __on_unsubscribe_with_tap(::TapAfterUnsubscription, subscription::TapOnUnsubscribeSubscription)
    result = unsubscribe!(subscription.subscription)
    if !subscription.is_unsubscribed
        subscription.tapFn()
    end
    return result
end

Base.show(io::IO, ::TapOnUnsubscribeOperator)          = print(io, "TapOnUnsubscribeOperator()")
Base.show(io::IO, ::TapOnUnsubscribeProxy)             = print(io, "TapOnUnsubscribeProxy()")
Base.show(io::IO, ::TapOnUnsubscribeSource{L}) where L = print(io, "TapOnUnsubscribeSource($L)")
Base.show(io::IO, ::TapOnUnsubscribeSubscription)      = print(io, "TapOnUnsubscribeSubscription()")


text:
abstract type Var end   # variable

struct DiscreteVariable{N,D,B,M} <: Var
    name::N
    domain::D
    booleans::B
    varmap::M  # forward dictionary from domain to variables
end

# Base.show(io::IO, x::DiscreteVariable) = print(io, "Variable $(x.name) ∈ $(x.domain)")
Base.show(io::IO, x::DiscreteVariable) = print(io, "variable $(x.name)")

DiscreteVariable(name, D) = DiscreteVariable(name, sort(collect(D)))

# index integers as themselves; otherwise just 1..n:
indices(v::Vector{Int}) = v
indices(v) = collect(eachindex(v))

function DiscreteVariable(name, domain::Vector)
    booleans = [Variable(name, i) for i ∈ indices(domain)]
    varmap = Dict(i => v for (i, v) in zip(domain, booleans))

    return DiscreteVariable(name, domain, booleans, varmap)
end

not_both(x, y) = ¬x ∨ ¬y
implies(x, y) = ¬x ∨ y
and(x, y) = ¬(¬x ∨ ¬y)

"Given a set of boolean variables, gives clauses specifying that exactly one of them is true"
function exactly_one(vars)
    clauses = [∨(vars...)]   # one must be true

    # two cannot be true:
    for i in 1:length(vars)
        for j in i+1:length(vars)
            push!(clauses, not_both(vars[i], vars[j]))
        end
    end

    return clauses
end

clauses(x::DiscreteVariable) = exactly_one(x.booleans)

domain(x::Var) = x.domain
# vars(x::DiscreteVariable) = [x]

function decode(varmap, x::Var)
    values = [varmap[v] for v in x.booleans]

    num_true = count(values)

    if num_true ≠ 1
        error("Variable $x has not been successfully solved: values $(x.vars .=> values)")
    end

    return x.domain[findfirst(values)]
end


text:
abstract type AbstractTransform end

struct Output{T}
    z::Vector{T}
    rz::Vector{T}
end

function display(o::Output)
    println("    z: Transformed Observations")
    display(o.z)
    println("\n    rz: Reduced Transformed Observations")
    display(o.rz)
    return nothing
end

function output(W, g, λ, y, p)
    z = (x -> transform(g, λ, x)).(y)
    @views rz = z[p+1:end] .- W' * z[1:p]
    return Output(z, rz)
end
function output(T, n, m)
    z = Matrix{T}(undef, n)
    rz = Matrix{T}(undef, m)
    return Output(z, rz)
end
function set!(o::Output, W, g, λ, y)
    p = length(o.z) - length(o.rz)
    o.z .= (x -> radial(g, λ, x)).(y)
    @views o.rz .= o.z[:, p+1:end] .- W * o.z[:, 1:p]
    return o
end

logjacobian(g::AbstractTransform, λ::Tuple, y) = sum(x -> log∂(g, λ, x), y)
transform(g::AbstractTransform, λ::Tuple, y) = transform(g, λ..., y)
log∂(g::AbstractTransform, λ::Tuple, y) = log∂(g, λ..., y)

struct Identity <: AbstractTransform end

transform(::Identity, y) = y
log∂(::Identity, _) = 0

struct BoxCox <: AbstractTransform end

transform(::BoxCox, λ::Number, y) = λ == 0 ? log(y) : expm1(log(y) * λ) / λ
log∂(::BoxCox, λ::Number, y) = (λ - 1) * log(y)

struct YeoJohnson <: AbstractTransform end

function transform(::YeoJohnson, λ::Number, y)
    λ′ = y >= 0 ? λ : 2 - λ
    return sign(y) * (λ′ == 0 ? log1p(abs(y)) : expm1(log1p(abs(y)) * λ′) / λ′)
end
log∂(::YeoJohnson, λ::Number, y) = sign(y) * (λ - 1) *  log(abs(y))

struct SumTanh <: AbstractTransform end

function transform(::SumTanh, α::Vector, β::Vector, c::Vector, y)
    return sum(α .* tanh.(β .* (y .+ c)))
end
function log∂(::SumTanh, α::Vector, β::Vector, c::Vector, y)
    return sum(α .* β .* sech.(β .* (y .+ c)) .^ 2)
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
##########################################################
###################### BUILD DEVICE ######################
using IntegratedOptics
######################
# Plot Options
ϵ_axes = Axes(title = "'ϵᵣ'",xlabel = "'x'", ylabel = "'y'", palette = :tempo, auto = "fix",size =" ratio -1")
Ez_axes = Axes(title = "'Ez'",xlabel = "'x'", ylabel = "'y'", palette = :balance, auto = "fix",size =" ratio -1")
######################
######################## MATERIALS ########################
Si =   Material(ϵᵣ = 3.64^2 )
SiO2 = Material(ϵᵣ = 1.45^2 )
Fg = Material(ϵᵣ = .8*(1.45^2 + 3.46^2) )
##########################################################
######################## SIMULATION #######################
λ₀ = 1.55
Lx , Ly =  λ₀ .* (4 , 3)
n_max = 3.64
n_min = 1.45
dx , dy = (λ₀/18/n_max) .* (1, 1)
sim = Simulation(λ₀ = 1.55 ;  grid = Grid(extent = (Lx,Ly) , spacing =  (dx,dy) ))
sim(setpml!,1.0)
##########################################################
###################### BUILD DEVICE ######################
#--------------------------------------------------------#
sim(setbackground!,SiO2)
#--------------------------------------------------------#
w_d = 1.2
h_d = 1.5
design_region = Cuboid(((.5(Lx - w_d) , .5(Ly - h_d)),(.5(Lx + w_d) , .5(Ly + h_d))))
sim(setmaterial!,Fg, design_region)
#--------------------------------------------------------#
h_wg_in = 0.5
input_waveguide  = Cuboid(((0 , .5(Ly - h_wg_in)),(.45(Lx) , .5(Ly + h_wg_in))))
sim(setmaterial!,Si,input_waveguide)
#--------------------------------------------------------#
h_wg_out = 0.5
output_waveguide  =  Cuboid(((.55Lx , .5(Ly - h_wg_out)),(Lx , .5(Ly + h_wg_out))))
sim(setmaterial!,Si,output_waveguide)
#--------------------------------------------------------#
##########################################################
sim(HeatMap,:ϵᵣ,ẑ, axes = ϵ_axes)
#sim(saveplot , output = "./device.png")

################### PREPARE TARGET FIELD #################
### INPUT
########################################################
x_cut = 1.2
sim_wg_in = sim(slice, x̂ , x_cut )
eig_vals , eig_vecs = sim_wg_in(solve_for_modes,solver = EigArpack())
mode_te_in = sim_wg_in(extractreshape, eig_vecs[:,1] ,x̂)
#mode_tm_in = sim_wg(Photon.extractreshape, eig_vecs[:,2] ,ŷ)
sim_wg_in(LinePlot,mode_te_in,lw = 2,title = "'mode profile'")
input_mode_sim = deepcopy(sim)
input_mode_sim(setbackground!,SiO2)
input_mode_sim(setmaterial!,Si,Cuboid(((0 , .5(Ly - h_wg_in)),(Lx , .5(Ly + h_wg_in)))))
input_mode_sim(HeatMap,:ϵᵣ,ẑ, axes = ϵ_axes)
input_mode_sim(attachmode!,mode_te_in, input_mode_sim(:J,ẑ), x̂ , x_cut )
input_mode_sim(HeatMap,:J,ẑ)
input_mode_sim(solve!)
input_mode_sim(HeatMap,:E,ẑ, axes = Ez_axes)
###########################################
## OUTPUT
##########################################################
x_cut = 1.2
sim_wg_out = sim(slice, x̂ , x_cut )
eig_vals , eig_vecs = sim_wg_out(solve_for_modes)
mode_te_out = sim_wg_out(extractreshape, eig_vecs[:,3] ,x̂)
sim_wg_out(LinePlot,mode_te_out,lw = 2,title = "'mode profile'")
output_mode_sim = deepcopy(sim)
output_mode_sim(setbackground!,SiO2)
output_mode_sim(setmaterial!,Si,Cuboid(((0 , .5(Ly - h_wg_in)),(Lx , .5(Ly + h_wg_in)))))
output_mode_sim(HeatMap,:ϵᵣ,ẑ,axes =ϵ_axes)
output_mode_sim(attachmode!,mode_te_out, output_mode_sim(:J,ẑ), x̂ , x_cut )
output_mode_sim(HeatMap,:J,ẑ)
output_mode_sim(solve!)
output_mode_sim(HeatMap,:E,ẑ,axes = Ez_axes)
#########################################################

target_x = zeros(ComplexF64,size(sim.grid))
target_y = zeros(ComplexF64,size(sim.grid))
target_z = zeros(ComplexF64,size(sim.grid))

sim(add!,target_z ,input_mode_sim(:E,ẑ),  x-> .4Lx > x[1] > x_cut)
sim(add!,target_z ,output_mode_sim(:E,ẑ), x-> .6Lx < x[1] <= Lx-x_cut)
sim(HeatMap,target_z, axes = Ez_axes)

sim(add!,sim(:J,ẑ) ,input_mode_sim(:E,ẑ), EVERYWHERE)
#sim(attachmode!,mode_te_in, sim(:J,ẑ), x̂ , x_cut )
#sim(HeatMap,  :J, ẑ)
sim(setTFSF!,x̂,LOW,x_cut + dx)
#sim.ϵᵣ_zz = input_mode_sim.ϵᵣ_zz
sim(solve!)
#############################################################
# Toplogy Optimization
#######################################
# Optimization Settings
obj =TargetObjective(target_x,target_y,target_z)
optimizer = ADAGrad(1)
using ImageFiltering
kernel = Kernel.gaussian(3)
ϵ_limits = (1.45^2,3.64^2)
Nt = 100
binarization_rate(i) =  0.5* i * -(-(ϵ_limits...))/Nt
###########################################
# Optimization Loop
for i in 1:Nt
    f , ∇f  =  f_∇f(obj,sim,solver = LU())
    sim(update!,optimizer, sim(:ϵᵣ,ẑ), imfilter(∇f[ẑ], kernel),design_region)

    sim(boxconstraint!,sim(:ϵᵣ,ẑ),ϵ_limits,design_region)
    sim(gradualbinarize!,sim(:ϵᵣ,ẑ),ϵ_limits,binarization_rate(i),design_region)
   ####
    iseven(i) ? display(sim(HeatMap,  :E, ẑ,axes =Ez_axes )) : display(sim(HeatMap,  :ϵᵣ, ẑ,axes =ϵ_axes ))
    trans = sum(real.(sim(:S,x̂,((x̂,4.5,4.5),(ŷ,1.5,3))))) / sum(real.(sim(:S,x̂,((x̂,2,2),(ŷ,1.5,3)))))
    println("iter:$(i)",": ",trans)
   ####
end

max_stars_repo_path: examples/integrated optics/two ports/spatial mode profile converters/TO_spatial_mode_converter.jl


text:
#### Riemann Matrix
import Base: size
export Riemann
"""
`Riemann` matrix

Riemann matrix is defined as `A = B[2:N+1, 2:N+1]`, where
`B[i,j] = i-1` if `i` divides `j`, and `-1` otherwise.
[Riemann hypothesis](http://en.wikipedia.org/wiki/Riemann_hypothesis) holds
if and only if `det(A) = O( N! N^(-1/2+epsilon))` for every `epsilon > 0`.

```julia
julia> Riemann(7)
7x7 Riemann{Int64}:
  1  -1   1  -1   1  -1   1
 -1   2  -1  -1   2  -1  -1
 -1  -1   3  -1  -1  -1   3
 -1  -1  -1   4  -1  -1  -1
 -1  -1  -1  -1   5  -1  -1
 -1  -1  -1  -1  -1   6  -1
 -1  -1  -1  -1  -1  -1   7
```

For more details see Friedrich Roesler, Riemann's hypothesis as an eigenvalue problem,
Linear Algebra and its Applications, Vol. 81, (1986)
http://www.sciencedirect.com/science/article/pii/0024379586902557
"""
struct Riemann{Int} <: AbstractMatrix{Int}
    n::Int
end # immutable

# Define its size

size(A::Riemann, dim::Integer) = A.n
size(A::Riemann)= size(A,1), size(A,1)

# Index into a Riemann
function getindex(A::Riemann,i::Integer,j::Integer)
#    return (i+1)%(j+1)==0 ? i : -1
    if i<=j && (j+1)%(i+1)==0
        return i
    else
        return -1
    end
end # getindex

# Dense version of Riemann
Matrix(A::Riemann) =[A[i,j] for i=1:size(A,1), j=1:size(A,2)]


text:
using DashBootstrapComponents, DashHtmlComponents

function make_popover(placement)
    return dbc_popover(
        [dbc_popoverheader("Header"), dbc_popoverbody("This is a $placement popover")],
        id = "popover-$placement",
        target = "popover-$placement-target",
        placement = placement,
        is_open = false,
    )
end;

function make_button(placement)
    return dbc_button(
        "Popover on $placement",
        id = "popover-$placement-target",
        className = "mx-2",
        n_clicks = 0,
    )
end;

popovers = html_div(
    [
        [make_button(p) for p in ["top", "left", "right", "bottom"]]
        [make_popover(p) for p in ["top", "left", "right", "bottom"]]
    ],
);

for p in ["top", "left", "right", "bottom"]
    callback!(
        app,
        Output("popover-$p", "is_open"),
        Input("popover-$p-target", "n_clicks"),
        State("popover-$p", "is_open"),
    ) do n, is_open
        return n > 0 ? is_open == 0 : is_open
    end
end;

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
module test_julia_bottom

using Test

@test Base.Bottom === Union{}
@test Core.TypeofBottom === typeof(Union{})

end # module test_julia_bottom