text:
@testset "size" begin
    x = Axis(UnitSRange(1, 3), UnitSRange(1, 3))
    #@test StaticRanges.Size(typeof(x)) === StaticRanges.Size{(3,)}()

    @test size(x) == (3,)
    @test size(x, 1) == 3
end

@testset "insert!" begin
    v = AxisVector{Int}()
    push!(v, 1)
    insert!(v, 1, 1)
    insert!(v, 1, 1.0)
    @test v == [1, 1, 1]

    v = AxisVector([1,2,3])
    @test keys.(axes(@inferred(reverse(v)))) == (3:-1:1,)

    v = AxisVector([1,2,3], ["a", "b", "c"])
    @test keys.(axes(@inferred(reverse(v)))) == (["c", "b", "a"],)

    # TODO should this be a strict key array still?
    v = OffsetVector([1,2,3], 1)
    @test keys.(axes(@inferred(reverse(v)))) == (4:-1:2,)

    @test_throws MethodError insert!(AxisArray(1:2), 2, 2)
end

@testset "AxisArray constructors" begin

    A_fixed = ones(2, 2)

    # TODO properly import StaticArrays
    A_static = SMatrix{2,2}(A_fixed)

    @testset "AxisArray(::AbstractArray)" begin
        A_fixed_axes = @inferred(AxisArray(A_fixed));
        A_static_axes = @inferred(AxisArray(A_static));

        @test @inferred(all(i -> length(i) == 2, axes(A_fixed_axes)))
        @test @inferred(all(i -> known_length(i) == 2, axes(A_static_axes)))
    end

    @testset "AxisArray(::AbstractArray, ::Tuple{Keys...})" begin
        A_fixed_axes = @inferred(AxisArray(A_fixed, (["a", "b"], [:one, :two])));
        A_static_axes = @inferred(AxisArray(A_static));

        @test @inferred(all(i -> length(i) == 2, axes(A_fixed_axes)))
        @test @inferred(all(i -> known_length(i) == 2, axes(A_static_axes)))
    end

    @testset "AxisArray(::Array{T,0})" begin
        A = AxisArray(Array{Int,0}(undef, ()))
        @test A isa AxisArray{Int,0}
    end

    # FIXME
    @testset "AxisArray{T,N}(::AbstractArray...)" begin
        @test parent_type(@inferred(AxisArray{Int,2}(ones(2,2), (["a", "b"], [:one, :two])))) <: Array{Int,2}
        @test parent_type(@inferred(AxisArray{Int,2}(ones(2,2), ["a", "b"], [:one, :two]))) <: Array{Int,2}
        # TODO why would I need these constructors?
        #@test parent_type(@inferred(AxisArray{Int,2}(ones(2,2), (2, 2)))) <: Array{Int,2}
        #@test parent_type(@inferred(AxisArray{Int,2}(ones(2,2), 2, 2))) <: Array{Int,2}

    end

    @testset "AxisArray(undef, ::Tuple{Keys...})" begin
        @test parent_type(@inferred(AxisArray{Int}(undef, (["a", "b"], [:one, :two])))) <: Array{Int,2}
        @test parent_type(@inferred(AxisArray{Int}(undef, ["a", "b"], [:one, :two]))) <: Array{Int,2}
        @test parent_type(@inferred(AxisArray{Int}(undef, (2, 2)))) <: Array{Int,2}
        @test parent_type(@inferred(AxisArray{Int}(undef, 2, 2))) <: Array{Int,2}

        @test parent_type(@inferred(AxisArray{Int,2}(undef, (["a", "b"], [:one, :two])))) <: Array{Int,2}
        @test parent_type(@inferred(AxisArray{Int,2}(undef, ["a", "b"], [:one, :two]))) <: Array{Int,2}
        @test parent_type(@inferred(AxisArray{Int,2}(undef, (2, 2)))) <: Array{Int,2}
        @test parent_type(@inferred(AxisArray{Int,2}(undef, 2, 2))) <: Array{Int,2}
    end

    @testset "AxisArray{T,N,P}" begin
        A = AxisArray(reshape(1:4, 2, 2))
        @test typeof(@inferred(convert(AxisArray{Int32,2,Array{Int32,2}}, A))) <: AxisArray{Int32,2,Array{Int32,2}}
    end

    @testset "collect(::AxisArray)" begin
        x = AxisArray(reshape(1:10, 2, 5))
        @test typeof(parent(collect(x))) <: typeof(collect(parent(x)))
        @test x isa AxisArray
    end
end

@testset "matrices" begin
    @testset "matmul" begin
    @testset "Matrix-Matrix" begin
            a = AxisArray(ones(2, 3), (3:4, 1:3));
            b = AxisArray(ones(3, 2), (2:4, 2:3));

            @testset "standard case" begin
                @test a * b == 3ones(2, 2)
                @test keys.(axes(a * b)) == (3:4, 2:3)

                # TODO errors
                @test ones(4, 3) * b == 3ones(4, 2)
                @test keys.(axes(ones(4, 3) * b)) == (1:4, 2:3)

                @test a * ones(3, 7) == 3ones(2, 7)
                @test keys.(axes(a * ones(3,7))) == (3:4, 1:7)
            end
        end

        m = AxisArray(ones(1, 1), (Axis(2:2), Axis(3:3),));
        v = AxisArray(ones(1), (Axis(4:4),));

        @testset "Matrix-Vector" begin
            @test m * v == ones(1)
            @test keys.(axes(m * v)) == (2:2,)
        end

        @testset "Vector-Matrix" begin
            @test v * m == ones(1, 1)
            @test keys.(axes(v * m)) == (4:4, 3:3)
        end

        @testset "Vector-Vector" begin
            v = [1, 2, 3]
            av = AxisArray(v, (Axis(2:4),))
            @test_throws MethodError av * av
            @test av' * av == 14
            @test av' * av == adjoint(av) * v == transpose(av) * v
            @test av' * av == adjoint(v) * av == transpose(v) * av
            @test av * av' == [1 2 3; 2 4 6; 3 6 9]
        end

        @testset "Transpose/Adjoint" begin
            A = ones(3, 3)
            @test A * AxisArray(ones(3)) == A * ones(3)
            @test A' * AxisArray(ones(3)) == A' * ones(3)
            @test A' * AxisArray(ones(3)) == A' * ones(3)
            @test A' * AxisArray(ones(3,3)) == A' * ones(3, 3)
            @test ones(3)' * A == AxisArray(ones(3))' * A
            @test ones(1, 3) * A == AxisArray(ones(1, 3)) * A
            @test AxisArray(ones(3)) * ones(3)' == ones(3) * ones(3)'

            @test transpose(A) * AxisArray(ones(3)) == transpose(A) * ones(3)
            @test transpose(A) * AxisArray(ones(3)) == transpose(A) * ones(3)
            @test transpose(A) * AxisArray(ones(3,3)) == transpose(A) * ones(3, 3)
            @test transpose(ones(3)) * A == transpose(AxisArray(ones(3))) * A
            @test AxisArray(ones(3)) * transpose(ones(3)) == ones(3) * transpose(ones(3))
        end

        @testset "Hermitian" begin
            A = Hermitian([1 0 2+2im 0 3-3im; 0 4 0 5 0; 6-6im 0 7 0 8+8im; 0 9 0 1 0; 2+2im 0 3-3im 0 4])
            @test A * AxisArray(ones(5)) == A * ones(5)
            @test A' * AxisArray(ones(5)) == A' * ones(5)
            @test A' * AxisArray(ones(5)) == A' * ones(5)
            @test A' * AxisArray(ones(5,5)) == A' * ones(5, 5)
            @test AxisArray(ones(5,5)) * A' == ones(5,5) * A'
            @test ones(5)' * A == AxisArray(ones(5))' * A
            @test ones(1, 5) * A == AxisArray(ones(1, 5)) * A

            @test AxisArray(ones(5,5)) * transpose(A) == ones(5,5) * transpose(A)
            @test transpose(A) * AxisArray(ones(5)) == transpose(A) * ones(5)
            @test transpose(A) * AxisArray(ones(5)) == transpose(A) * ones(5)
            @test transpose(A) * AxisArray(ones(5,5)) == transpose(A) * ones(5, 5)
            @test transpose(ones(5)) * A == transpose(AxisArray(ones(5))) * A
        end

        @testset "UpperTriangular" begin
            A = UpperTriangular([1.0 2.0 3.0;
                                4.0 5.0 6.0;
                                7.0 8.0 9.0])
            @test A * AxisArray(ones(3)) == A * ones(3)
            @test A' * AxisArray(ones(3)) == A' * ones(3)
            @test A' * AxisArray(ones(3)) == A' * ones(3)
            @test A' * AxisArray(ones(3,3)) == A' * ones(3, 3)
            @test ones(3)' * A == AxisArray(ones(3))' * A
            @test ones(1, 3) * A == AxisArray(ones(1, 3)) * A

            @test transpose(A) * AxisArray(ones(3)) == transpose(A) * ones(3)
            @test transpose(A) * AxisArray(ones(3)) == transpose(A) * ones(3)
            @test transpose(A) * AxisArray(ones(3,3)) == transpose(A) * ones(3, 3)
            @test transpose(ones(3)) * A == transpose(AxisArray(ones(3))) * A
        end

    end

end

@testset "Base" begin
    m = [10 20; 31 40]
    am = AxisArray(m, (2:3, 3:4))

    @testset "$f" for f in (sum, prod, maximum, minimum, extrema)
        @test f(am) == f(m)
        @test f(am; dims=1) == f(m; dims=1)

        @test keys.(axes(f(am; dims=1))) == (2:2, 3:4)
    end

    @testset "$f" for f in (cumsum, cumprod, sort)
        @test f(am; dims=1) == f(m; dims=1)

        @test keys.(axes(f(am; dims=1))) == (2:3, 3:4)

        @test f([1, 4, 3]) == f(AxisArray([1, 4, 3]))
    end

    #= TODO
    @testset "sort!" begin
        a = [1 9; 7 3]
        am = AxisArray(a, (2:3, 3:4))

        # Vector case
        veca = [1, 9, 7, 3]
        sort!(AxisArray(veca, :vec); order=Base.Reverse)
        @test issorted(veca; order=Base.Reverse)

        # Higher-dim case: `dims` keyword in `sort!` requires Julia v1.1+
        if VERSION > v"1.1-"
            sort!(nda, dims=:y)
            @test issorted(a[2, :])
            @test_throws UndefKeywordError sort!(nda)

            sort!(nda; dims=:x, order=Base.Reverse)
            @test issorted(a[:, 1]; order=Base.Reverse)
        end
    end
    =#

    @testset "$f!" for (f,f!) in zip((sum, prod, maximum, minimum), (sum!, prod!, maximum!, minimum!))
        a = [10 20; 31 40]
        am = AxisArray(a, (2:3, 3:4)) # size (2,2)

        a1 = sum(am, dims=1)           # size (1,2)
        a2 = sum(am, dims=2)           # size (2,1)
        @testset "ndims==2" begin
            @test f!(a1, am) == f!(a1, a) == f(a, dims=1)
            @test f!(a2, am) == f!(a2, a) == f(a, dims=2)

            @test keys.(axes(f!(a1, am))) == (2:2, 3:4) == keys.(axes(f!(a1, a)))
            @test keys.(axes(f!(a2, am))) == (2:3, 3:3) == keys.(axes(f!(a2, a)))
        end
        @testset "ndims==1 too" begin
            x = AxisArray([3, 4], (2:3,))
            y = AxisArray([5, 6], (3:4,))
            v = AxisArray([7, 8])

            @test f!(parent(x), am) == f!([0,0], a) == dropdims(f(a, dims=2), dims=2)
            @test f!(v, am) == f!(x, a)
            @test f!(y', am) == f!([0 0], am) == f(a, dims=1)
        end
    end

    # FIXME
    #=
    @testset "eachslice" begin
        if VERSION > v"1.1-"
            slices = [[111 121; 211 221], [112 122; 212 222]]
            cat_slices = cat(slices...; dims=3)
            a = AxisArray(cat_slices, (2:3, 3:4, 4:5))

            @test sum(eachslice(a; dims=3)) ==
                  sum(eachslice(cat_slices; dims=3)) ==
                  slices[1] + slices[2]
            #@test_throws ArgumentError eachslice(nda; dims=(1, 2))
            #@test_throws ArgumentError eachslice(a; dims=(1, 2))

            #@test_throws UndefKeywordError eachslice(nda)
            #@test_throws UndefKeywordError eachslice(cat_slices)

            @test axes_keys(first(eachslice(a; dims=2))) == (2:3, 4:5)
        end
    end
    =#

    @testset "mapslices" begin
        m = [10 20; 31 40]
        am = AxisArray(m, (2:3, 3:4))

        @test mapslices(join, am; dims=1) ==
              mapslices(join, m; dims=1) ==
              ["1031" "2040"]

        @test mapslices(join, am; dims=2) ==
              mapslices(join, m, dims=2) ==
              reshape(["1020", "3140"], Val(2))

        @test mapslices(join, am; dims=(1, 2)) ==
              mapslices(join, m; dims=(1, 2)) ==
              reshape(["10312040"], (1, 1))

       @test keys.(axes(mapslices(join, am; dims=2))) == (2:3, 3:3)
    end

    @testset "mapreduce" begin
        m = [10 20; 31 40]
        maxes = AxisArray(m, (2:3, 3:4))

        @test mapreduce(isodd, |, maxes) == true == mapreduce(isodd, |, m)
        @test mapreduce(isodd, |, maxes; dims=1) == [true false]
        @test mapreduce(isodd, |, maxes; dims=2) == [false true]'
        @test keys.(axes(mapreduce(isodd, |, maxes; dims=2))) == (2:3, 3:3)
    end

    @testset "zero" begin
        m = [10 20; 31 40]
        maxes = AxisArray(m, (2:3, 3:4))

        @test zero(maxes) == [0 0; 0 0] == zero(m)
        @test keys.(axes(zero(maxes))) == (2:3, 3:4)
    end

    @testset "count" begin
        m = [true false; true true]
        maxes = AxisArray(m, (2:3, 3:4))

        @test count(maxes) == count(m) == 3
        @test_throws Exception count(a; dims=1)
    end

    # TODO test warnings for immutable axes
    @testset "push!, pop!, etc" begin
        v = AxisArray([10, 20, 30], (Axis(UnitMRange(2, 4)),))

        @test length(push!(v, 40)) == 4
        @test keys.(axes(pushfirst!(v, 0))) == (1:5,)
        @test v == [0, 10, 20, 30, 40]

        @test pop!(v) == 40
        @test popfirst!(v) == 0
        @test v == [10, 20, 30]
    end

    @testset "append!, empty!" begin
        v = AxisArray([10, 20, 30], (UnitMRange(2, 4),))
        v45 = AxisArray([40, 50], (UnitMRange(3, 4),))
        v0 = AxisArray([0, 0], (UnitMRange(4, 5),))

        append!(v, v45)
        @test length(v) == 5
        append!(v, [60,70])
        @test keys(axes(v, 1)) == 2:8

        #@test_throws DimensionMismatch append!(ndv, ndv0)
        @test v == 10:10:70 # error was thrown before altering

        @test keys(axes(empty!(v), 1)) == 2:1
        @test length(v) == 0
    end

    @testset "map, map!" begin
        m = AxisArray([11 12; 21 22], (2:3, 3:4))

        @test keys.(axes(map(+, m, m, m))) == (2:3, 3:4)
        @test keys.(axes(map(+, m, parent(m), m))) == (2:3, 3:4)
        @test keys.(axes(map(+, parent(m), m))) == (2:3, 3:4)

        # this method only called based on first two arguments:
        #@test dimnames(map(+, parent(maxes), parent(maxes), maxes)) == (:_, :_)

        # one-arg forms work without adding anything... except on 1.0...
        @test keys.(axes(map(sqrt, m))) == (2:3, 3:4)
        @test foreach(sqrt, m) === nothing

        #= TODO is this something we actually want?
        # map! may return a different wrapper of the same data, like sum!
        semi = AxisArray(rand(2,2), (:x, :_))
        @test dimnames(map!(sqrt, rand(2,2), maxes)) == (:x, :y)
        @test dimnames(map!(sqrt, semi, maxes)) == (:x, :y)
        =#

        zed = similar(m, Float64)
        @test map!(sqrt, zed, m) == sqrt.(m)
        @test zed[begin,begin] == sqrt(m[begin,begin])
    end

    #= TODO
    @testset "filter" begin
        nda = AxisArray([11 12; 21 22], (:x, :y))
        ndv = AxisArray(1:7, (:z,))

        @test dimnames(filter(isodd, ndv)) == (:z,)
        @test dimnames(filter(isodd, nda)) == (:_,)
    end
    =#

    #= TODO
    @testset "collect(generator)" begin
        nda = AxisArray([11 12; 21 22], (:x, :y))
        ndv = AxisArray([10, 20, 30], (:z,))

        @test dimnames([sqrt(x) for x in nda]) == (:x, :y)

        @test dimnames([x^i for (i,x) in enumerate(ndv)]) == (:z,)
        @test dimnames([x^i for (i,x) in enumerate(nda)]) == (:x, :y)

        # Iterators.product -- has all names
        @test dimnames([x+y for x in nda, y in ndv]) == (:x, :y, :z)
        @test dimnames([x+y for x in nda, y in 1:5]) == (:x, :y, :_)
        @test dimnames([x+y for x in 1:5, y in ndv]) == (:_, :z)
        four = [x*y/z^p for p in 1:2, x in ndv, y in 1:2, z in nda]
        @test dimnames(four) == (:_, :z, :_, :x, :y)

        # Iterators.flatten -- no obvious name to use
        @test dimnames([x+y for x in nda for y in ndv]) == (:_,)

        if VERSION >= v"1.1"
            # can't see inside eachslice generators, until:
            # https://github.com/JuliaLang/julia/pull/32310
            @test dimnames([sum(c) for c in eachcol(nda)]) == (:_,)
        end
    end
    =#

    @testset "equality" begin
        a = AxisArray([10 20; 30 40])
        a2 = AxisArray([10 20; 30 40])
        a3 = AxisArray([10 20; 30 40])
        a4 = AxisArray([11 22; 33 44])
        v = AxisArray([10, 20, 30])

        @testset "$eq" for eq in (Base.:(==), isequal, isapprox)
            @test eq(a, a) == eq(parent(a), a) == eq(a, parent(a))
            @test eq(a, a2)
            @test eq(a, a3)
            @test !eq(a, a4)
        end
        @test isapprox(a, a4; atol=2π)
        @test isapprox(a, parent(a4); atol=2π)
        @test isapprox(parent(a), a4; atol=2π)
        @test AxisArray(1:2) == SimpleAxis(1:2)
    end

end  # Base

@testset "Statistics" begin
    m = [10 20; 30 40]
    maxes = AxisArray(m, (2:3, 3:4))
    @testset "$f" for f in (mean, std, var, median)
        @test f(maxes) == f(m)
        @test f(maxes; dims=1) == f(m; dims=1)

        @test keys.(axes(f(maxes; dims=1))) == (2:2, 3:4)
    end
end

@testset "Array Interface" begin
    x = AxisArray([1 2; 3 4]);
    @test parent_type(typeof(x)) == parent_type(x) == Array{Int,2}
    @test size(x) == (2, 2)
    @test parent(x) == parent(parent(x))

    @test x[CartesianIndex(2,2)] == 4
    x[CartesianIndex(2,2)] = 5
    @test x[CartesianIndex(2,2)] == 5

    @test eltype(similar(x, Float64, (2:3, 4:5))) <: Float64
    @test_throws DimensionMismatch AxisArray(rand(2,2), (2:9,2:1))

    x = AxisArray(reshape(1:8, 2, 2, 2));
    @test @inferred(x[CartesianIndex(1,1), 1]) == 1
    @test @inferred(x[CartesianIndex(1,1), :]) == parent(parent(x)[CartesianIndex(1,1), :])
    @test @inferred(x[:, CartesianIndex(1,1)]) == parent(parent(x)[:, CartesianIndex(1,1)])
    @test @inferred(x[[true,true], CartesianIndex(1,1)]) == parent(parent(x)[[true,true], CartesianIndex(1,1)])
end

@testset "I/O" begin
    io = IOBuffer()
    x = AxisArray([1 2; 3 4])
    write(io, x)

    seek(io, 0)
    y = AxisArray(Array{Int}(undef, 2, 2))
    read!(io, y)

    @test y == x
end

@testset "resize!" begin
    x = AxisArray(ones(3))
    resize!(x, 2)
    @test x == [1, 1]
end

@testset "view" begin
    A = [1 2; 3 4]
    Aaxes = AxisArray(A, ["a", "b"], 2.0:3.0)

    A_view = @inferred(view(A, :, 1))
    Aaxes_view = @inferred(view(Aaxes, :, 1))
    @test A_view == Aaxes_view

    fill!(A_view, 0)
    fill!(Aaxes_view, 0)
    @test A_view == Aaxes_view
    @test A_view == Aaxes_view
end

@testset "push!" begin
    x = AxisArray([1], [:a])
    push!(x, :b => 2)
    @test keys(axes(x, 1)) == [:a, :b]
    @test x == [1, 2]
    pushfirst!(x, :pre_a => 0)
    @test keys(axes(x, 1)) == [:pre_a, :a, :b]
    @test x == [0, 1, 2]
end

@testset "cat_axis" begin
    @test @inferred(cat_axis(Axis(UnitMRange(1, 10)), SimpleAxis(UnitMRange(1, 10)), OneTo(20))) == 1:20
    @test @inferred(cat_axis(SimpleAxis(UnitMRange(1, 10)), SimpleAxis(UnitMRange(1, 10)), OneTo(20))) == 1:20
    @test @inferred(cat_axis(SimpleAxis(UnitMRange(1, 10)), Base.OneTo(10), OneTo(20))) == 1:20
    @test @inferred(cat_axis(SimpleAxis(Base.OneTo(10)), UnitMRange(1, 10), OneTo(20))) == 1:20
    @test @inferred(cat_axis(Axis(Base.OneTo(10)), UnitMRange(1, 10), OneTo(20))) == 1:20
end

#= FIXME cat tests
#
@testset "cat axes" begin
    @test @inferred(cat_axis(SimpleAxis(1:2), 2:4, 1:5)) === SimpleAxis(1:5)
    a, b = [1; 2; 3; 4; 5], [6 7; 8 9; 10 11; 12 13; 14 15];
    c, d = CartesianAxes((Axis(1:5),)), CartesianAxes((Axis(1:5), Axis(1:2)));
    #hcat_axes((Axis(1:4), Axis(1:2)), (Axis(1:4), Axis(1:2)))
    @test length.(@inferred(hcat_axes(c, d))) == length.(hcat_axes(a, b))
    @test length.(@inferred(hcat_axes(d, c))) == length.(hcat_axes(a, b))
    @test length.(@inferred(hcat_axes(CartesianAxes((10,)), CartesianAxes((10,))))) == (10, 2)
end
=#

@testset "hcat" begin
    a = AxisArray([1; 2; 3; 4; 5], (["a", "b", "c", "d", "e"],));
    b = [6 7; 8 9; 10 11; 12 13; 14 15];
    #@test keys.(@inferred(hcat_axes(a, b))) == (["a", "b", "c", "d", "e"], OneToMRange(3))

    @test keys.(axes(@inferred(hcat(a, b)))) == (["a", "b", "c", "d", "e"], OneToMRange(3))
    @test keys.(axes(@inferred(hcat(b, a)))) == (["a", "b", "c", "d", "e"], OneToMRange(3))

    @test keys.(axes(@inferred(hcat(a, a)))) == (["a", "b", "c", "d", "e"], 1:2)
    @test @inferred(hcat(a)) isa AbstractMatrix
    @test @inferred(hcat(hcat(a))) isa AbstractMatrix
end

@testset "vcat" begin
    a = AxisArray([1 2 3 4 5], (1:1, ["a", "b", "c", "d", "e"],));
    b = [6 7 8 9 10; 11 12 13 14 15];

    @test keys.(axes(vcat(a, b))) == keys.(axes(vcat(b, a)))
    @test keys.(axes(vcat(a, a))) == (1:2, ["a", "b", "c", "d", "e"])
    @test vcat(a) == a
end

@testset "cat" begin
    a = AxisArray(reshape(1:12, (3, 4)), (["a", "b", "c"], 2:5))
    b = parent(a)
    @test keys.(axes(cat(a, a, dims=3))) == (["a", "b", "c"], 2:5, Base.OneTo(2))
    @test keys.(axes(cat(b, a, dims=3))) == (["a", "b", "c"], 2:5, Base.OneTo(2))
    @test keys.(axes(cat(a, b, dims=3))) == (["a", "b", "c"], 2:5, Base.OneTo(2))

    @test keys.(axes(cat(a, a, a, dims=3))) == (["a", "b", "c"], 2:5, Base.OneTo(3))
    @test keys.(axes(cat(a, a, b, dims=3))) == (["a", "b", "c"], 2:5, Base.OneTo(3))
    @test keys.(axes(cat(a, b, a, dims=3))) == (["a", "b", "c"], 2:5, Base.OneTo(3))
    @test keys.(axes(cat(b, a, a, dims=3))) == (["a", "b", "c"], 2:5, Base.OneTo(3))

    # TODO this involves combining strings
    # @test keys(cat(a, a, dims=(1, 2))) == (['a','b','c', 'a','b','c'], [2,3,4,5, 2,3,4,5])
end

@testset "dropdims" begin
    axs = (Axis(["a", "b"]), Axis([:a]), Axis([1.0]), Axis(1:2))
    @test map(keys, AxisIndices.drop_axes(axs, 2)) == (["a", "b"], [1.0], 1:2)
    @test map(keys, AxisIndices.drop_axes(axs, (2, 3))) == (["a", "b"], 1:2)

    a = AxisArray(ones(10, 1, 1, 20), (2:11, [:a], 4:4, 5:24));

    @test dropdims(a; dims=2) == ones(10, 1, 20)
    @test keys.(axes(dropdims(a; dims=2))) == (2:11, 4:4, 5:24)

    @test dropdims(a; dims=(2, 3)) == ones(10, 20)
    @test keys.(axes(dropdims(a; dims=(2, 3)))) == (2:11, 5:24)
end

@testset "promote_shape" begin
    @test promote_shape(CartesianAxes((3, 4)), CartesianAxes((3,4,1,1,1))) ==
            (Base.OneTo(3), Base.OneTo(4), Base.OneTo(1), Base.OneTo(1), Base.OneTo(1))

    @test promote_shape(CartesianAxes((3,4,1,1,1)), CartesianAxes((3, 4))) ==
            (Base.OneTo(3), Base.OneTo(4), Base.OneTo(1), Base.OneTo(1), Base.OneTo(1))

    @test promote_shape(CartesianIndices((3,4,1,1,1)), CartesianAxes((3, 4))) ==
            (Base.OneTo(3), Base.OneTo(4), Base.OneTo(1), Base.OneTo(1), Base.OneTo(1))

    @test promote_shape(CartesianAxes((3,4,1,1,1)), CartesianIndices((3, 4))) ==
            (Base.OneTo(3), Base.OneTo(4), Base.OneTo(1), Base.OneTo(1), Base.OneTo(1))
end

@testset "+" begin
    a = AxisArray(ones(3), (2:4,))

    @testset "standard case" begin
        @test +(a) == ones(3)
        @test keys(axes(+(a), 1)) == 2:4

        @test +(a, a) == 2ones(3)
        @test keys(axes(+(a, a), 1)) == 2:4

        @test +(a, a, a) == 3ones(3)
        @test keys(axes(+(a, a, a), 1)) == 2:4
    end

    #= TODO decide what to do with axes that don't match
    @testset "Dimension disagreement" begin
        @test_throws DimensionMismatch +(
            NamedDimsArray{(:a, :b, :c, :d)}(zeros(3, 3, 3, 3)),
            NamedDimsArray{(:w, :x, :y, :z)}(ones(3, 3, 3, 3))
        )

        @test_throws DimensionMismatch +(
            NamedDimsArray{(:time,)}(zeros(3,)),
            NamedDimsArray{(:time, :value)}(ones(3, 3))
        )
    end
    =#

    @testset "Mixed array types" begin
        axs = (2:4, 3:5, 4:6, 5:7)
        lhs_sum = +(AxisArray(zeros(3, 3, 3, 3), axs), ones(3, 3, 3, 3));
        @test lhs_sum == ones(3, 3, 3, 3)
        @test keys.(axes(lhs_sum)) == axs

        rhs_sum = +(zeros(3, 3, 3, 3), AxisArray(ones(3, 3, 3, 3), axs))
        @test rhs_sum == ones(3, 3, 3, 3)
        @test keys.(axes(rhs_sum)) == axs

        #= TODO
        casts = (AxisArray{(:foo, :bar)}, identity)
        for (T1, T2, T3, T4) in Iterators.product(casts, casts, casts, casts)
            all(isequal(identity), (T1, T2, T3, T4)) && continue
            total = T1(ones(3, 6)) + T2(2ones(3, 6)) + T3(3ones(3, 6)) + T4(4ones(3, 6))
            @test total == 10ones(3, 6)
            @test dimnames(total) == (:foo, :bar)
        end
        =#
    end
end

@testset "-" begin
    # This is actually covered by the tests for + above, since that uses the same code
    # just one extra as a sensability check
    a = AxisArray(ones(3, 100), (2:4, 3:102))
    @test a - a == zeros(3, 100)
    @test keys.(axes(a - a)) == (2:4, 3:102)
end

@testset "scalar product" begin
    ax = (Axis(1:10), Axis(2:21), Axis(3:32), Axis(4:43), Axis(5:54))
    a = AxisArray(ones(10, 20, 30, 40, 50), ax);
    @test 10a == 10ones(10, 20, 30, 40, 50)
    @test keys.(axes(10a)) == keys.(ax)
end

@testset "reshape" begin
    V = AxisArray(Vector(1:8), [:a, :b, :c, :d, :e, :f, :g, :h]);

    A = @inferred(reshape(V, 4, 2));
    @test axes(A) == (Axis([:a, :b, :c, :d] => Base.OneTo(4)), SimpleAxis(Base.OneTo(2)))

    @test axes(@inferred(reshape(A, 2, :))) == (Axis([:a, :b] => Base.OneTo(2)), SimpleAxis(Base.OneTo(4)))
end

#= TODO

dims = 2
a = AxisArray(ones(10, 1, 1, 20), (2:11, 3:3, 4:4, 5:24))

p = dropdims(parent(a); dims=dims)
axs = AxisIndices.drop_axes(a, dims)

AxisIndices.unsafe_reconstruct(a, p, axs)
@testset "reshape" begin
    a = NamedDimsArray(rand(2, 3), (:r, :c))

    @test reshape(nda, 3, 2) isa Array
    @test reshape(nda, 1, :) isa Array
    @test reshape(nda, :) isa Array
    @test vec(nda) isa Array
end
=#

@testset "selectdim" begin
    a = AxisArray(reshape(1:6, (2, 3)), (2:3, 2:4));
    p = parent(a)

    @test selectdim(a, 1, 2) == a[2, :]
    @test keys.(axes(selectdim(a, 1, 2))) == (2:4,)

    @test vec(selectdim(a, 1, 2:2)) == a[2, :]
    @test keys.(axes(selectdim(a, 1, 2:2))) == (2:2, 2:4)
end

for f in (adjoint, transpose)
    @testset "$f" begin
        @testset "Vector $f" begin
            v = AxisArray([10, 20, 30], (2:4,))
            @test f(v) == [10 20 30]
            @test keys.(axes(f(v))) == (1:1, 2:4)

            # Make sure vector double adjoint gets you back to the start.
            @test f(f(v)) == [10, 20, 30]
            @test keys.(axes(f(f(v)))) == (2:4,)
        end

        @testset "Matrix $f" begin
            m = AxisArray([10 20 30; 11 22 33], (2:3, 2:4))
            @test f(m) == [10 11; 20 22; 30 33]
            @test keys.(axes(f(m))) == (2:4, 2:3)

            # Make sure implementation of matrix double adjoint is correct
            # since it is easy for the implementation of vector double adjoint broke it
            @test f(f(m)) == [10 20 30; 11 22 33]
            @test keys.(axes(f(f(m)))) == (2:3, 2:4)
        end
    end
end

# We test pinv here as it is defined in src/function_dims.jl
# using the same logic as permutedims, transpose etc
@testset "pinv" begin
    @testset "Matrix" begin
        a = AxisArray([1.0 2 3; 4 5 6], (2:3, 4:6))
        @test keys.(axes(pinv(a))) == (4:6, 2:3)
        @test a * pinv(a) ≈ [1.0 0; 0 1]
        @test keys.(axes(a * pinv(a))) == (2:3, 2:3)
    end

    @testset "Vector" begin
        v = AxisArray([1.0, 2, 3], (2:4,))
        @test keys.(axes(pinv(v))) == (1:1, 2:4)

        @test keys.(axes(pinv(pinv(v)))) == (2:4,)
        @test pinv(pinv(v)) ≈ v
    end
end

@testset "ReinterpretAxisArray" begin
    x = [1.0 2 3; 4 5 6]
    rx = reinterpret(Float32, x)
    ax = AxisArray(x, (2:3, 4:6))
    rax = reinterpret(Float32, ax)
    @test axes(rax) == (2:5, 4:6)
    @test rax == rx
end


text:
function tics(radec_min::T, radec_max::T, numx::T, ticsize::T, ra::Bool) where {T<:AbstractFloat}

    numtics = numx/ticsize
    if ra
        mul = T(4)
    else
        mul = T(60)
    end
    mins = abs(radec_min - radec_max)*mul
    rapix = (numx - 1)/mins
    incr = mins/numtics

    if incr >= 120
        incr = T(480)
    elseif incr >= 60
        incr = T(120)
    elseif incr >= 30
        incr = T(60)
    elseif incr >= 15
        incr = T(30)
    elseif incr >= 10
        incr = T(15)
    elseif incr >= 5
        incr = T(10)
    elseif incr >= 2
        incr = T(5)
    elseif incr >= 1
        incr = T(2)
    elseif incr >= 1//2
        incr = T(1)
    elseif incr >= 1//4
        incr = T(1//2)
    elseif incr >= 10//60
        incr = T(1//4)
    elseif incr >= 5//60
        incr = T(10//60)
    elseif incr >= 2//60
        incr = T(5//60)
    elseif incr >= 1//60
        incr = T(2//60)
    elseif incr >= 5//600
        incr = T(1//60)
    elseif incr >= 2//600
        incr = T(5//600)
    elseif incr >= 1//600
        incr = T(2//600)
    elseif incr >= 5//6000
        incr = T(1//600)
    elseif incr >= 2//6000
        incr = T(5//6000)
    elseif incr >= 1//6000
        incr = T(2//6000)
    else incr >= 0
        incr = T(1//6000)
    end

    ticsize = rapix * incr

    if ra && incr == 480
        incr = T(240)
    end
    if radec_min > radec_max
        incr = -incr
    end
    return ticsize, incr
end

"""
    tics(radec_min, radec_max, numx, ticsize[, ra=true]) -> ticsize, incr

### Purpose ###

Compute a nice increment between tic marks for astronomical images.

### Explanation ###

For use in labelling a displayed image with right ascension or declination
axes.  An approximate distance between tic marks is input, and a new value
is computed such that the distance between tic marks is in simple increments
of the tic label values.

### Arguements ###

* `radec_min` : minimum axis value (degrees).
* `radec_min` : maximum axis value (degrees).
* `numx` : number of pixels in x direction.
* `ticsize` : distance between tic marks (pixels).
* `ra` (optional boolean keyword): if true, incremental value would
  be in minutes of time. Default is false.

### Output ###

A 2-tuple `(ticsize, incr)`:
* `ticsize` : distance between tic marks (pixels).
* `incr` : incremental value for tic labels.  The format is dependent
  on the optional keyword. If true (i.e for right ascension), it's in minutes
  of time, else it's in minutes of arc (for declination).

### Example ###

```jldoctest
julia> using AstroLib

julia> tics(55, 60, 100.0, 1/2)
(0.66, 2.0)

julia> tics(30, 60, 12, 2, true)
(2.75, 30.0)
```

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
tics(radec_min::Real, radec_max::Real, numx::Real,
     ticsize::Real, ra::Bool=false) =
         tics(promote(float(radec_min), float(radec_max), float(numx),
                      float(ticsize))..., ra)

ext: jl


text:
using Pkg

index_template = """
using Interact, Mux

b = button("Click me")
ui = hbox(b, observe(b))

port = parse(Int64, ARGS[1])
wait(WebIO.webio_serve(page("/", req -> ui), port))
"""

procfile_template = """
web: julia --project index.jl \$PORT
"""

function create_index(path)
    open(joinpath(path, "index.jl"), "w") do f
        write(f, index_template)
    end
end

function create_procfile(path)
    open(joinpath(path, "Procfile"), "w") do f
        write(f, procfile_template)
    end
end

function create_project_toml(path)
    Pkg.activate(path)
    Pkg.add("Interact")
    Pkg.add("Mux")
end

function check_if_exists(path, cwd; file = true)
    if file
        return isfile(joinpath(cwd, path))
    end
    return isdir(joinpath(cwd, path))
end

function main()
    directory_name = ""
    app_name = "my-default-app-name"
    cwd = pwd()
    try
        directory_name = ARGS[1]
    catch
        println("No directory provided.")
        exit(1)
    end
    try
        app_name = ARGS[2]
    catch
        println("No appname provided.")
    end
    if !check_if_exists(directory_name, cwd, file = false)
        mkdir(joinpath(cwd, directory_name))
    end
    if !check_if_exists("$directory_name/index.jl", cwd)
        create_index(directory_name)
    end
    if !check_if_exists("$directory_name/Procfile", cwd)
        create_procfile(directory_name)
    end
    if !check_if_exists("$directory_name/Project.toml", cwd)
        create_project_toml(directory_name)
    end

    cd(joinpath(cwd, "$directory_name"))
    run(`git init`)
    run(`heroku create "$app_name" --buildpack https://github.com/Optomatica/heroku-buildpack-julia.git`)
    run(`git add -A`)
    run(`git commit -m "Initial commit"`)
    run(`git push heroku master`)
end

main()

max_stars_repo_name: abhishalya/interact-deploy

max_issues_repo_name: abhishalya/interact-deploy


text:
default_masktxt = "topo4mask.txt"
default_maskgrd = "topomask.grd"
tmp4grd = "tmp4grd.txt"

###################################################
"""
    filename = landmask_asc(topo::VisClaw.Topo, filename::AbstractString=default_masktxt)

output [x y z] data in txt for masking
"""
function landmask_asc(topo::VisClaw.Topo, filename::AbstractString=default_masktxt)
    xv = vec(repeat(topo.x, inner=(topo.nrows,1)))
    yv = vec(repeat(topo.y, outer=(topo.ncols,1)))
    topov = vec(topo.elevation);

    inds = topov .< 0.0 # ocean
    deleteat!(xv, inds)
    deleteat!(yv, inds)
    deleteat!(topov, inds)

    open(filename, "w") do file
        Base.print_array(file, [xv yv topov])
    end
    return filename
end
###################################################

###################################################
"""
    G = landmask_grd(txtfile::AbstractString=default_masktxt; grdfile::AbstractString="", kwargs...)

output masking grid
"""
function landmask_grd(txtfile::AbstractString=default_masktxt;
                      grdfile::AbstractString="", kwargs...)
    # check
    if !isfile(txtfile); error("Not found: $txtfile"); end

    # keyword args
    d = KWARG(kwargs)

    # (part of GMT.jl surface.jl)
    cmd = GMT.parse_common_opts(d, "", [:R :V_params :a :bi :di :e :f :h :i :r :yx])
    cmd = *(cmd...)
    #println(cmd)

    # (part of GMT.jl psmask.jl)
    cmd = GMT.parse_common_opts(d, cmd, [:I :UVXY :JZ :c :e :p :r :t :yx :params], true)
    cmd = *(cmd...)
    cmd = GMT.parse_these_opts(cmd, d, [[:C :end_clip_path], [:D :dump], [:F :oriented_polygons],
                    [:L :node_grid], [:N :invert], [:Q :cut_number], [:S :search_radius], [:T :tiles]])
    cmd = *(cmd...)
    #println(cmd)

    if isempty(grdfile)
        grdfile = default_maskgrd
        Gout = true
    else
        Gout = false
    end

    # grid
    GMT.gmt("grdmask \"$txtfile\" $cmd -G\"$grdfile\" ")

    # return
    if Gout
        G = GMT.gmt("read -Tg $grdfile")
        rm(grdfile, force=true)
        return G
    else
        return nothing
    end
end
###################################################

###################################################
"""
    G = tilegrd(tile::VisClaw.AMRGrid; length_unit::AbstractString="", kwargs...)

make a grid file of VisClaw.AMRGrid with landmask
"""
function tilegrd(tile::VisClaw.AMRGrid; length_unit::AbstractString="", kwargs...)
    # var
    var = VisClaw.keytile(tile)
    # prameters & options
    R = VisClaw.getR_tile(tile)
    Δ = tile.dx

    xvec, yvec, zdata = VisClaw.tilez(tile, var)
    xmat = repeat(xvec, inner=(length(yvec),1))
    ymat = repeat(yvec, outer=(length(xvec),1))

    if !isempty(length_unit)
        Δ = "$(Δ)"*length_unit
        # +ue?
        R = length_unit*R
    end

    # if NaN in all
    if !any(.!isnan.(zdata[:]))
        tmp_eta = "eta_tile.grd"
        faint="tmp.txt"
        open(faint, "w") do file
            Base.print_array(file, [xmat[:] ymat[:]])
        end
        GMT.gmt("grdmask $faint -R$R -I$Δ -S$Δ -NNaN/NaN/NaN -G$tmp_eta ")
        G = GMT.gmt("read -Tg $tmp_eta")
        rm(faint, force=true)
    else
        # eta grid
        G = GMT.surface([xmat[:] ymat[:] zdata[:]]; R=R, I=Δ, kwargs...)
    end

    # return value (GMT.GMTgrid)
    return G
end
###################################################
tilegrd(amrs::VisClaw.AMR, istep::Integer; length_unit::AbstractString="", kwargs...) =
tilegrd.(amrs.amr[istep]; length_unit=length_unit, kwargs...)
###################################################

###################################################
"""
    G = tilegrd_xyz(tile::VisClaw.AMRGrid; kwargs...)

make a grid file of VisClaw.AMRGrid
"""
function tilegrd_xyz(tile::VisClaw.AMRGrid; kwargs...)
    # var
    var = VisClaw.keytile(tile)
    # prameters & options
    R = VisClaw.getR_tile(tile)
    Δ = tile.dx

    xvec, yvec, zdata = VisClaw.tilez(tile, var)
    nx = length(xvec)
    ny = length(yvec)
    xvec = repeat(xvec, inner=(ny,1)) |> vec
    yvec = repeat(yvec, outer=(nx,1)) |> vec
    zvec = vec(zdata)
    inds = isnan.(zvec)
    deleteat!(xvec, inds)
    deleteat!(yvec, inds)
    deleteat!(zvec, inds)

    tmp4grd = "tmp4grd.txt"
    f = open(tmp4grd, "w"); Base.print_array(f, [xvec yvec zvec]); close(f)
    G = GMT.gmt("xyz2grd $(tmp4grd)  -R$(R) -I$(Δ)")
    rm(tmp4grd, force=true)

    # return value (GMT.GMTgrid)
    return G
end
###################################################
tilegrd_xyz(amrs::VisClaw.AMR, istep::Integer; kwargs...) = tilegrd_xyz.(amrs.amr[istep]; kwargs...)
###################################################

###################################################
"""
    G = tilegrd_mask(tile::VisClaw.AMRGrid, maskfile::AbstractString=""; length_unit::AbstractString="", kwargs...)

make a grid file of VisClaw.AMRGrid with landmask
"""
function tilegrd_mask(tile::VisClaw.AMRGrid, maskfile::AbstractString=""; length_unit::AbstractString="", kwargs...)
    # var
    var = VisClaw.keytile(tile)
    # prameters & options
    R = VisClaw.getR_tile(tile)
    Δ = tile.dx
    r = sqrt(2.0)Δ

    xvec, yvec, zdata = VisClaw.tilez(tile, var)
    xmat = repeat(xvec', inner=(length(yvec),1))
    ymat = repeat(yvec, outer=(length(xvec),1))

    tmp_mask = "mask_tile.grd"
    tmp_eta = "eta_tile.grd"
    eta_masked = "eta_masked.grd"

    if !isempty(length_unit)
        Δ = "$(Δ)"*length_unit
        r = "$(r)"*length_unit
        # +ue?
        R = length_unit*R
    end

    # makegrd
    # land mask grid
    VisClaw.landmask_grd(maskfile; grdfile=tmp_mask, R=R, I=Δ, S=r, N="0/0/NaN", kwargs...)

    # if NaN in all
    if !any(.!isnan.(zdata[:]))
        faint="tmp.txt"
        open(faint, "w") do file
            Base.print_array(file, [xmat[:] ymat[:]])
            #Base.print_array(file, [xmat[:] reverse(ymat, dims=1)[:]])
        end
        GMT.gmt("grdmask $faint -R$R -I$Δ -S$Δ -NNaN/NaN/NaN -G$tmp_eta ")
        rm(faint, force=true)
    else
        # eta grid
        GMT.surface([xmat[:] ymat[:] zdata[:]]; R=R, I=Δ, G=tmp_eta)
    end

    # masking
    GMT.gmt("grdmath $tmp_eta $tmp_mask OR = $eta_masked ")
    # read
    G = GMT.gmt("read -Tg $eta_masked ")

    rm(tmp_mask)
    rm(tmp_eta)
    rm(eta_masked)

    # return value (GMT.GMTgrid)
    return G
end
###################################################
tilegrd_mask(amrs::VisClaw.AMR, istep::Integer, maskfile::AbstractString=""; length_unit::AbstractString="", kwargs...) =
tilegrd_mask.(amrs.amr[istep], maskfile; length_unit=length_unit, kwargs...)


text:
using Gurobi
import MathProgBase
using Pkg

function mathprogbase_file(file::String)
    return joinpath(dirname(dirname(pathof(MathProgBase))), "test", file)
end

include(mathprogbase_file("linprog.jl"))
linprogtest(GurobiSolver(OutputFlag=0, InfUnbdInfo=1))

s = GurobiSolver()
MathProgBase.setparameters!(s, Silent=true, TimeLimit=100.0)

include(mathprogbase_file("linproginterface.jl"))
linprogsolvertest(s)

include(mathprogbase_file("mixintprog.jl"))
mixintprogtest(GurobiSolver(OutputFlag=0))

include(mathprogbase_file("quadprog.jl"))
quadprogtest(GurobiSolver(OutputFlag=0))
socptest(GurobiSolver(OutputFlag=0))
qpdualtest(GurobiSolver(OutputFlag=0, QCPDual=1))

include(mathprogbase_file("conicinterface.jl"))
coniclineartest(GurobiSolver(OutputFlag=0))
# The following tests are not passing on Gurobi 8.0.0 due to known bug
# on infeasibility check of SOC models.
# https://github.com/JuliaOpt/Gurobi.jl/pull/123
if Gurobi.version < v"8.0.0"
    conicSOCtest(GurobiSolver(OutputFlag=0))
    conicSOCRotatedtest(GurobiSolver(OutputFlag=0))
end
conicSOCINTtest(GurobiSolver(OutputFlag=0))


text:
######################
# Trait implementation
######################
#
# Provides @traitimpl macro which will do:
#  - check supertraits are implemented
#  - find module where the generic function is defined
#  - make methods
#
# Example syntax:
# @traitimpl Cmp{MyType1, MyType1} begin
#     isless(x::MyType1, y::MyType1) = x.t<y.t
# end
#
# ToDo:
# - implement parametric methods
# - more checks?

function get_fname(ex::Expr)
    # :(sin{T1<:Interger}(x::T1) = sin(x.t1))
    # -> :sin{T1<:Interger}
    out = ex.args[1].args[1]
    if typeof(out)==Symbol || out.head==:curly
        return out
    else
        error("Could not parse function definition: $ex")
    end
end

function get_fname_only(ex)
    # :(sin{T1<:Interger}(x::T1) = sin(x.t1))
    # -> :sin
    ex = get_fname(ex)
    isa(ex,Symbol) ? ex : ex.args[1]
end

function get_fsig(ex::Expr)
    ex.args[1].args[2:end]
end

function check_macro_body(bodyargs, implfs, trait)
    # TODO: this check is broken

    # check that there are no double defs
    length(trait().methods)==length(implfs) || error("Duplicate method definition(s)")
    # check right number of defs
    length(trait().methods)==length(implfs) || error("Not right number of method definitions")
    # check that the signature of implfs agrees with trait().methods
    for (f,sig) in trait().methods
        # for now just check length.
        if length(sig[1])!=length(get_fsig(implfs[f]))
            error("""Method definition:
                     $f  $sig
                     does not match implementation:
                     $(implfs[f])""")
        end
    end
    nothing
end

function parse_body(body::Expr)
    # TODO try and remove eval:
    implfs = Dict()
    for ln in Lines(body)
        if !isdefined(get_fname_only(ln))
            # define a standard generic function:
            eval_curmod(:($(get_fname_only(ln))() = error("Not defined")))
        end
        implfs[eval_curmod(get_fname_only(ln))] = ln
    end
    return implfs
end

function prefix_module!(ex::Expr, modname::Symbol)
    # Prefix the module of the function, so import is not necessary:
    # :(+(x::T1, y::T2) = x.t1 + y.t2) -> :(Base.+(x::T1, y::T2) = x.t1 + y.t2)
    if ex.head== :(=)
        if ex.args[1].head!=:call
            error("Not a function definition :\n$ex")
        end
    elseif ex.head!= :function
        error("Not a function definition:\n$ex")
    end

    fnname = get_fname(ex)
    fnname_only = get_fname_only(ex)
    if isa(fnname, Symbol)
        ex.args[1].args[1] = :($modname.$fnname)
    elseif fnname.head==:curly
        if isa(fnname.args[1], Symbol)
            ex.args[1].args[1].args[1] = :($modname.$fnname_only)
        end
    else
        error("something went wrong, $fnname is not a Symbol")
    end
    nothing
end

@doc """The `@traitimpl` macro can be used to implement a trait for a
     set of types.  Note however, that traits are also be implemented
     implicitly, i.e. any set of types is part of a trait if it
     fulfills it.

     Example continuing from the documentation of `@traitdef`, implementing the
     `MyArith` trait:

     ```
     type A; a end; type AB; b end
     @traitimpl MyArith{A,AB} begin
         +(x::A,y::AB) = A(x.a+y.b)
         -(x::A,y::AB) = A(x.a-y.b)
         *(x::A,y::AB) = A(x.a*y.b)
         /(x::A,y::AB) = A(x.a/y.b)
     end
     istrait(MyArith{A, AB}) # -> true
     ```

     Notes

     - the type annotations are mandatory.  No parameterized methods
       are allowed (for now).
     """ ->
macro traitimpl(head, body)
    ## Parse macro header
    name, paras, trait_expr = parsecurly(head)

    ## Parse macro body
    implfs = parse_body(body)
    #check_macro_body(body.args, implfs, trait) # doesn't work with associated types

    out = quote
        ## Check supertraits are implemented:
        if !istrait(traitgetsuper($trait_expr))
            istrait(traitgetsuper($trait_expr); verbose=true)
            throw(TraitException("""Not all supertraits of $($trait_expr) are implemented.
                                 Implement them first."""))
        end
    end
    ## Make methods
    for (fn, fndef) in implfs
        modname = module_name(Base.function_module(fn, Tuple{Vararg{Any}}))
        prefix_module!(fndef, modname)
        push!(out.args,fndef)
    end

    ## Assert that the implementation went smoothly
    push!(out.args, :(istrait($trait_expr) ? nothing :  @assert istrait($trait_expr, verbose=true)))
    return esc(out)
end


text:
export VarFunction

struct VarFunction{T} <: SingleVarMathmaticalFunction{T}
    var_index::UInt16

    VarFunction{T}(var_index) where {T} = new{T}(UInt16(var_index))

    (f::VarFunction{T})(call::MFCall) where {T} = call(call.args[f.var_index])

    Functions.string(f::VarFunction, emitter::FunctionEmitter) = "x_" * string(f.var_index)
    Base.isequal(f::VarFunction, f2::VarFunction) = f.value == f2.value
    Functions.emit(f::VarFunction{T}, emitter::FunctionEmitter{T}) where {T} = Expr(:ref, :ARGS, Int64(f.var_index))
end


text:
using PkgSRA, Test, Random
Random.seed!(13)

include("../utils/sindy_turing_utils.jl")

# Generate test data: Van der Pol (2-d)
include("../examples/example_vanDerPol.jl")
dat = Array(solve_vdp_system())
numerical_grad = numerical_derivative(dat, ts)

# Initialize main SRA object
this_model = sra_stateful_object(ts, tspan, dat, u0, numerical_grad)
prams = this_model.parameters
prams.sindyc_ensemble_parameters[:selection_criterion] =
    sindy_cross_validate;
# Several changes are required because there are only 2 variables here
prams.variable_names = ["x", "y"];
prams.sindy_terms_list = Iterators.product(1:3, 1:3)
prams.sindy_library["cross_terms"] = [2, 3] # Also include cubic terms

# Fit; iteration not necessary here
fit_first_model(this_model, 1);
calculate_subsampled_ind(this_model);
all_models = fit_model(this_model);

@testset "Exact Nonzero terms" begin
    names = get_nonzero_term_names(all_models[1]);
    @test names == ["dxdt_x", "dydt_x"]
    names = get_nonzero_term_names(all_models[2]);
    @test names == ["dxdt_x", "dydt_x", "dxdt_y"]
    names = get_nonzero_term_names(all_models[3]);
    @test names == ["dxdt_x", "dydt_x", "dxdt_y", "dxdt_xxx"]
    # TODO: finish
end;

@testset "Correct Regression" begin
    names = get_nonzero_term_names(this_model.sindy_model);
    @test names == ["dxdt_x",
                    "dydt_x",
                    "dxdt_y",
                    "dydt_y", # Spurious term
                    "dydt_", # Spurious term
                    "dxdt_xxx"]
end;


text:
######### CmdStan program example  ###########

using CmdStan

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  bernoullimodel = "
  data {
    int<lower=1> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    theta ~ beta(1,1);
    y ~ bernoulli(theta);
  }
  ";

  observeddata = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

  global stanmodel, rc, chn, nt, cnames, sdf

  tmpdir = ProjDir*"/tmp"

  stanmodel = Stanmodel(Sample(save_warmup=false, num_warmup=1000,
    num_samples=2000, thin=1), name="bernoulli", model=bernoullimodel,
    printsummary=false, tmpdir=tmpdir, output_format=:namedtuple);

  rc, nt, cnames = stan(stanmodel, observeddata, ProjDir, diagnostics=false,
    CmdStanDir=CMDSTAN_HOME);

  if rc == 0

    display(nt)

    # Ceate a ChainDataFrame
    sdf = read_summary(stanmodel)
    sdf[sdf.parameters .== :theta, [:mean, :ess]]
  end

end # cd


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    RSOCtoSOCBridge{T} <: Bridges.Variable.AbstractBridge

`RSOCtoSOCBridge` implements the following reformulation:

 * ``||x||_2^2 \\le 2tu`` into ``||v|| \\le w ``, with the substitution rules
   ``t = \\frac{w}{\\sqrt 2} + \\frac{v_1}{\\sqrt 2}``,
   ``u = \\frac{w}{\\sqrt 2} - \\frac{v_1}{\\sqrt 2}``, and
   ``x = (v_2,\\ldots,v_N)``.

## Source node

`RSOCtoSOCBridge` supports:

 * [`MOI.VectorOfVariables`](@ref) in [`MOI.RotatedSecondOrderCone`](@ref)

## Target node

`RSOCtoSOCBridge` creates:

 * [`MOI.VectorOfVariables`](@ref) in [`MOI.SecondOrderCone`](@ref)
"""
struct RSOCtoSOCBridge{T} <:
       SetMapBridge{T,MOI.SecondOrderCone,MOI.RotatedSecondOrderCone}
    variables::Vector{MOI.VariableIndex}
    constraint::MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.SecondOrderCone}
end

const RSOCtoSOC{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{RSOCtoSOCBridge{T},OT}

function MOI.Bridges.map_set(
    ::Type{<:RSOCtoSOCBridge},
    set::MOI.SecondOrderCone,
)
    return MOI.RotatedSecondOrderCone(MOI.dimension(set))
end

function MOI.Bridges.inverse_map_set(
    ::Type{<:RSOCtoSOCBridge},
    set::MOI.RotatedSecondOrderCone,
)
    return MOI.SecondOrderCone(MOI.dimension(set))
end

function MOI.Bridges.map_function(
    ::Type{<:RSOCtoSOCBridge{T}},
    func,
    i::MOI.Bridges.IndexInVector,
) where {T}
    scalars = MOI.Utilities.eachscalar(func)
    if i.value > 2
        return scalars[i.value]
    end
    t, u = scalars[1], scalars[2]
    ts = MOI.Utilities.operate!(/, T, t, sqrt(T(2)))
    us = MOI.Utilities.operate!(/, T, u, sqrt(T(2)))
    if i.value == 1
        return MOI.Utilities.operate!(+, T, ts, us)
    else
        return MOI.Utilities.operate!(-, T, ts, us)
    end
end

function MOI.Bridges.map_function(::Type{<:RSOCtoSOCBridge{T}}, func) where {T}
    scalars = MOI.Utilities.eachscalar(func)
    t, u, x = scalars[1], scalars[2], scalars[3:end]
    ts = MOI.Utilities.operate!(/, T, t, sqrt(T(2)))
    us = MOI.Utilities.operate!(/, T, u, sqrt(T(2)))
    return MOI.Utilities.operate(vcat, T, ts + us, ts - us, x)
end

# The map is an involution
function MOI.Bridges.inverse_map_function(BT::Type{<:RSOCtoSOCBridge}, func)
    return MOI.Bridges.map_function(BT, func)
end

# The map is symmetric
function MOI.Bridges.adjoint_map_function(BT::Type{<:RSOCtoSOCBridge}, func)
    return MOI.Bridges.map_function(BT, func)
end

# The map is a symmetric involution
function MOI.Bridges.inverse_adjoint_map_function(
    BT::Type{<:RSOCtoSOCBridge},
    func,
)
    return MOI.Bridges.map_function(BT, func)
end


text:
export ESDGSEM

struct ESDGSEM{L, C, G, A1, A2, A3, NFV, NFS}
  law::L
  cell::C
  grid::G
  MJ::A1
  MJI::A2
  faceMJ::A3
  volume_numericalflux::NFV
  surface_numericalflux::NFS

  function ESDGSEM(; law, cell, grid, volume_numericalflux, surface_numericalflux)
    M = mass(cell)
    _, J = components(metrics(grid))
    MJ = M * J
    MJI = 1 ./ MJ

    faceM = facemass(cell)
    _, faceJ = components(facemetrics(grid))

    faceMJ = faceM * faceJ

    args = (law, cell, grid, MJ, MJI, faceMJ, volume_numericalflux, surface_numericalflux)
    new{typeof.(args)...}(args...)
  end
end

function (dg::ESDGSEM)(dq, q, time)
  cell = dg.cell
  grid = dg.grid
  device = Bennu.device(arraytype(cell))
  dim = ndims(cell)
  Nq = size(cell)[1]
  @assert all(size(cell) .== Nq)
  @assert(length(eltype(q)) == numberofstates(dg.law))

  comp_stream = Event(device)

  kernel_type = :per_dir
  if kernel_type == :naive
    workgroup = ntuple(i -> i <= min(2, dim) ? Nq : 1, 3)
    ndrange = (length(grid) * workgroup[1], Base.tail(workgroup)...)
    comp_stream = esvolumeterm!(device, workgroup)(
      dg.law,
      dq,
      q,
      derivatives_1d(cell)[1],
      dg.volume_numericalflux,
      metrics(grid),
      dg.MJ,
      dg.MJI,
      points(grid),
      Val(dim),
      Val(Nq),
      Val(numberofstates(dg.law));
      ndrange,
      dependencies = comp_stream
    )
  elseif kernel_type == :per_dir
    workgroup = ntuple(i -> i <= dim ? Nq : 1, 3)
    ndrange = (length(grid) * workgroup[1], Base.tail(workgroup)...)
    for dir in 1:dim
      comp_stream = esvolumeterm_dir!(device, workgroup)(
        dg.law,
        dq,
        q,
        derivatives_1d(cell)[1],
        dg.volume_numericalflux,
        metrics(grid),
        dg.MJ,
        dg.MJI,
        points(grid),
        Val(dir),
        Val(dim),
        Val(Nq),
        Val(numberofstates(dg.law));
        ndrange,
        dependencies = comp_stream
      )
    end
  else
    error("Unknown kernel type $kernel_type")
  end

  Nfp = Nq ^ (dim  - 1)
  workgroup_face = (Nfp,)
  ndrange = (Nfp * length(grid),)
  faceix⁻, faceix⁺ = faceindices(grid)
  facenormal, _ = components(facemetrics(grid))
  comp_stream = surfaceterm!(device, workgroup_face)(
    dg.law,
    dq,
    q,
    Val(Bennu.connectivityoffsets(cell, Val(2))),
    dg.surface_numericalflux,
    dg.MJI,
    faceix⁻,
    faceix⁺,
    dg.faceMJ,
    facenormal,
    boundaryfaces(grid),
    points(grid),
    Val(dim);
    ndrange,
    dependencies = comp_stream
  )

  wait(comp_stream)
end

@kernel function esvolumeterm!(law,
                               dq,
                               q,
                               D,
                               volume_numericalflux,
                               metrics,
                               MJ,
                               MJI,
                               points,
                               ::Val{dim},
                               ::Val{Nq},
                               ::Val{Ns}) where {dim, Nq, Ns}
  @uniform begin
    FT = eltype(law)
    Nq1 = Nq
    Nq2 = dim > 1 ? Nq : 1
    Nq3 = dim > 2 ? Nq : 1

    q1 = MVector{Ns, FT}(undef)
    q2 = MVector{Ns, FT}(undef)
    x⃗1 = MVector{dim, FT}(undef)
    x⃗2 = MVector{dim, FT}(undef)
  end

  dqijk = @private FT (Ns,)

  pencil_q = @private FT (Ns, Nq3)
  pencil_x⃗ = @private FT (dim, Nq3)
  pencil_g3 = @private FT (3, Nq3)
  pencil_MJ = @private FT (Nq3,)

  l_q = @localmem FT (Nq1, Nq2, Ns)
  l_x⃗ = @localmem FT (Nq1, Nq2, dim)
  l_g = @localmem FT (Nq1, Nq2, min(2, dim), dim)

  e = @index(Group, Linear)
  i, j = @index(Local, NTuple)

  @inbounds begin
    @unroll for k in 1:Nq3
      ijk = i + Nq * (j - 1 + Nq * (k - 1))

      pencil_MJ[k] = MJ[ijk, e]
      @unroll for s in 1:Ns
        pencil_q[s, k] = q[ijk, e][s]
      end
      @unroll for d in 1:dim
        pencil_x⃗[d, k] = points[ijk, e][d]
      end
      if dim > 2
        @unroll for d in 1:dim
          pencil_g3[d, k] = metrics[ijk, e].g[3, d]
          pencil_g3[d, k] *= pencil_MJ[k]
        end
      end
    end

    @unroll for k in 1:Nq3
      @synchronize
      ijk = i + Nq * (j - 1 + Nq * (k - 1))

      @unroll for s in 1:Ns
        l_q[i, j, s] = pencil_q[s, k]
      end
      @unroll for d in 1:dim
        l_x⃗[i, j, d] = pencil_x⃗[d, k]
      end

      MJk = pencil_MJ[k]
      @unroll for d in 1:dim
        l_g[i, j, 1, d] = MJk * metrics[ijk, e].g[1, d]
        if dim > 1
          l_g[i, j, 2, d] = MJk * metrics[ijk, e].g[2, d]
        end
      end

      @synchronize

      fill!(dqijk, -zero(FT))

      @unroll for s in 1:Ns
        q1[s] = l_q[i, j, s]
      end
      @unroll for d in 1:dim
        x⃗1[d] = l_x⃗[i, j, d]
      end

      source!(law, dqijk, q1, x⃗1)

      MJIijk = 1 / pencil_MJ[k]
      @unroll for n in 1:Nq
        @unroll for s in 1:Ns
          q2[s] = l_q[n, j, s]
        end
        @unroll for d in 1:dim
          x⃗2[d] = l_x⃗[n, j, d]
        end

        f = twopointflux(volume_numericalflux, law, q1, x⃗1, q2, x⃗2)
        @unroll for s in 1:Ns
          Din = MJIijk * D[i, n]
          Dni = MJIijk * D[n, i]
          @unroll for d in 1:dim
            dqijk[s] -= Din * l_g[i, j, 1, d] * f[d, s]
            dqijk[s] += f[d, s] * l_g[n, j, 1, d] * Dni
          end
        end

        if dim > 1
          @unroll for s in 1:Ns
            q2[s] = l_q[i, n, s]
          end
          @unroll for d in 1:dim
            x⃗2[d] = l_x⃗[i, n, d]
          end
          f = twopointflux(volume_numericalflux, law, q1, x⃗1, q2, x⃗2)
          @unroll for s in 1:Ns
            Djn = MJIijk * D[j, n]
            Dnj = MJIijk * D[n, j]
            @unroll for d in 1:dim
              dqijk[s] -= Djn * l_g[i, j, 2, d] * f[d, s]
              dqijk[s] += f[d, s] * l_g[i, n, 2, d] * Dnj
            end
          end
        end

        if dim > 2
          @unroll for s in 1:Ns
            q2[s] = pencil_q[s, n]
          end
          @unroll for d in 1:dim
            x⃗2[d] = pencil_x⃗[d, n]
          end
          f = twopointflux(volume_numericalflux, law, q1, x⃗1, q2, x⃗2)
          @unroll for s in 1:Ns
            Dkn = MJIijk * D[k, n]
            Dnk = MJIijk * D[n, k]
            @unroll for d in 1:dim
              dqijk[s] -= Dkn * pencil_g3[d, k] * f[d, s]
              dqijk[s] += f[d, s] * pencil_g3[d, n] * Dnk
            end
          end
        end
      end

      ijk = i + Nq * (j - 1 + Nq * (k - 1))
      dq[ijk, e] += dqijk
    end
  end
end

@kernel function esvolumeterm_dir!(law,
                                   dq,
                                   q,
                                   D,
                                   volume_numericalflux,
                                   metrics,
                                   MJ,
                                   MJI,
                                   points,
                                   ::Val{dir},
                                   ::Val{dim},
                                   ::Val{Nq},
                                   ::Val{Ns}) where {dir, dim, Nq, Ns}
  @uniform begin
    FT = eltype(law)
    Nq1 = Nq
    Nq2 = dim > 1 ? Nq : 1
    Nq3 = dim > 2 ? Nq : 1

    q2 = MVector{Ns, FT}(undef)
    x⃗2 = MVector{dim, FT}(undef)
  end

  dqijk = @private FT (Ns,)

  q1 = @private FT (Ns,)
  x⃗1 = @private FT (dim,)

  l_g = @localmem FT (Nq ^ 3, 3)

  e = @index(Group, Linear)
  i, j, k = @index(Local, NTuple)

  @inbounds begin
    ijk = i + Nq * (j - 1 + Nq * (k - 1))

    MJijk = MJ[ijk, e]
    @unroll for d in 1:dim
      l_g[ijk, d] = MJijk * metrics[ijk, e].g[dir, d]
    end

    fill!(dqijk, -zero(FT))

    @unroll for s in 1:Ns
      q1[s] = q[ijk, e][s]
    end
    @unroll for d in 1:dim
      x⃗1[d] = points[ijk, e][d]
    end

    source!(law, dqijk, q1, x⃗1)

    @synchronize

    ijk = i + Nq * (j - 1 + Nq * (k - 1))

    MJIijk = MJI[ijk, e]
    @unroll for n in 1:Nq
      if dir == 1
        id = i
        ild = n + Nq * ((j - 1) + Nq * (k - 1))
      elseif dir == 2
        id = j
        ild = i + Nq * ((n - 1) + Nq * (k - 1))
      elseif dir == 3
        id = k
        ild = i + Nq * ((j - 1) + Nq * (n - 1))
      end

      @unroll for s in 1:Ns
        q2[s] = q[ild, e][s]
      end
      @unroll for d in 1:dim
        x⃗2[d] = points[ild, e][d]
      end

      f = twopointflux(volume_numericalflux, law, q1, x⃗1, q2, x⃗2)
      @unroll for s in 1:Ns
        Ddn = MJIijk * D[id, n]
        Dnd = MJIijk * D[n, id]
        @unroll for d in 1:dim
          dqijk[s] -= Ddn * l_g[ijk, d] * f[d, s]
          dqijk[s] += f[d, s] * l_g[ild, d] * Dnd
        end
      end
    end

    dq[ijk, e] += dqijk
  end
end


text:
#==============================================================================
    Convert variable name to string

    author:     longemen3000@stackoverflow
 =============================================================================#

"""
    @tostring

Return the variable, function, Symbol, object, et al name as a string.

@tostring(arg::Any)

# Examples
```jldoctest
julia> a = 55;
julia> bingo = "I'm the new sherrif";
julia> ch = 'd';
julia> hippy(x) = x^2
julia> var1 = @tostring a
"a"
julia> var2 = @tostring bingo
"bingo"
julia> var3 = @tostring(ch)     # macro as function call
"ch"
julia> var4 = @tostring hippy
"hippy"
# note: function names can be returned with string()
julia> string(hippy)
"hippy"
```
"""
macro tostring(aeg::Any)
    x = string(aeg)
    quote
        $x
    end
end

max_stars_repo_name: mahiki/DesertIslandDisk

max_issues_repo_name: mahiki/DesertIslandDisk

max_forks_repo_name: mahiki/DesertIslandDisk


text:
cd(@__DIR__); import Pkg; Pkg.activate(".")
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
Pkg.develop(path="$(homedir())/github/SimpleSDMLayers.jl")
using Revise
using SimpleSDMLayers
using Test

## Step 9 : Minor improvements ####

## 1. geotiff(file; coords)
file = "assets/wc2.1_10m_bio_1.tif"
methods(geotiff)
geotiff(SimpleSDMPredictor, file, 1; coords...)

geotiff(SimpleSDMPredictor, file, 1, coords)
geotiff(SimpleSDMPredictor, file, 1, coords)
# Whatever, too many methods to define and too hard to maintain

## 2. ==
l1, l2 = SimpleSDMPredictor(WorldClim, BioClim, 1:2)
l3 = copy(l1)
l4 = similar(l1)
replace!(l4, nothing => NaN)
l5 = SimpleSDMPredictor(replace(l1.grid, nothing => missing), l1)

@test l1 == l1
@test l1 === l1
@test l2 != l1
@test l3 == l1
@which l3 !== l1

@test l4 != l1
@test l4 != l4
@test !isequal(l4, l1)
@test isequal(l4, l4)

@test ismissing(l5 == l1)
@test ismissing(l5 == l5)
@test !isequal(l5, l1)
@test isequal(l5, l5)

# basepath = expanduser("~/julia/julia-1.6.1/share/julia/base/")
# @which isequal(l3, l1) # Base.operators.jl:123
# run(`code $(basepath)/operators.jl`)
# @which l3 == l1 # Base.Base.jl:87
# run(`code $(basepath)/Base.jl`)
# @which isequal(l3.grid, l1.grid) # Base.abstractarray.jl:1962
# @which l3.grid == l1.grid # Base.abstractarray.jl:1991
# run(`code $(basepath)/abstractarray.jl`)

# 3. replace!
l1 = SimpleSDMPredictor(WorldClim, BioClim, 1)
l2 = convert(SimpleSDMResponse, l2)
@which replace!(l1, nothing => NaN)
l1 isa SimpleSDMPredictor
l2 isa SimpleSDMPredictor
replace!(l1, nothing => NaN)
replace!(l2, nothing => NaN)
l2.grid

## 4. layer1[layer2]

l1, l2 = SimpleSDMPredictor(WorldClim, BioClim, 1:2)
l3 = SimpleSDMPredictor(WorldClim, BioClim, 1; left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
stride(l1) == stride(l3) # same stride

l4 = l1[l2] # works
l4 == l1 # but it's exactly the same

l5 = l1[l3] # doesn't work as layers have different sizes, but should work
l5 == l3 # should return true

l3[l1]

## 5. boundingbox(layer)

## 6. hash() method for custom type

l1 = SimpleSDMPredictor(WorldClim, BioClim, 1; left = 0.0, right = 10.0, bottom = 0.0, top = 10.0)
l2 = copy(l1)

l1 == l2
hash(l1) == hash(l2)
hash(l1.grid) == hash(l2.grid)

## 7. WorldClim no data
using ArchGDAL
import SimpleSDMLayers: _find_span
import SimpleSDMLayers: _prepare_layer_for_burnin
coords = (left = 0.0, right = 10.0, bottom = 0.0, top = 10.0)
layer = SimpleSDMPredictor(WorldClim, BioClim, 1; coords...)

wcfile = "assets/wc2.1_10m_bio_1.tif"
tempfile = tempname()
bandnumber = 1
left, right, bottom, top = fill(nothing, 4)

geotiff(tempfile, layer)
geotiff(SimpleSDMPredictor, tempfile)
geotiff(SimpleSDMPredictor, wcfile)

dataset = ArchGDAL.read(tempfile)
band = ArchGDAL.getband(dataset, bandnumber)
ArchGDAL.getnodatavalue(band)

dataset = ArchGDAL.read(wcfile)
band = ArchGDAL.getband(dataset, bandnumber)
nd = ArchGDAL.getnodatavalue(band)
nd2 = Float32(nd)

geotiff(tempfile, layer; nodata = Float32(nd))
geotiff(tempfile, layer; nodata = Float32(nd2))
geotiff(SimpleSDMPredictor, tempfile)

l = SimpleSDMPredictor(WorldClim, BioClim, 1)

f2 = tempname()
geotiff(f2, l; nodata=-3.4f38)
mp2 = geotiff(SimpleSDMResponse, f2)

@test typeof(mp2) <: SimpleSDMResponse
@test size(mp2) == size(l)
@test mp2 == l

## Issue comment
l = SimpleSDMPredictor(WorldClim, BioClim, 1)

# With -9999 as nodata
f1 = tempname()
geotiff(f1, l) # use -9999 as nodata by default
mp1 = geotiff(SimpleSDMResponse, f1)

length(mp1) == length(l) # false, but should be true
mp1 == l # false, same thing
replace(mp1, NaN => nothing) == l # true, but could be avoided

# With custom nodata
f2 = tempname()
geotiff(f2, l; nodata=-3.4f38) # this is the nodata from the WorldClim files
mp2 = geotiff(SimpleSDMPredictor, f2)

length(mp2) == length(l) # false
mp2 == l # false
replace(mp2, NaN => nothing) == l # true


text:
type Config
    runtime::Runtime
    agent::Agent
    queue::Queue
    is_alive::Bool
end

const DefaultQueueSize = 32

const debug_output = false

function debug_print(args...)
    if debug_output
        println(args...)
    end
end

global Defaults = Nullable{Config}()

export set_defaults
function set_defaults(rt::Runtime;
    agent = nothing,
    queue = nothing)
    global Defaults

    if agent == nothing
        # Prefer GPU Agents
        agents = HSA.all_agents(dev = HSA.DeviceTypeGpu)
        if isempty(agents)
            agents = HSA.all_agents()
        end

        agent = first(agents)
    end

    if queue == nothing
        max_q = HSA.agent_info_queue_max_size(agent)

        queue = Queue(agent, min(max_q, DefaultQueueSize), typ = HSA.QueueTypeSingle, error_callback = (s, q) -> begin
            str = HSA.status_string(s)
            println("QueueErrorHandler: ($s) -> $str")
        end)
    end

    cfg = Config(
            rt,agent,queue,true
        )

    Defaults = Nullable{Config}(cfg)

    finalizer(cfg, shutdown_managed)

    return cfg
end

function set_defaults()
    rt = Runtime()

    set_defaults(rt)
end

function get_defaults()
    return get(Defaults)
end

export get_or_init_defaults
function get_or_init_defaults()
    if isnull(Defaults) || !get_defaults().is_alive
        return set_defaults()
    else
        return get_defaults()
    end
end

export init_managed
init_managed() = get_or_init_defaults()

export clear_defaults
function clear_defaults()
    global Defaults
    shutdown_managed(Defaults)
    Defaults = Nullable{Config}()
end

function shutdown_managed(cfg)
    if isa(cfg, Nullable)
        if isnull(cfg)
            return
        else
            cfg = get(cfg)
        end
    elseif !isa(cfg, Config)
        error("cfg must be an HSA Config object to shut down")
    end

    if !cfg.is_alive
        return
    end

    # perform orderly shutdown of the referenced objects
    finalize(cfg.agent)
    finalize(cfg.queue)
    finalize(cfg.runtime)
    cfg.is_alive = false
end


text:
# Bertrand Nash example
# Generates random equilibrium data
using Distributions

indxID(ID) = (ID == :1) ? (1, 2) : (2, 1)

function bestResponsePrice(thetas, prices, error, pmax, ID)
	me, you = indxID(ID)
	b = thetas[you]*prices[you] + thetas[3]*error + thetas[4]
	p = (-b - sqrt(b^2 - 8 * thetas[me]))/ 4. / thetas[me]

	assert( p>= 0)
	return min(p, pmax)
end

#solve for equilibrium prices.  True model is log with a random, firm dependent effect
function solveNashPrices(thetas1::Array{Float64, 1}, thetas2::Array{Float64, 1},
						errors::Array{Float64, 2}, pmax; MAX_ITERS = 20, TOL=1e-6, trace=false,
		prices=fill(convert(Float64, pmax / 2), 2))
	for i = 1:MAX_ITERS
		prices_old = copy(prices)
		prices[1] = bestResponsePrice(thetas1, prices_old, errors[1], pmax, :1)
		prices[2] = bestResponsePrice(thetas2, prices_old, errors[2], pmax, :2)
		dist = norm(prices - prices_old)
		if trace
			println(prices', dist)
		end
		if dist < TOL
			return prices
		end
	end
	error("Max Iters reached")
end

function margRev(thetas, prices, error, ID)
	me, you = indxID(ID)
	1./prices[me] + 2 * thetas[me] * prices[me] + thetas[you] * prices[you] + thetas[3] * error + thetas[4]
end

max_stars_repo_name: nlperic/InverseVIs

max_issues_repo_name: nlperic/InverseVIs


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ContainerServiceDiagnosticsProfile <: SwaggerModel
    vmDiagnostics::Any # spec type: Union{ Nothing, ContainerServiceVMDiagnostics } # spec name: vmDiagnostics

    function ContainerServiceDiagnosticsProfile(;vmDiagnostics=nothing)
        o = new()
        validate_property(ContainerServiceDiagnosticsProfile, Symbol("vmDiagnostics"), vmDiagnostics)
        setfield!(o, Symbol("vmDiagnostics"), vmDiagnostics)
        o
    end
end # type ContainerServiceDiagnosticsProfile

const _property_map_ContainerServiceDiagnosticsProfile = Dict{Symbol,Symbol}(Symbol("vmDiagnostics")=>Symbol("vmDiagnostics"))
const _property_types_ContainerServiceDiagnosticsProfile = Dict{Symbol,String}(Symbol("vmDiagnostics")=>"ContainerServiceVMDiagnostics")
Base.propertynames(::Type{ ContainerServiceDiagnosticsProfile }) = collect(keys(_property_map_ContainerServiceDiagnosticsProfile))
Swagger.property_type(::Type{ ContainerServiceDiagnosticsProfile }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ContainerServiceDiagnosticsProfile[name]))}
Swagger.field_name(::Type{ ContainerServiceDiagnosticsProfile }, property_name::Symbol) =  _property_map_ContainerServiceDiagnosticsProfile[property_name]

function check_required(o::ContainerServiceDiagnosticsProfile)
    (getproperty(o, Symbol("vmDiagnostics")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ContainerServiceDiagnosticsProfile }, name::Symbol, val)
end


text:
export SimplePolyhedraLibrary, SimplePolyhedron

"""
    SimplePolyhedraLibrary{T}

Default library for polyhedra of dimension larger than 1 ([`IntervalLibrary`](@ref) is the default for polyhedra of dimension 1).
The library implements the bare minimum and uses the fallback implementation for all operations.
"""
struct SimplePolyhedraLibrary{T} <: PolyhedraLibrary
    solver::Union{Nothing, MPB.AbstractMathProgSolver}
    function SimplePolyhedraLibrary{T}(solver=nothing) where T
        new{T}(solver)
    end
end

similar_library(lib::SimplePolyhedraLibrary, d::FullDim, ::Type{T}) where T = default_library(d, T) # default_library allows to fallback to Interval if d is FullDim{1}

mutable struct SimplePolyhedron{N, T, HRepT<:HRepresentation{N, T}, VRepT<:VRepresentation{N, T}} <: Polyhedron{N, T}
    hrep::Union{HRepT, Nothing}
    vrep::Union{VRepT, Nothing}
    solver::Union{Nothing, MPB.AbstractMathProgSolver}
    function SimplePolyhedron{N, T, HRepT, VRepT}(hrep::Union{HRepT, Nothing}, vrep::Union{VRepT, Nothing}, solver::Union{Nothing, MPB.AbstractMathProgSolver}) where {N, T, HRepT<:HRepresentation{N, T}, VRepT<:VRepresentation{N, T}}
        new{N, T, HRepT, VRepT}(hrep, vrep, solver)
    end
end
function SimplePolyhedron{N, T, HRepT, VRepT}(hrep::HRepT, solver::Union{Nothing, MPB.AbstractMathProgSolver}) where {N, T, HRepT<:HRepresentation{N, T}, VRepT<:VRepresentation{N, T}}
    SimplePolyhedron{N, T, HRepT, VRepT}(hrep, nothing, solver)
end
function SimplePolyhedron{N, T, HRepT, VRepT}(vrep::VRepT, solver::Union{Nothing, MPB.AbstractMathProgSolver}) where {N, T, HRepT<:HRepresentation{N, T}, VRepT<:VRepresentation{N, T}}
    SimplePolyhedron{N, T, HRepT, VRepT}(nothing, vrep, solver)
end
function SimplePolyhedron{N, T, HRepT, VRepT}(hrep::HRepresentation{N}, solver::Union{Nothing, MPB.AbstractMathProgSolver}) where {N, T, HRepT<:HRepresentation{N, T}, VRepT<:VRepresentation{N, T}}
    SimplePolyhedron{N, T, HRepT, VRepT}(HRepT(hrep), solver)
end
function SimplePolyhedron{N, T, HRepT, VRepT}(vrep::VRepresentation{N}, solver::Union{Nothing, MPB.AbstractMathProgSolver}) where {N, T, HRepT<:HRepresentation{N, T}, VRepT<:VRepresentation{N, T}}
    SimplePolyhedron{N, T, HRepT, VRepT}(VRepT(vrep), solver)
end

library(::Union{SimplePolyhedron{N, T}, Type{<:SimplePolyhedron{N, T}}}) where {N, T} = SimplePolyhedraLibrary{T}()
default_solver(p::SimplePolyhedron) = p.solver
supportssolver(::Type{<:SimplePolyhedron}) = true

hvectortype(::Type{<:SimplePolyhedron{N, T, HRepT}}) where {N, T, HRepT} = hvectortype(HRepT)
vvectortype(::Type{SimplePolyhedron{N, T, HRepT, VRepT}}) where {N, T, HRepT, VRepT} = vvectortype(VRepT)

similar_type(::Type{<:SimplePolyhedron{M, S, HRepT, VRepT}}, d::FullDim{N}, ::Type{T}) where {M, S, HRepT, VRepT, N, T} = SimplePolyhedron{N, T, similar_type(HRepT, d, T), similar_type(VRepT, d, T)}

function SimplePolyhedron{N, T, HRepT, VRepT}(hits::HIt...; solver=nothing) where {N, T, HRepT, VRepT}
    SimplePolyhedron{N, T, HRepT, VRepT}(HRepT(hits...), solver)
end
function SimplePolyhedron{N, T, HRepT, VRepT}(vits::VIt...; solver=nothing) where {N, T, HRepT, VRepT}
    SimplePolyhedron{N, T, HRepT, VRepT}(VRepT(vits...), solver)
end

# Need fulltype in case the use does `intersect!` with another element
SimplePolyhedron{N, T}(rep::Representation{N}, solver::Union{Nothing, MPB.AbstractMathProgSolver}) where {N, T} = SimplePolyhedron{N, T}(MultivariatePolynomials.changecoefficienttype(rep, T), solver)
function SimplePolyhedron{N, T}(rep::HRepresentation{N, T}, solver::Union{Nothing, MPB.AbstractMathProgSolver}) where {N, T}
    HRepT = fulltype(typeof(rep))
    VRepT = dualtype(HRepT)
    SimplePolyhedron{N, T, HRepT, VRepT}(rep, solver)
end
function SimplePolyhedron{N, T}(rep::VRepresentation{N, T}, solver::Union{Nothing, MPB.AbstractMathProgSolver}) where {N, T}
    VRepT = fulltype(typeof(rep))
    HRepT = dualtype(VRepT)
    SimplePolyhedron{N, T, HRepT, VRepT}(rep, solver)
end

function polyhedron(rep::Representation{N}, lib::SimplePolyhedraLibrary{T}) where {N, T}
    SimplePolyhedron{N, polytypefor(T)}(rep, lib.solver)
end

function Base.copy(p::SimplePolyhedron{N, T}) where {N, T}
    if p.hrep !== nothing
        SimplePolyhedron{N, T}(p.hrep, p.solver)
    else
        SimplePolyhedron{N, T}(p.vrep, p.solver)
    end
end

hrepiscomputed(p::SimplePolyhedron) = p.hrep !== nothing
function computehrep!(p::SimplePolyhedron)
    # vrep(p) could trigger an infinite loop if both vrep and hrep are null
    p.hrep = doubledescription(p.vrep)
end
function hrep(p::SimplePolyhedron)
    if !hrepiscomputed(p)
        computehrep!(p)
    end
    return p.hrep
end
vrepiscomputed(p::SimplePolyhedron) = p.vrep !== nothing
function computevrep!(p::SimplePolyhedron)
    # hrep(p) could trigger an infinite loop if both vrep and hrep are null
    p.vrep = doubledescription(p.hrep)
end
function vrep(p::SimplePolyhedron)
    if !vrepiscomputed(p)
        computevrep!(p)
    end
    return p.vrep
end

function sethrep!(p::SimplePolyhedron, h::HRepresentation)
    p.hrep = h
end
function setvrep!(p::SimplePolyhedron, v::VRepresentation)
    p.vrep = v
end
function resethrep!(p::SimplePolyhedron, h::HRepresentation)
    p.hrep = h
    p.vrep = nothing
end
function resetvrep!(p::SimplePolyhedron, v::VRepresentation)
    p.vrep = v
    p.hrep = nothing
end


text:
using Laplace
using Documenter

makedocs(
    modules=[Laplace],
    authors="Michael Helton",
    sitename="Laplace.jl",
    format = Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ]
)

deploydocs(
    repo="github.com/heltonmc/Laplace.jl.git",
    devbranch = "main"
)


text:
using FFTW
using LinearAlgebra
using SparseArrays
using Random
using Distributions
using GaussianRandomFields

struct OfflineMatrices{T<:AbstractArray, S<:AbstractArray, U<:AbstractArray}

    rho_bar::S
    R12::T
    R21_bar::T
    R22::T # Covariance between observation stations
    R22_inv::T # Inverse of R22
    R12_invR22::T
    Lambda::S
    K::U
    L::T
    mu20::T

    buf1::T
    buf2::T

end

struct OnlineMatrices{T<:AbstractVector, S<:AbstractArray, U<:AbstractArray}

    z1_bar::T
    z2::T

    Z1::S
    Z2::S

    samples::U
    mean::U

end

struct Grid{T}
    nx::Int
    ny::Int
    dx::T
    dy::T
    x_length::T
    y_length::T
end

# Covariance function r(x,y), equation 1 in Dietrich & Newsam 1996
function covariance(x::T, y::T, covariance_structure::IsotropicCovarianceStructure{T}) where T
    return covariance_structure.σ^2 * apply(covariance_structure, [x, y])
end

# Extended covariance function ̄r(x,y), equation 8 of Dietrich & Newsam 1996
function extended_covariance(x::T, y::T, grid::Grid, covariance_structure::IsotropicCovarianceStructure{T}) where T

    if 0 <= x <= grid.x_length
        x_ext = x
    elseif grid.x_length <= x <= 2 * grid.x_length
        x_ext = 2 * grid.x_length - x
    else
        @error "value of x is out of bounds"
    end

    if 0 <= y <= grid.y_length
        y_ext = y
    elseif grid.y_length <= y <= 2 * grid.y_length
        y_ext = 2 * grid.y_length - y
    else
        @error "value of y is out of bounds"
    end

    return covariance(x_ext, y_ext, covariance_structure)

end

# Covariance between observations and the extended grid R̅₂₁, from equation 11 in Dietrich & Newsam 1996
function covariance_stations_extended_grid!(
    cov::AbstractMatrix{T},
    grid::Grid,
    grid_ext::Grid,
    stations::NamedTuple,
    covariance_structure::IsotropicCovarianceStructure{T}
) where T

    xid = 1:grid_ext.nx
    yid = 1:grid_ext.ny

    c = CartesianIndices((xid, yid))[:]

    @assert size(cov) == (stations.nst, length(c))

    for (i,ist,jst) in zip(1:stations.nst, stations.ist, stations.jst)
        cov[i,:] .= extended_covariance.(abs.(getindex.(c, 1) .- ist) .* grid_ext.dx,
                                         abs.(getindex.(c, 2) .- jst) .* grid_ext.dy,
                                         (grid,), (covariance_structure,))
    end

end

# Covariance between stations and the original grid R₂₁, from equations 3 and 4 in Dietrich & Newsam 1996
function covariance_grid_stations!(cov::AbstractMatrix{T}, grid::Grid, stations::NamedTuple, covariance_structure::IsotropicCovarianceStructure{T}) where T

    xid = 1:grid.nx
    yid = 1:grid.ny

    c = CartesianIndices((xid, yid))[:]

    @assert size(cov) == (length(c), stations.nst)

    for (i,ist,jst) in zip(1:stations.nst, stations.ist, stations.jst)
        cov[:,i] .= extended_covariance.(abs.(getindex.(c, 1) .- ist) .* grid.dx,
                                         abs.(getindex.(c, 2) .- jst) .* grid.dy,
                                         (grid,), (covariance_structure,))
    end

end

# Covariance between observations R̅₂₂ = R₂₂ + Σ, where R₂₂ is as defined in equations 3 and 4 in
# Dietrich & Newsam 96 and Σ is the observation noise covariance (see note at end of
# Section 2 in Dietrich & Newsam 1996)
function covariance_stations!(cov::AbstractMatrix{T}, grid::Grid, stations::NamedTuple, covariance_structure::IsotropicCovarianceStructure{T}, std::T) where T

    cov .= covariance.(abs.(stations.ist .- stations.ist') .* grid.dx,
                       abs.(stations.jst' .- stations.jst) .* grid.dy,
                       (covariance_structure,)) .+ I(stations.nst) .* std.^2

end

# First column vector rho_bar of covariance matrix among pairs of points of the extended grid R̅₁₁,
# from Dietrich & Newsam 1996 described in text between equations 11 and 12
function first_column_covariance_grid!(rho::AbstractVector{T}, grid::Grid, grid_ext::Grid, covariance_structure::IsotropicCovarianceStructure{T}) where T

    xid = 1:grid_ext.nx
    yid = 1:grid_ext.ny

    c = CartesianIndices((xid, yid))[:]

    rho .= extended_covariance.((getindex.(c, 1) .- 1) .* grid_ext.dx,
                                (getindex.(c, 2) .- 1) .* grid_ext.dy,
                                (grid,), (covariance_structure,))

end

# Normalized two-dimension discrete Fourier transofrm normalized by sqrt(̄n₁).
# Operates on the 2d data stored as a matrix.
# The argument `f` lets you switch between forward transform (`f=identity`) and
# inverse transform (`f=inv`).
# From Dietrich & Newsam 96 in text following equation 12
function normalized_2d_fft!(transformed_array::AbstractMatrix{<:Complex}, array::AbstractMatrix{S}, fft_plan::FFTW.FFTWPlan, fft_plan!::FFTW.FFTWPlan, grid_ext, f=identity) where {S}

    normalization_factor = f(sqrt(grid_ext.nx * grid_ext.ny))
    if pointer(transformed_array) == pointer(array)
        mul!(transformed_array, f(fft_plan!), array)
    else
        mul!(transformed_array, f(fft_plan), array)
    end
    transformed_array ./= normalization_factor

end

# Normalized two-dimension discrete Fourier transofrm normalized by sqrt(̄n₁).
# Operates on the 2d data stored as a vector.
# The argument `f` lets you switch between forward transform (`f=identity`) and
# inverse transform (`f=inv`).
# From Dietrich & Newsam 96 in text following equation 12
function normalized_2d_fft!(transformed_vector::AbstractVector{<:Complex}, vector::AbstractVector{S}, fft_plan::FFTW.FFTWPlan, fft_plan!::FFTW.FFTWPlan, grid_ext, f=identity) where {S}

    normalization_factor = f(sqrt(grid_ext.nx * grid_ext.ny))
    tmp_array = complex(reshape(vector, grid_ext.nx, grid_ext.ny))
    mul!(tmp_array, f(fft_plan!), tmp_array)
    transformed_vector .= @view(tmp_array[:]) ./ normalization_factor

end

# Decomposition of R̅₁₁, equation 12 of Deitrich & Newsam 1996
function WΛWH_decomposition!(transformed_array::AbstractMatrix{T},
                             array::AbstractMatrix{T},
                             offline_matrices::OfflineMatrices,
                             grid::Grid,
                             grid_ext::Grid,
                             fft_plan::FFTW.FFTWPlan,
                             fft_plan!::FFTW.FFTWPlan,
                             ) where T

    @assert size(array) == (grid.nx, grid.ny)

    extended_array = zeros(ComplexF64, grid_ext.nx, grid_ext.ny)

    extended_array[1:grid.nx, 1:grid.ny] .= array

    normalized_2d_fft!(extended_array, extended_array, fft_plan, fft_plan!, grid_ext)

    # Here we do an element-wise multiplication of the extended_array with the vector Lambda. This is identical to
    # Diagonal(Lambda) * extended_array[:], but avoids flattening and reshaping extended_array.
    normalized_2d_fft!(extended_array, reshape(offline_matrices.Lambda, grid_ext.nx, grid_ext.ny).*extended_array, fft_plan, fft_plan!, grid_ext, inv)

    transformed_array .= real.(@view(extended_array[1:grid.nx, 1:grid.ny]))

end

# Get a vector of values of field at positions of [stations.ist, stations.jst]
function get_values_at_stations(field::AbstractMatrix{T}, stations) where T

    values = zeros(T, length(stations.ist))

    for (num,(i,j)) in enumerate(zip(stations.ist, stations.jst))
        values[num] = field[i,j]
    end

    return values

end

# Allocate and compute matrices that do not depend on time-dependent variables (height and observations).
function init_offline_matrices(grid::Grid,
                               grid_ext::Grid,
                               stations::NamedTuple,
                               covariance_structure::IsotropicCovarianceStructure{T},
                               obs_noise_std::T,
                               fft_plan::FFTW.FFTWPlan,
                               fft_plan!::FFTW.FFTWPlan,
                               F::Type,
                               ) where T

    n1 = grid.nx * grid.ny # number of elements in original grid
    n1_bar = grid_ext.nx * grid_ext.ny # number of elements in extended grid

    C = complex(F)

    matrices = OfflineMatrices(Vector{F}(undef, n1_bar),                   #rho_bar
                               Matrix{F}(undef, n1, stations.nst),          #R12
                               Matrix{F}(undef, stations.nst, n1_bar),      #R21_bar
                               Matrix{F}(undef, stations.nst, stations.nst), #R22
                               Matrix{F}(undef, stations.nst, stations.nst), #R22_inv
                               Matrix{F}(undef, n1, stations.nst),          #R12_invR22
                               Vector{F}(undef, n1_bar),                   #Lambda (diagonal elements)
                               Matrix{C}(undef, stations.nst, n1_bar),      #K
                               Matrix{F}(undef, stations.nst, stations.nst), #L
                               Matrix{F}(undef, stations.nst, stations.nst), #mu20

                               Matrix{F}(undef, grid.nx, grid.ny), #buf1
                               Matrix{F}(undef, grid.nx, grid.ny)  #buf2
                               )

    first_column_covariance_grid!(matrices.rho_bar, grid, grid_ext, covariance_structure)
    covariance_grid_stations!(matrices.R12, grid, stations, covariance_structure)
    covariance_stations_extended_grid!(matrices.R21_bar, grid, grid_ext, stations, covariance_structure)
    covariance_stations!(matrices.R22, grid, stations, covariance_structure, obs_noise_std)
    matrices.R22_inv .= inv(matrices.R22)
    matrices.R12_invR22 .= matrices.R12 * matrices.R22_inv

    fourier_coeffs = Vector{C}(undef, n1_bar)
    normalized_2d_fft!(fourier_coeffs, matrices.rho_bar, fft_plan, fft_plan!, grid_ext, inv)
    matrices.Lambda .= sqrt(n1_bar) .* real.(fourier_coeffs)

    WHbar_R12 = Matrix{C}(undef, n1_bar, stations.nst)
    for i in 1:stations.nst
        normalized_2d_fft!(@view(WHbar_R12[:,i]), @view(matrices.R21_bar[i,:]), fft_plan, fft_plan!, grid_ext)
    end
    KH = Diagonal(matrices.Lambda)^(-1/2)*WHbar_R12
    matrices.K .= KH'

    A = real.(matrices.R22 .- matrices.K*KH)

    if A ≈ Hermitian(A)
        matrices.L .= cholesky(Hermitian(A)).L
    else
        error("R22 - K*KH is not hermitian! Increasing the grid size may help. See Dietrich and Newsam '96 for details.")
    end

    matrices.mu20 .= obs_noise_std^(-2) .* (matrices.R22 .- I(stations.nst) .* obs_noise_std^2)

    return matrices

end

# Allocate memory for matrices that will be updated during the time stepping loop.
function init_online_matrices(grid::Grid, grid_ext::Grid, stations::NamedTuple, nprt_per_rank::Int, T::Type)

    n1 = grid.nx * grid.ny # number of elements in original grid
    n1_bar = grid_ext.nx * grid_ext.ny # number of elements in extended grid

    C = complex(T)

    matrices = OnlineMatrices(Vector{C}(undef, n1_bar),
                              Vector{C}(undef, stations.nst),
                              Matrix{C}(undef, grid.nx, grid.ny),
                              Matrix{C}(undef, grid.nx, grid.ny),
                              Array{T}(undef, grid.nx, grid.ny, nprt_per_rank),
                              Array{T}(undef, grid.nx, grid.ny, nprt_per_rank)
                              )

    return matrices

end

# Calculate the mean for the optimal proposal of the height
function calculate_mean_height!(mean::AbstractArray{T,3},
                                height::AbstractArray{T,3},
                                offline_matrices::OfflineMatrices,
                                observations::AbstractVector{T},
                                stations::NamedTuple,
                                grid::Grid,
                                grid_ext::Grid,
                                fft_plan::FFTW.FFTWPlan,
                                fft_plan!::FFTW.FFTWPlan,
                                nprt_per_rank::Int,
                                obs_noise_std::T,
                                ) where T

    # The arguments for the WΛWH decompositions are matrices that only have nonzero values
    # at the indices of the stations. Store them as sparse arrays to save space.
    mu21 = sparse(stations.ist, stations.jst,
                  offline_matrices.R22_inv * (offline_matrices.mu20 * observations),
                  grid.nx, grid.ny)
    mu22 = obs_noise_std^(-2) * sparse(stations.ist, stations.jst, observations, grid.nx, grid.ny)

    # Compute WΛWH decompositions, results are dense matrices, store them in buffers
    # These correspond to mu21 and mu22 in Alex's code
    WΛWH_decomposition!(offline_matrices.buf1, mu21, offline_matrices, grid, grid_ext, fft_plan, fft_plan!)
    WΛWH_decomposition!(offline_matrices.buf2, mu22, offline_matrices, grid, grid_ext, fft_plan, fft_plan!)

    # Compute the difference of the decomposition results, store in offline_matrices.buf1
    # This corresponds to mu2 in Alex's code.
    # TODO: Check if WΛWH is linear and you could swap the order
    offline_matrices.buf2 .-= offline_matrices.buf1

    mu10 = Vector{T}(undef, stations.nst)

    # Loop over particles
    for iprt = 1:nprt_per_rank

        mul!(mu10, offline_matrices.R22_inv, get_values_at_stations(@view(height[:,:,iprt]), stations))
        mu10_sparse = sparse(stations.ist, stations.jst, mu10, grid.nx, grid.ny)

        # Compute decomposition of height values at stations times the inverse covariance matrix
        # The argument corresponds to mu10 and the outcome to mu11 in Alex's code
        WΛWH_decomposition!(offline_matrices.buf1, mu10_sparse, offline_matrices, grid, grid_ext, fft_plan, fft_plan!)

        # Compute the mean for the ith particle using mu2 and mu11
        # Skip storing the temporary mu1 in Alex's code
        mean[:,:,iprt] .= @view(height[:,:,iprt]) .- offline_matrices.buf1 .+ offline_matrices.buf2

    end

end

function sample_height_proposal!(height::AbstractArray{T,3},
                                 offline_matrices::OfflineMatrices,
                                 online_matrices::OnlineMatrices,
                                 observations::AbstractVector{T},
                                 stations::NamedTuple,
                                 grid::Grid,
                                 grid_ext::Grid,
                                 fft_plan::FFTW.FFTWPlan,
                                 fft_plan!::FFTW.FFTWPlan,
                                 nprt_per_rank::Int,
                                 rng::Random.AbstractRNG,
                                 obs_noise_std::T,
                                 ) where T

    @assert iseven(nprt_per_rank) "Number of particles per rank must be even to use the Optimal Filter"

    calculate_mean_height!(online_matrices.mean, height, offline_matrices, observations, stations, grid, grid_ext, fft_plan, fft_plan!, nprt_per_rank, obs_noise_std)

    i_n1 = LinearIndices((grid.nx, grid.ny))
    i_n1_bar = LinearIndices((grid_ext.nx, grid_ext.ny))

    e1 = Vector{ComplexF64}(undef, grid_ext.nx * grid_ext.ny)
    e2 = Vector{ComplexF64}(undef, stations.nst)

    for iprt in 1:2:nprt_per_rank

        @. e1 = complex(randn(rng), randn(rng))
        @. e2 = complex(randn(rng), randn(rng))

        # This gives the vector z1_bar
        normalized_2d_fft!(online_matrices.z1_bar, Diagonal(offline_matrices.Lambda)^(1/2) * e1, fft_plan, fft_plan!, grid_ext, inv)

        # This is the vector z2
        mul!(online_matrices.z2, offline_matrices.K, e1)
        mul!(online_matrices.z2, offline_matrices.L, e2, 1, 1)

        # Restrict z1_bar to Omega1 and reshape into an array
        online_matrices.Z1 .= online_matrices.z1_bar[i_n1_bar[1:grid.nx,1:grid.nx]]
        # Multiply z2 with R12*R22^-1 and reshape the result into an array
        online_matrices.Z2 .= (offline_matrices.R12_invR22 * online_matrices.z2)[i_n1]

        online_matrices.samples[:,:,iprt    ] .= @view(online_matrices.mean[:,:,iprt    ]) .+ real.(online_matrices.Z1 .- online_matrices.Z2)
        online_matrices.samples[:,:,iprt + 1] .= @view(online_matrices.mean[:,:,iprt + 1]) .+ imag.(online_matrices.Z1 .- online_matrices.Z2)

    end

end

function get_log_weights!(log_weights::AbstractVector{T},
                          obs::AbstractVector{T},
                          obs_model::AbstractMatrix{T},
                          matrices::OfflineMatrices) where T

    #TODO: Is this just the same as get_log_weights! in ParticleDA with R22_inv as cov_obs?
    # DONE: No, ParticleDA.get_log_weights! fails with
    # ERROR: PosDefException: matrix is not Hermitian; Cholesky factorization failed.

    nprt = size(obs_model,2)

    for iprt in 1:nprt

        log_weights[iprt] = -0.5 * (obs - obs_model[:,iprt])' * matrices.R22_inv * (obs - obs_model[:,iprt])

    end

    #Normalization is done in a separate function due to parallelism optimisation

end


text:
# This file was generated, do not modify it. # hide
for editcol in 1:k
    for refcol in (editcol+1):k
        m = -y[refcol, editcol]
        for i in refcol:n
            y[i,editcol] += m * y[i,refcol]
        end
    end
end

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
## Closeness and harmonic centrality ##

# TODO: verify this against a correct implementation of Dijkstra's algorithm
# TODO: set up array reuse for `closeness` and `harmonic_centrality`

function distances(g::Graph, n::Integer,
                   dists=fill(typemax(Int), nodecount(g)),
                   visited=falses(nodecount(g)),
                   stack=zeros(Int, nodecount(g)))
    head, tail = 1, 2
    dists[n] = 0
    stack[head] = n
    visited[n] = true
    @inbounds while head < tail
        i = stack[head]
        @inbounds for j = neighbors(g, i)
            if !visited[j]
                if dists[j] > dists[i] + 1
                    dists[j] = dists[i] + 1
                end
                stack[tail] = Int(j)
                tail += 1
                visited[j] = true
            end
        end
        head += 1
    end
    @inbounds for i = 1:tail - 1
        visited[stack[i]] = false
    end
    return dists
end

"""
    closeness(g::Graph[, v::Integer])

Closeness centrality of node `v`.

# References
[Tore Opsahl's blog](https://toreopsahl.com/2010/03/20/closeness-centrality-in-networks-with-disconnected-components/)
"""
function closeness(g::Graph, v::Integer)
    n, s = 0, 0
    for d = distances(g, v)
        if d < typemax(Int)
            n += 1
            s += d
        end
    end
    return n / s
end

closeness(g::Graph) = [closeness(g, i) for i = nodes(g)]

"""
    harmonic_centrality(g::Graph[, v::Integer])

Harmonic centrality for node `v`, a better-behaved version of closeness centrality

# References
Rochat, Yannick. 2009. “Closeness Centrality Extended to Unconnected Graphs:
The Harmonic Centrality Index.” In ASNA.

[Tore Opsahl's blog](https://toreopsahl.com/2010/03/20/closeness-centrality-in-networks-with-disconnected-components/)
"""
function harmonic_centrality(g::Graph, v::Integer)
    s = 0.
    for d = distances(g, v)
        if 0 < d < typemax(Int)
            s += 1 / d
        end
    end
    return s
end

harmonic_centrality(g::Graph) = [harmonic_centrality(g, i) for i = nodes(g)]