text:
using Base.Test, Base.Random

using Main.BayesianNonparametricStatistics

srand(123)

sqrttwo = sqrt(2.0)

@testset "BayesianNonparametricStatistics package" begin
    @testset "model.jl" begin
        @test supertype(SDEModel) == AbstractModel
        @test supertype(AbstractModel) == Any
        for σ in -2.0:0.1:2.0
            for beginvalue in -2.0:0.1:2.0
                for endtime in 0.1:0.1:2.0
                    for Δ in 0.01:0.01:0.99
                        M = SDEModel(σ, beginvalue, endtime, Δ)
                        @test M.σ == σ
                        @test M.beginvalue == beginvalue
                        @test M.endtime == endtime
                        @test M.Δ==Δ
                    end
                end
            end
        end
        σ = identity
        M = SDEModel(σ, 0.0,10.0, 0.01)
        @test M.σ == identity
        σ2(x) = 1.0
        M = SDEModel(σ2, 0.0,10.0,0.01)
        @test M.σ == σ2

        σ3(x) = 2+sin(x)
        M = SDEModel(σ3, 1.3, 101.0, 0.5)
        @test M.σ == σ3
        @test M.beginvalue == 1.3
        @test M.endtime == 101.0
        @test M.Δ == 0.5

        @test_throws ArgumentError SDEModel(sin, 0.0, -1.0, 0.1)
        @test_throws ArgumentError SDEModel(sin, 0.0, 1.0, -0.1)
        @test_throws ArgumentError SDEModel(sin, 0.0, -1.0, -0.1)
        @test_throws ArgumentError SDEModel(1.0, 0.0, 0.0, 0.1)
        @test_throws ArgumentError SDEModel(1.0, 0.0, 10.0, 0.0)
    end

    #Test functions are correct, complete and the same in both versions.
    @testset "samplepath.jl" begin
        @test supertype(SamplePath) == AbstractSamplePath

        @test supertype(AbstractSamplePath) == Any

        X = SamplePath(0.0:0.1:1.0, sin)
        @test X.timeinterval == 0.0:0.1:1.0
        @test X.samplevalues == sin.(0.0:0.1:1.0)

        X = SamplePath(0.0:0.1:1.0, 1.0:11.0)

        @test X.timeinterval == 0.0:0.1:1.0
        @test X.samplevalues == 1.0:11.0

        X = SamplePath(0.0:0.1:1.0, x->2*x)
        @test X.timeinterval == 0.0:0.1:1.0
        @test X.samplevalues == collect(0.0:0.2:2.0)

        @test step(X) == 0.1
        @test length(X) == 11

        X = SamplePath([1., 2.], [1.1, 2.1])
        @test X.timeinterval == [1., 2.]
        @test X.samplevalues == [1.1, 2.1]

        @test length(X) == 2

        X = SamplePath([1.0, 3.0, 5.0, 9.0], x->x^2)
        @test X.timeinterval == [1.0, 3.0, 5.0, 9.0]
        @test X.samplevalues == [1.0, 9.0, 25.0, 81.0]

        X = SamplePath(collect(-2.0:0.01:2.0), sin)
        X.timeinterval == collect(-2.0:0.01:2.0)
        X.samplevalues == sin.(collect(-2.0:0.01:2.0))

        X = SamplePath(-2.0:0.01:2.0, sin)
        X.timeinterval == collect(-2.0:0.01:2.0)
        X.samplevalues == sin.(collect(-2.0:0.01:2.0))

        X = SamplePath(0.0:0.1:2π, sin)
        @test length(X) == length(0.0:0.1:2π)

        @test_throws DimensionMismatch SamplePath(0.0:2.0, [1.0, 2.0])

        @test_throws ArgumentError SamplePath(0.0:-0.1:-0.1, [0.0, 2.0])

        @test_throws DimensionMismatch SamplePath([0.1, 2.0],[3.0])

        @test_throws ArgumentError SamplePath([0.1, 2.0, 1.0], [1.0,2.0,3.0])

        @test BayesianNonparametricStatistics.isincreasing(0:0.1:1)
        @test !BayesianNonparametricStatistics.isincreasing(0:-0.1:-1.0)
        @test !BayesianNonparametricStatistics.isincreasing([0.0, 0.1, 0.2, 0.1, 0.3, 0.4])
        @test !BayesianNonparametricStatistics.isincreasing([0.0, -0.1, 0.0, 0.2, 0.3])
        @test !BayesianNonparametricStatistics.isincreasing([1, 2, 3, 4, -10, 11, 12])
        @test BayesianNonparametricStatistics.isincreasing([1, 2, 3, 4, 500, 600, 2000, 10000])
        @test BayesianNonparametricStatistics.isincreasing([1., 2., 2.00000001, 2.00000002])
        @test !BayesianNonparametricStatistics.isincreasing([2.,2.,2.,2.,2.])
        @test !BayesianNonparametricStatistics.isincreasing([2.0,2.0,2.1,2.2,2.3])
        @test !BayesianNonparametricStatistics.isincreasing([2.0,2.1,2.2,2.2,2.3])
        @test !BayesianNonparametricStatistics.isincreasing([2.1,2.2,2.3,2.4,2.4])
        @test !BayesianNonparametricStatistics.isincreasing([2.,2., 3.0, 3.0001])
    end

    #Test functions are correct, complete and the same in both versions.
    @testset "basisfunctions.jl" begin
        @test_throws AssertionError fourier(-1)
        @test_throws AssertionError fourier(-10)

        numberofbasisfunctionstested = 100
        numberoffaberschauderlevelstested = convert(Int64, round(log2(numberofbasisfunctionstested), RoundUp)) - 1
        numberoffourierfunctionstested = numberofbasisfunctionstested
        #The Fourier functions are orthogonal.
        t = 0.0:1/(4*ceil(numberoffourierfunctionstested/2)):1.0
        tprecise = 0.0:10.0^-1/(4*ceil(numberoffourierfunctionstested/2)):1.0
        for k in 0:numberoffourierfunctionstested-1
            for l in k+1:numberoffourierfunctionstested
                phikvalues = fourier(k).(tprecise[1:end-1])
                philvalues = fourier(l).(tprecise[1:end-1])
                @test abs(sum(step(tprecise)*phikvalues.*philvalues)) < 10.0^-10
            end
        end
        # And the square integrates to one.
        for k in 0:numberoffourierfunctionstested
            values = fourier(k).(tprecise)
            @test abs(step(tprecise)*sum(values.*values)-1.0) < 10.0^-2
        end

        # They are one-periodic.
        for k in 0:numberoffourierfunctionstested
            for y in t
                value = fourier(k)(y)
                for m in -10.0:10.0
                    valuetwo = fourier(k)(y+m)
                    @test abs(value - valuetwo) < 10.0^-10
                end
            end
        end

        # The function fourier(0) is constant 1.0
        for x in -10.0:0.01:10.0
            @test fourier(0)(x) == 1.0
        end
        #
        # For k≥1, the functions have minimum -sqrt(2) and maximum sqrt(2).
        for k in 1:numberoffourierfunctionstested
            if k % 2 == 1
                t = 0.0:0.1*1/(2k+2):1.0
            else
                t = 0.0:0.1*1/k:1.0
            end
            values = fourier(k).(t)
            @test abs(minimum(values)+sqrttwo) < 10.0^-5
            @test abs(maximum(values)-sqrttwo) < 10.0^-5
        end
        # # For k≥1:
        # # Functions \phi_k with k odd are zero when x=m/(k+1), m=0,...,k.
        # # Functions \phi_k with k even are zero when x=(1+2m)/(2k), m=0,...,k-1.
        for k in 1:2:numberoffourierfunctionstested
            for m in 0:k
                @test abs(fourier(k)(m/(k+1))) < 10.0^-5
            end
        end
        #
        for k in 2:2:numberoffourierfunctionstested
            for m in 0:k-1
                @test abs(fourier(k)((1+2m)/(2k))) < 10.0^-5
            end
        end
        #
        # # When k is odd, then \phi_k(x)=\sqrt(2) when x = (4m+1)/(2k+2), m=0,..,(k-1)/2.
        for k in 1:2:numberoffourierfunctionstested
            for m in 0:convert(Int64,(k-1)/2)
                @test abs(fourier(k)((4m+1)/(2k+2))-sqrt(2)) < 10.0^-5
            end
        end
        #
        # When k is odd, then \phi_k(x)= -sqrt(2) when x = (3+4m)/(2k+2), m=0,..,(k-1)/2.
        for k in 1:2:numberoffourierfunctionstested
            for m in 0:convert(Int64, (k-1)/2)
                @test abs(fourier(k)((3+4m)/(2k+2))+sqrt(2)) < 10.0^-5
            end
        end
        #
        # When k is even, then \phi_k(x)=sqrt(2) when x = 2m/k, m=0,...,k/2-1.
        for k in 2:2:numberoffourierfunctionstested
            for m in 0:convert(Int64, k/2)-1
                @test abs(fourier(k)(2m/k)-sqrt(2)) < 10.0^-5
            end
        end
        #
        # when k is even, then \phi_k(x)=-sqrt(2) when x=(2m+1)/k, m = 0,...,k/2-1.
        for k in 2:2:numberoffourierfunctionstested
            for m in 0:convert(Int64, k/2)-1
                @test abs(fourier(k)((2m+1)/k)+sqrt(2)) < 10.0^-5
            end
        end
        #
        @test_throws AssertionError faberschauder(-1, 1)
        @test_throws AssertionError faberschauder(-3, 1)
        @test_throws AssertionError faberschauder(1, 0)
        @test_throws AssertionError faberschauder(1, -2)
        @test_throws AssertionError faberschauder(1, 3)

        # # Below we test Faber-Schauder functions on their mathematical properities, up
        # # to level numberoffaberschauderlevelstested, defined below.
        t = 0.:2.0^(-numberoffaberschauderlevelstested-2):1.0

        # Faber-Schauder functions are level-wise orthogonal, for levels j≥1.
        for j in 1:numberoffaberschauderlevelstested
            for k in 1:2^j-1
                for l in k+1:2^j
                    psijkvalues = faberschauder(j,k).(t[1:end-1])
                    psijlvalues = faberschauder(j,l).(t[1:end-1])
                    @test abs(step(t)*sum(psijkvalues.*psijlvalues))<10.0^-3
                end
            end
        end
        #
        # They integrate to 2^(-j-1)
        psionevalues = faberschauderone.(t[1:end-1])
        @test abs(step(t)*sum(psionevalues)-2.0^-1) < 10.0^-3
        for j in 1:numberoffaberschauderlevelstested
            for k in 1:2^j
                psijkvalues = faberschauder(j,k).(t[1:end-1])
                @test abs(step(t)*sum(psijkvalues)-2.0^(-j-1)) < 10.0^-3
            end
        end

        # psi_1 has maximum 1 at 0 and 1.
        @test abs(faberschauderone(0.0)-1.0)<10.0^-5
        @test abs(faberschauderone(1.0)-1.0)<10.0^-5

        # psi_1 has minimum 0 at 0.5.
        @test abs(faberschauderone(0.5)) < 10.0^-5

        # psi_{j,k} has maximum 1.0 at 2.0^(-j-1)+(k-1)*2.0^-j
        for j in 0:numberoffaberschauderlevelstested
            for k in 1:2^j
                @test abs(faberschauder(j,k)(2.0^(-j-1)+(k-1)*2.0^-j)-1.0)<10.0^-5
            end
        end

        #psi_{j,k} has minimum 0.0 at (k-1)*2.0^-j and k*2.0^-j
        for j in 0:numberoffaberschauderlevelstested
            for k in 1:2^j
                @test abs(faberschauder(j,k)((k-1)*2.0^-j))<10.0^-5
                @test abs(faberschauder(j,k)(k*2.0^-j))<10.0^-5
            end
        end

        # They are one-periodic.
        x = 0.0:0.01:0.99
        for y in x
            value = faberschauderone(y)
            for k in -10.0:10.0
                valuetwo = faberschauderone(y+k)
                @test abs(value-valuetwo) < 10.0^-5
            end
        end

        for j in 0:numberoffaberschauderlevelstested
            for k in 1:2^j
                for y in x
                    value = faberschauder(j,k)(y)
                    for m in -10.:10.
                        valuetwo = faberschauder(j,k)(y+m)
                        @test abs(value-valuetwo) < 10.0^-5
                    end
                end
            end
        end
        # End of tests for mathematical properities of Faber-Schauder functions.
    end

    #Test functions are correct, complete and the same in both versions.
    @testset "SDE.jl" begin
        @test supertype(SDE) == AbstractSDE

        @test supertype(AbstractSDE) == Any

        # The following should represent a Brownian motion.
        model = SDEModel(1.0, 0.0, 1.0, 0.001)
        sde = SDE(x->0.0, model)

        lengthvector = 1000
        x = Vector{Float64}(lengthvector)

        for k in 1:lengthvector
            X = rand(sde)
            x[k] = X.samplevalues[end]
        end

        @test abs(mean(x)) < 0.1
        @test 0.9 < var(x) < 1.1

        X = rand(sde)

        @test length(X) == length(0.0:sde.model.Δ:sde.model.endtime)
        @test X.timeinterval == 0.0:sde.model.Δ:sde.model.endtime

        model = SDEModel(identity, 1.0, 1.0, 0.001)
        geometricbrownianmotion = SDE(x->0.5*x, model)

        for k in 1:lengthvector
            X = rand(geometricbrownianmotion)
            x[k] = X.samplevalues[end]
        end

        @test abs(e^(0.5)-mean(x)) < 0.1

        X = rand(geometricbrownianmotion)

        @test length(X) == length(0.0:geometricbrownianmotion.model.Δ:geometricbrownianmotion.model.endtime)
        @test X.timeinterval ==  0.0:geometricbrownianmotion.model.Δ:geometricbrownianmotion.model.endtime

        sde = SDE(x->0.0, SDEModel(1.0, 0.0, 1.0, 0.001))

        for k in 1:lengthvector
            X = rand(sde)
            x[k] = X.samplevalues[end]
        end

        @test abs(mean(x)) < 0.1
        @test 0.9 < var(x) < 1.1

        X = rand(sde)
        @test X.timeinterval == 0.0:sde.model.Δ:sde.model.endtime
        @test length(X) == length(0.0:sde.model.Δ:sde.model.endtime)

        ornsteinuhlenbeckprocess = SDE(x-> -x, SDEModel(1.0, 2.72, 1.0, 0.001))

        for k in 1:lengthvector
            X = rand(ornsteinuhlenbeckprocess)
            x[k] = X.samplevalues[end]
        end

        @test abs(mean(x) - 1) < 0.1

        X = rand(ornsteinuhlenbeckprocess)

        @test X.timeinterval == 0.0:ornsteinuhlenbeckprocess.model.Δ:ornsteinuhlenbeckprocess.model.endtime
        @test length(X) == length(0.0:ornsteinuhlenbeckprocess.model.Δ:ornsteinuhlenbeckprocess.model.endtime)

        model = SDEModel(1.0, 0.0, 10.0, 1.0)
        prevXval = 1.0
        @test BayesianNonparametricStatistics.calculatenextsamplevalue(prevXval, identity, model, 1.0) == 3.0

        model = SDEModel(identity, 0.0, 10.0, 1.0)
        prevXval = 1.0
        @test BayesianNonparametricStatistics.calculatenextsamplevalue(prevXval, identity, model, 1.0) == 3.0
        # End tests "sde.jl".
    end

    #Test functions are correct, complete and the same in both versions.
    @testset "gaussianvectors.jl" begin
        m = [10.0, -3.5]
        A = [1.0 3.0; 5.0 7.0]
        B = [1.0 1.0 1.0; 2.0 2.0 2.0 ; 3.0 4.0 5.0]
        C = m
        GV = GaussianVector(A)
        @test mean(GV) == GV.mean == [0.0, 0.0]
        @test var(GV) == GV.var == A * A'
        @test GaussianVector(C).mean ≈ [0.0, 0.0]
        GV = GaussianVector(m, A)
        @test mean(GV) == GV.mean == m
        @test var(GV) == GV.var == A * A'
        @test length(GV) == GV.length == 2
        @test GV.numberofGaussians == 2
        @test_throws DimensionMismatch GaussianVector([0.0], A)
        @test_throws DimensionMismatch GaussianVector([3.4, 2.2, 9.5], A)
        @test_throws DimensionMismatch GaussianVector([0.0], [1.0, 2.0])
        @test_throws DimensionMismatch GaussianVector([0.0, 0.0], [2.0])
        GV = GaussianVector([0.0, 0.0], [1.0, 2.0])
        @test mean(GV) == GV.mean == [0.0, 0.0]
        @test var(GV) == GV.var == [1.0, 2.0] * [1.0, 2.0]'
        @test length(GV) == GV.length == 2
        @test GV.numberofGaussians == 1
        randomvector = rand(GV)
        @test 2.0 * randomvector[1] ≈ randomvector[2]
        m = [1.0, 2.0]
        Σ = [1.0 0.0; 0.0 2.0]
        n = 100000
        GV = GaussianVector(m, Σ)
        nrandomdraws = Array{Float64, 2}(n, 2)
        for k in 1:n
            nrandomdraws[k, :] = rand(GV)
        end
        mempirical = mean(nrandomdraws, 1)
        vempirical = var(nrandomdraws, 1)
        @test norm([1.0 2.0]-mempirical) < 0.1
        @test norm(vempirical-[1.0 4.0]) < 0.1
    end

    #Test functions are correct, complete and the same in both versions.
    @testset "gaussianprocess.jl" begin
        @test Set([GaussianProcess, FaberSchauderExpansionWithGaussianCoefficients]) ⊆ Set(subtypes(AbstractGaussianProcess))

        @test supertype(AbstractGaussianProcess) == Any

        @test_throws DimensionMismatch GaussianProcess([fourier(k) for k in 1:2], GaussianVector(sparse(Diagonal([1.0,1.0,1.0]))))

        X = GaussianProcess([sinpi, cospi], GaussianVector(sparse(Diagonal([1.0, 1.0]))))

        @test X.basis == [sinpi, cospi]

        @test typeof(X.distribution) <: GaussianVector
        @test X.distribution.Σ == [1.0 0.0; 0.0 1.0]
        @test mean(X.distribution) == [0.0,0.0]
        @test length(X) == 2
        @test typeof(rand(X)) <: Function

        n = 100000
        a = sinpi(1/4)
        x = Vector{Float64}(n)
        for k in 1:n
            f = rand(X)
            x[k] = f(1/4)
        end

        @test abs(mean(x))<0.01
        @test abs(var(x)-2*a^2)<0.1

        x = 0.0:0.1:1.0
        y = sumoffunctions([sin, cos], [1., 1.]).(x)
        z = map(x-> sin(x)+cos(x), x)
        @test y ≈ z

        A = sparse(Diagonal([1.0, 1.0]))

        distribution = GaussianVector(A)

        @test_throws AssertionError FaberSchauderExpansionWithGaussianCoefficients(1,distribution)

        @test_throws AssertionError FaberSchauderExpansionWithGaussianCoefficients(Vector{Float64}(0))

        Π = FaberSchauderExpansionWithGaussianCoefficients(0,distribution)
        @test length(Π) == 2

        @test_throws AssertionError sumoffunctions([sin], [1.0,2.0])

        vectoroffunctions = [sin, x -> cos(x+π/2)]
        vectorofscalars = [1.0, 1.0]
        f = sumoffunctions(vectoroffunctions, vectorofscalars)
        x = 0.0:0.01:1.0
        y = f.(x)
        @test maximum(abs.(y)) < 0.001

        n = 10000
        x = Vector{Float64}(n)

        for k in 1:n
            f = rand(Π)
            x[k] = f(0.5)
        end

        @test abs(mean(x))<0.1
        @test abs(var(x)-1.0) <0.1

        for k in 1:n
            f = rand(Π)
            x[k] = f(0.0)
        end

        @test abs(mean(x))<0.1
        @test abs(var(x)-1.0) <0.1

        for k in 1:n
            f = rand(Π)
            x[k] = f(1.0)
        end

        @test abs(mean(x))<0.1
        @test abs(var(x)-1.0) <0.1

        distribution = GaussianVector([1.0, -2.0], A)
        Π = FaberSchauderExpansionWithGaussianCoefficients(0,distribution)

        @test length(Π) == 2

        n = 10000
        x = Vector{Float64}(n)

        for k in 1:n
            f = rand(Π)
            x[k] = f(0.5)
        end

        @test abs(mean(x)+2.0)<0.1
        @test abs(var(x)-1.0) <0.1

        for k in 1:n
            f = rand(Π)
            x[k] = f(0.0)
        end

        @test abs(mean(x)-1.0)<0.1
        @test abs(var(x)-1.0) <0.1

        for k in 1:n
            f = rand(Π)
            x[k] = f(1.0)
        end

        @test abs(mean(x)-1.0)<0.1
        @test abs(var(x)-1.0) <0.1

        Π = FaberSchauderExpansionWithGaussianCoefficients([1.0, sqrt(2)])

        n = 10000
        x = Vector{Float64}(n)
        for k in 1:n
            f = rand(Π)
            x[k] = f(0.25)
        end

        @test abs(mean(x)) < 0.1
        @test abs(var(x) - 2.5) < 0.1

        A = sparse(eye(4))

        distribution = GaussianVector([1.0,1.0,1.0,1.0], A)

        Π = FaberSchauderExpansionWithGaussianCoefficients(1, distribution)

        n = 10000
        x = Vector{Float64}(n)
        for k in 1:n
            f = rand(Π)
            x[k] = f(0.25)
        end

        @test abs(mean(x)-2.0) < 0.1
        @test abs(var(x)-1.5) < 0.1
    end #testset "gaussianprocess.jl"

    #Test functions are correct, complete and the same in both versions.
    @testset "calculateposterior.jl" begin
        # Test for two cases whether calculatedependentfsfunctions
        # calculates the right dependency structure of the Faber-Schauder basis.
        A = [
        true true true true true true true true;
        true true true true true true true true;
        true true true false true true false false;
        true true false true false false true true;
        true true true false true false false false;
        true true true false false true false false;
        true true false true false false true false;
        true true false true false false false true
        ]

        A = Symmetric(A)
        f = BayesianNonparametricStatistics.calculatedependentfaberschauderfunctions
        lengthvectors, rowindices, columnindices = f(2)

        @test lengthvectors == length(rowindices) == length(columnindices)

        B = falses(8, 8)
        for i in 1:lengthvectors
            B[rowindices[i], columnindices[i]] = true
            B[columnindices[i], rowindices[i]] = true
        end
        for i in 1:8
            B[i,i] = true
        end

        @test A == B

        A = [
        true true true true true true true true true true true true true true true true; #1
        true true true true true true true true true true true true true true true true; #2
        true true true false true true false false true true true true false false false false; #3
        true true false true false false true true false false false false true true true true; #4
        true true true false true false false false true true false false false false false false; #5
        true true true false false true false false false false true true false false false false; #6
        true true false true false false true false false false false false true true false false; #7
        true true false true false false false true false false false false false false true true; #8
        true true true false true false false false true false false false false false false false; #9
        true true true false true false false false false true false false false false false false; #10
        true true true false false true false false false false true false false false false false; #11
        true true true false false true false false false false false true false false false false; #12
        true true false true false false true false false false false false true false false false; #13
        true true false true false false true false false false false false false true false false; #14
        true true false true false false false true false false false false false false true false; #15
        true true false true false false false true false false false false false false false true #16
        ]

        A = Symmetric(A)

        B = falses(16, 16)
        lengthvectors, rowindices, columnindices = f(3)

        @test lengthvectors == length(rowindices) == length(columnindices)

        for i in 1:lengthvectors
            B[rowindices[i], columnindices[i]] = true
            B[columnindices[i], rowindices[i]] = true
        end
        for i in 1:16
            B[i,i] = true
        end

        @test A == B

        f = BayesianNonparametricStatistics.calculategirsanovmatrixelement

        @test_throws DimensionMismatch f(collect(1.0:2.0), collect(1.0:3.0),
            collect(1.0:2.0), collect(1.0:2.0))

        @test_throws DimensionMismatch f(collect(1.0:2.0), collect(1.0:3.0),
            collect(1.0:2.0), 3.0)

        @test_throws DimensionMismatch f(collect(1.0:2.0), collect(1.0:3.0),
            2.5, collect(1.0:2.0))

        @test_throws DimensionMismatch f(collect(1.0:2.0), collect(1.0:3.0),
            0.9, 1.2)

        @test_throws DimensionMismatch f(collect(1.0:3.0), collect(1.0:2.0),
            collect(1.0:2.0), collect(1.0:2.0))

        @test_throws DimensionMismatch f(collect(1.0:3.0), collect(1.0:2.0),
            collect(1.0:2.0), 3.0)

        @test_throws DimensionMismatch f(collect(1.0:3.0), collect(1.0:2.0),
            2.5, collect(1.0:2.0))

        @test_throws DimensionMismatch f(collect(1.0:3.0), collect(1.0:2.0),
            0.9, 1.2)

        @test_throws DimensionMismatch f(collect(1.0:2.0), collect(1.0:2.0),
            collect(1.0:3.0), collect(1.0:2.0))

        @test_throws DimensionMismatch f(collect(1.0:2.0), collect(1.0:2.0),
            collect(1.0:3.0), collect(1.0:3.0))

        @test_throws DimensionMismatch f(collect(1.0:2.0), collect(1.0:2.0),
            collect(1.0:3.0), 1.2)

        @test f(collect(1.0:100.0),ones(Float64,100),1.0,1.0) == 5050.0

        @test f(collect(1.0:100.0),2.0*ones(Float64,100),1.0,1.0) == 10100.0

        @test f(1.0:100.0, 1.0:100.0, 1.0:100.0, ones(Float64, 100)) == 100.0

        @test f(1.0:100.0, 1.0:100.0, 1.0:100.0, 1.0) == 100.0

        @test f(1.0:100.0, 1.0:100.0, 1.0, 1.0 ./ (1.0:100.0)) == 5050.0

        @test f(Float64[], Float64[], Float64[], Float64[]) == 0.0

        @test f(Float64[], Float64[], Float64[], 10.0) == 0.0

        @test f(Float64[], Float64[], 2.5, Float64[]) == 0.0

        @test f(Float64[], Float64[], 2.0, 3.0) == 0.0

        samplevalueindices = vcat(trues(100), falses(100))

        @test_throws BoundsError f(samplevalueindices, 1.0:100.0, 1.0:100.0,
            1.0, 1.0)

        @test_throws BoundsError f(samplevalueindices, 1.0:200.0, 1.0:100.0,
            1.0, 1.0)

        @test_throws BoundsError f(samplevalueindices, 1.0:200.0, 1.0:200.0,
            1.0:100.0, 1.0)

        @test_throws BoundsError f(samplevalueindices, 1.0:200.0, 1.0:200.0,
            1.0:200.0, 1.0:100.0)

        @test_throws BoundsError f(samplevalueindices, 1.0:300.0, 1.0:300.0,
            1.0, 1.0)

        @test_throws BoundsError f(samplevalueindices, 1.0:200.0, 1.0:300.0,
            1.0, 1.0)

        @test_throws BoundsError f(samplevalueindices, 1.0:200.0, 1.0:200.0,
            1.0:300.0, 1.0)

        @test_throws BoundsError f(samplevalueindices, 1.0:200.0, 1.0:200.0,
            1.0:200.0, 1.0:300.0)

        @test f(samplevalueindices, 1.0:200.0, ones(Float64, 200), 1.0, 1.0) ==
            5050.0

        @test f(samplevalueindices, 1.0:200.0, ones(Float64, 200),
            ones(Float64,200), 1.0) == 5050.0

        @test f(samplevalueindices, 1.0:200.0, ones(Float64, 200), 1.0,
            ones(Float64,200)) == 5050.0

        @test f(samplevalueindices, 1.0:200.0, ones(Float64, 200),
            ones(Float64, 200), ones(Float64,200)) == 5050.0

        @test f(samplevalueindices, ones(Float64, 200), 1.0:200.0,
            ones(Float64, 200), ones(Float64,200)) == 5050.0

        @test f(samplevalueindices, ones(Float64, 200), 1.0:200.0, 1.0, 1.0) ==
            5050.0

        @test f(samplevalueindices, ones(Float64, 200), 1.0:200.0,
            ones(Float64,200), 1.0) == 5050.0

        @test f(samplevalueindices, ones(Float64, 200), 1.0:200.0, 1.0,
            ones(Float64,200)) == 5050.0

        f = BayesianNonparametricStatistics.calculateΔt

        @test f(0.0:0.1:100.0) == 0.1

        @test f(0.0:0.1:-0.1) == 0.1

        @test f([0.0, 1.0, 100.0]) == [1.0, 99.0]

        f = BayesianNonparametricStatistics.calculategirsanovmatrix

        # bla.

        f = BayesianNonparametricStatistics.calculategirsanovvectorelement

        @test f(1.0:100.0, ones(Float64, 100), ones(Float64,100)) == 5050.0

        @test f(1.0:100.0, ones(Float64, 100), 1.0) == 5050.0

        @test_throws DimensionMismatch f(1.0:100.0, 1.0:50.0, 2.0)

        @test_throws DimensionMismatch f(1.0:66.0, 1.0:100.0, 2.0)

        @test_throws DimensionMismatch f(1.0:100.0, 1.0:50.0, 1.0:50.0)

        @test_throws DimensionMismatch f(1.0:50.0, 1.0:100.0, 1.0:50.0)

        @test_throws DimensionMismatch f(1.0:100.0, 1.0:100.0, 1.0:50.0)

        @test f(1.0:10.0, 1.0:10.0, 0.0) == Inf

        @test f(1.0:2.0, 1.0:2.0, [0.0, 1.0]) == Inf

        samplevalueindices = vcat(trues(100), falses(100))

        @test f(samplevalueindices, collect(1.0:200.0), collect(1.0:200.0).^2, 1.0) == sum(k^3
            for k in 1:100)

        @test f(samplevalueindices, collect(1.0:200.0), collect(1.0:200.0).^2, collect(1.0/k
            for k in 1.0:200.0)) == sum(k^5 for k in 1:100)

        @test_throws BoundsError f(samplevalueindices, 1.0:100.0, 1.0:100.0,
            1.0)

        @test_throws BoundsError f(samplevalueindices, 1.0:100.0, 1.0:100.0,
            1.0:100.0)

        f = BayesianNonparametricStatistics.calculategirsanovvector

        # lengthvector = 2
        samplevalues = [1.0, 2.0]
        ψXt = [[1.0], [1.0]]
        σXt = [1.0]
        @test f(2, samplevalues, ψXt, σXt) == [1.0, 1.0]

        # lengthvector = 2
        samplevalues = [1.0, 2.0]
        ψXt = [[1.0], [1.0]]
        σXt = 1.0
        @test f(2, samplevalues, ψXt, σXt) == [1.0, 1.0]

        samplevalueindices = [BitArray([true, false]), BitArray([false, true])]

        # lengthvector = 2
        samplevalues = [1.0, 2.0, 3.0]
        ψXt = [[1.0, 1.0], [1.0, 1.0]]
        σXt = [1.0, 1.0]
        @test f(2, samplevalueindices, samplevalues, ψXt, σXt) == [1.0, 1.0]

        # lengthvector = 2
        samplevalues = [1.0, 2.0, 3.0]
        ψXt = [[1.0, 1.0], [1.0, 1.0]]
        σXt = 1.0
        @test f(2, samplevalueindices, samplevalues, ψXt, σXt) == [1.0, 1.0]

        f = BayesianNonparametricStatistics.calculateσXt

        @test f(1.0, 1.0:10.0) == 1.0

        @test f(identity, 1.0:10.0) == collect(1.0:10.0)

        α = 0.5
        Π = FaberSchauderExpansionWithGaussianCoefficients([2^(α*j) for j in 1:5])

        model = SDEModel(1.0, 0.0, 10000.0, 0.01)
        sde = SDE(x->0.0, model)
        X = rand(sde)

        postΠ = calculateposterior(Π, X, model)

        n = 1000
        y = 0.0:0.01:1.0
        x = Array{Float64}(length(y),n)
        for k in 1:n
            f = rand(postΠ)
            x[:,k] = f.(y)
        end
        @test maximum(abs.(mean(x,1))) < 0.1

        Π = GaussianProcess([fourier(k) for k in 1:40], GaussianVector(sparse(Diagonal([k^(-1.0) for k in 1:40]))))

        postΠ = calculateposterior(Π, X, model)

        n = 1000
        y = 0.0:0.01:1.0
        x = Array{Float64}(length(y),n)
        for k in 1:n
            f = rand(postΠ)
            x[:,k] = f.(y)
        end

        @test maximum(abs.(mean(x,1))) < 0.1
    end # testset calculateposterior.jl
end # testset.


text:
include("EDPolicy.jl")

export EDManager

## definition
"""
    EDManager(agents::Dict{Any, EDPolicy})

A special MultiAgentManager in which all agents use Exploitability Descent(ED) algorithm to play the game.
"""
mutable struct EDManager <: AbstractPolicy
    agents::Dict{Any, EDPolicy}
end

## interactions
function (π::EDManager)(env::AbstractEnv)
    player = current_player(env)
    if player == chance_player(env)
        rand(legal_action_space(env))
    else
        π.agents[player](env)
    end
end

RLBase.prob(π::EDManager, env::AbstractEnv, args...) = prob(π.agents[current_player(env)], env, args...)

## run function
function Base.run(
    π::EDManager,
    env::AbstractEnv,
    stop_condition = StopAfterEpisode(1),
    hook::AbstractHook = EmptyHook(),
)
    @assert NumAgentStyle(env) == MultiAgent(2) "ED algorithm only support 2-players games."
    @assert UtilityStyle(env) isa ZeroSum "ED algorithm only support zero-sum games."

    is_stop = false

    while !is_stop
        RLBase.reset!(env)
        hook(PRE_EPISODE_STAGE, π, env)

        for (player, policy) in π.agents
            # construct opponent's best response policy.
            oppo_best_response = BestResponsePolicy(π, env, policy.opponent)
            # update player's policy by using policy-gradient.
            update!(policy, oppo_best_response, env, player)
        end

        # run one episode for update stop_condition
        RLBase.reset!(env)
        while !is_terminated(env)
            π(env) |> env
        end

        if stop_condition(π, env)
            is_stop = true
            break
        end
        hook(POST_EPISODE_STAGE, π, env)
    end
    hook(POST_EXPERIMENT_STAGE, π, env)
    hook
end

ext: jl


text:
# This file is part of the ValidatedNumerics.jl package; MIT licensed

# Decorated intervals, following the IEEE 1758-2015 standard

"""
    DECORATION

Enumeration constant for the types of interval decorations.
The nomenclature of the follows the IEEE-1788 (2015) standard
(sect 11.2):

- `com -> 4`: common: bounded, non-empty
- `dac -> 3`: defined (nonempty) and continuous
- `def -> 2`: defined (nonempty)
- `trv -> 1`: always true (no information)
- `ill -> 0`: nai ("not an interval")
"""
@enum DECORATION ill=0 trv=1 def=2 dac=3 com=4
# Note that `isless`, and hence ``<` and `min`, are automatically defined for enums

"""
    DecoratedInterval

A `DecoratedInterval` is an interval, together with a *decoration*, i.e.
a flag that records the status of the interval when thought of as the result
of a previously executed sequence of functions acting on an initial interval.
"""
type DecoratedInterval{T<:Real} <: AbstractInterval
    interval::Interval{T}
    decoration::DECORATION

    function DecoratedInterval(I::Interval, d::DECORATION)
        dd = decoration(I)
        dd <= trv && return new(I, dd)
        d == ill && return new(nai(I), d)
        return new(I, d)
    end
end

DecoratedInterval{T<:AbstractFloat}(I::Interval{T}, d::DECORATION) =
    DecoratedInterval{T}(I, d)
function DecoratedInterval{T<:Real}(a::T, b::T, d::DECORATION)
    a > b && return DecoratedInterval(nai(T), ill)
    DecoratedInterval(Interval(a,b), d)
end
DecoratedInterval{T<:Real}(a::T, d::DECORATION) = DecoratedInterval(Interval(a,a), d)
DecoratedInterval(a::Tuple, d::DECORATION) = DecoratedInterval(a..., d)
DecoratedInterval{T<:Real, S<:Real}(a::T, b::S, d::DECORATION) =
    DecoratedInterval(promote(a,b)..., d)

# Automatic decorations for an interval
DecoratedInterval(I::Interval) = DecoratedInterval(I, decoration(I))
function DecoratedInterval{T<:Real}(a::T, b::T)
    a > b && return DecoratedInterval(nai(T), ill)
    DecoratedInterval(Interval(a,b))
end
DecoratedInterval{T<:Integer}(a::T, b::T) = DecoratedInterval(float(a),float(b))
DecoratedInterval{T<:Real}(a::T) = DecoratedInterval(Interval(a,a))
DecoratedInterval{T<:Real, S<:Real}(a::T, b::S) = DecoratedInterval(promote(a, b)...)
DecoratedInterval(a::Tuple) = DecoratedInterval(a...)

DecoratedInterval(I::DecoratedInterval, dec::DECORATION) = DecoratedInterval(I.interval, dec)

interval_part(x::DecoratedInterval) = x.interval
decoration(x::DecoratedInterval) = x.decoration

# Automatic decorations for an Interval
function decoration(I::Interval)
    isnai(I) && return ill           # nai()
    isempty(I) && return trv         # emptyinterval
    isunbounded(I) && return dac     # unbounded
    com                              # common
end

# Promotion and conversion, and other constructors

promote_rule{T<:Real, N, R<:Real}(::Type{DecoratedInterval{T}},
    ::Type{ForwardDiff.Dual{N,R}}) = ForwardDiff.Dual{N, DecoratedInterval{promote_type(T,R)}}

promote_rule{T<:Real, S<:Real}(::Type{DecoratedInterval{T}}, ::Type{S}) =
    DecoratedInterval{promote_type(T, S)}
promote_rule{T<:Real, S<:Real}(::Type{DecoratedInterval{T}}, ::Type{DecoratedInterval{S}}) =
    DecoratedInterval{promote_type(T, S)}

convert{T<:Real, S<:Real}(::Type{DecoratedInterval{T}}, x::S) =
    DecoratedInterval( Interval(T(x, RoundDown), T(x, RoundUp)) )
convert{T<:Real, S<:Integer}(::Type{DecoratedInterval{T}}, x::S) =
    DecoratedInterval( Interval(T(x), T(x)) )
# function convert{T<:AbstractFloat}(::Type{DecoratedInterval{T}}, x::Float64)
#     convert(DecoratedInterval{T}, rationalize(x))
# end
function convert{T<:Real}(::Type{DecoratedInterval{T}}, xx::DecoratedInterval)
    x = interval_part(xx)
    x = convert(Interval{T},x)
    DecoratedInterval( x, decoration(xx) )
end

# show(io::IO, x::DecoratedInterval) = print(io, x.interval, "_", x.decoration)

macro decorated(ex...)
    local x

    if length(ex) == 1
        x = :(@interval($(esc(ex[1]))))
    else
        x = :($(esc(ex[1])), $(esc(ex[2])))
    end

    :(DecoratedInterval($x))
end


text:
struct InconsistentVRep{T, AT, D<:Polyhedra.FullDim} <: VRepresentation{T}
    points::Polyhedra.PointsHull{T, AT, D}
    rays::Polyhedra.RaysHull{T, AT, D}
    function InconsistentVRep{T, AT, D}(d::Polyhedra.FullDim, points, lines,
                                        rays) where {T, AT, D}
        new{T, AT, D}(Polyhedra.PointsHull(d, points),
                      Polyhedra.RaysHull(d, lines, rays))
    end
end
Polyhedra.FullDim(rep::InconsistentVRep) = Polyhedra.FullDim(rep.points)
Polyhedra.dualtype(::Type{InconsistentVRep{T, AT, D}}, ::Type{AT}) where {T, AT, D} = Polyhedra.Intersection{T, AT, D}
Polyhedra.hvectortype(::Type{<:InconsistentVRep{T, AT}}) where {T, AT} = AT
Polyhedra.vvectortype(::Type{<:InconsistentVRep{T, AT}}) where {T, AT} = AT
Polyhedra.similar_type(PT::Type{<:InconsistentVRep}, d::Polyhedra.FullDim, ::Type{T}) where {T} = InconsistentVRep{T, Polyhedra.similar_type(Polyhedra.hvectortype(PT), d, T), typeof(d)}
Polyhedra.fulltype(::Type{InconsistentVRep{T, AT, D}}) where {T, AT, D} = InconsistentVRep{T, AT, D}
#Polyhedra.@subrepelem InconsistentVRep SymPoint points
Polyhedra.@subrepelem InconsistentVRep Point points
Polyhedra.@subrepelem InconsistentVRep Line rays
Polyhedra.@subrepelem InconsistentVRep Ray rays


text:
using DFTK
using Test

@testset "split_evenly" begin
    function run_test(itr, n)
        splitted = DFTK.split_evenly(itr, n)
        @test collect(itr) == vcat(splitted...)
    end

    run_test(1:12, 4)
    run_test(1:31, 7)
    run_test(1:168, 16)
    run_test(1:14, 16)
end


text:
include("m_step.jl")

function calibrate_item_mmle!(
    item::AbstractItem,
    examinees::Vector{<:AbstractExaminee},
    responses::Vector{Response};
    clean = false,
    int_opt_x_tol_rel::Float64 = 0.0001,
    int_opt_max_time::Float64 = 10.0,
    int_opt_f_tol_rel::Float64 = 0.00001,
)
    if !clean
        responses = sort(filter(r -> r.item_idx == item.idx, responses), by = r -> r.examinee_idx)
        examinees_idx = map(r -> r.examinee_idx, responses)
        examinees = sort(filter( e -> e.idx in examinees_idx, examinees), by = e -> e.idx)
    end
    item.parameters = m_step(
        item.parameters,
        get_latents(examinees),
        map(r -> r.val, responses),
        [Float64(int_opt_max_time), int_opt_x_tol_rel, int_opt_f_tol_rel]
        )
    #pmap( i -> calibrate_item_mmle!(i, examinees, responses; int_opt_x_tol_rel =int_opt_x_tol_rel, int_opt_max_time = int_opt_max_time, int_opt_f_tol_rel = int_opt_f_tol_rel), items)
    return nothing
end

function calibrate_item_mmle!(
    items::Vector{AbstractItem},
    examinees::Vector{<:AbstractExaminee},
    responses::Vector{Response};
    clean = false,
    int_opt_x_tol_rel::Float64 = 0.0001,
    int_opt_max_time::Float64 = 10.0,
    int_opt_f_tol_rel::Float64 = 0.00001,
)
    pmap( i -> calibrate_item_mmle!(
            i,
            examinees,
            responses;
            int_opt_max_time = int_opt_max_time,
            int_opt_x_tol_rel = int_opt_x_tol_rel,
            int_opt_f_tol_rel = int_opt_f_tol_rel,
            clean = clean,
        ),
        items
    )
    return nothing
end


text:
mutable struct inner_problem
	# structure
	block_dims::Vector{Int64}
	m::Int64
	# Data
	A::SparseMatrixCSC{Float64, BlasInt}
	b::Vector{Float64}
	C::Vector{Float64}

	# Solutions
	X::Vector{Float64}
	y::Vector{Float64}
	S::Vector{Float64}

	# Norms
	normsA::Vector{Float64}
	normb::Float64
	normC::Float64

    num_iters::Int64

	# Params
	verbose::Bool
	maxIter::Int64
	ε_relgap::Float64
	ε_φ::Float64

	function inner_problem(block_dims_::Vector{Int64},
						   A_::SparseMatrixCSC{Float64, Int64},
						   b_::Vector{Float64},
						   c_::Vector{Float64})
		prb = new()
		prb.block_dims = block_dims_
		prb.m = length(b_)
		prb.A = A_
		prb.b = b_
		prb.C = c_
		prb.X = Vector{Float64}(undef, 0)
		prb.y = Vector{Float64}(undef, 0)
		prb.S = Vector{Float64}(undef, 0)
		prb.normsA = Vector{Float64}(undef, 0)
		prb.normb = 1.0
		prb.normC = 1.0
		prb.num_iters = zero(Int64)
		prb.verbose = false
		prb.maxIter = 50
		prb.ε_relgap = 1e-6
		prb.ε_φ = 1e-6
		return prb
	end
end

function scale_inner_problem!(prb::inner_problem)
	nblocks = length(prb.block_dims)
	blk_sizes = [(elem > 0) ? elem*(elem+1)÷2 : abs(elem) for elem in [0; prb.block_dims]]
	blk_indices = cumsum(blk_sizes)
	blocks = [blk_indices[j]+1:blk_indices[j+1] for j = 1:nblocks]
	n = sum(blk_sizes)
	N = sum(abs.(prb.block_dims))

    BLAS.scal!(length(prb.b), 1.0 / prb.normb, prb.b, 1)

	@views for j = 1:nblocks
		BLAS.scal!(blk_sizes[j+1], 1.0 / (prb.normC * prb.normsA[j]) , prb.C[blocks[j]], 1)
		BLAS.scal!(blk_sizes[j+1], prb.normsA[j], prb.X[blocks[j]], 1)
		BLAS.scal!(blk_sizes[j+1], 1.0 / (prb.normC * prb.normsA[j]), prb.S[blocks[j]], 1)
	end

	# Scale A
	for j = 1:nblocks
		@. prb.A[blocks[j], :] /= prb.normsA[j]
	end
end

function generate_init!(prb::inner_problem)
	nblocks = length(prb.block_dims)
	ξ = zeros(Float64, nblocks)
	η = zeros(Float64, nblocks)

	temp1 = zeros(Float64, length(prb.b))
	temp2 = zeros(Float64, length(prb.b))

	blk_sizes = [(elem > 0) ? elem*(elem+1)÷2 : abs(elem) for elem in [0; prb.block_dims]] # size of each block, n*(n+1)/2 if nn diag, |n| else
	blk_indices = cumsum(blk_sizes)
	blocks = [blk_indices[j]+1:blk_indices[j+1] for j = 1:nblocks]
	n_elems = sum(blk_sizes)

	#=
	Compute the actual scaling coefficients.
	=#
	@views for j = 1:nblocks
		if prb.block_dims[j] < 0
			for k = 1:length(prb.b)
				temp1[k] = (1.0 + abs(prb.b[k])) / (1.0 + norm(prb.A[blocks[j], k]))
				temp2[k] = norm(prb.A[blocks[j], k])
			end
			ξ[j] = max(1.0, maximum(temp1))
			η[j] = max(1.0, 1.0/sqrt(abs(prb.block_dims[j])) * (1.0 + max(maximum(temp2), norm(prb.C[blocks[j]]))))
		else
			for k = 1:length(prb.b)
				temp1[k] = (1.0 + abs(prb.b[k])) / (1.0 + norm(prb.A[blocks[j], k]))
				temp2[k] = norm(prb.A[blocks[j], k])
			end
			sqrt_s = sqrt(prb.block_dims[j])
			ξ[j] = sqrt_s*max(1.0, sqrt_s*maximum(temp1))
			η[j] = sqrt_s*max(1.0, (1.0/sqrt_s) * (1.0 + max(maximum(temp2), norm(prb.C[blocks[j]]))))
		end
	end

	X0 = zeros(Float64, n_elems)
	S0 = zeros(Float64, n_elems)

	#=
	Initialize the iterates.
	TODO: not the most efficient way, but performed only once.
	=#
	@views for j = 1:nblocks
		if prb.block_dims[j] > 0
			X0[blocks[j]] = svec(ξ[j]*Matrix{Float64}(I, abs(prb.block_dims[j]), abs(prb.block_dims[j])))
			S0[blocks[j]] = svec(η[j]*Matrix{Float64}(I, abs(prb.block_dims[j]), abs(prb.block_dims[j])))
		else
			X0[blocks[j]] = diag(ξ[j]*Matrix{Float64}(I, abs(prb.block_dims[j]), abs(prb.block_dims[j])))
			S0[blocks[j]] = diag(η[j]*Matrix{Float64}(I, abs(prb.block_dims[j]), abs(prb.block_dims[j])))
		end
	end

	prb.X = X0
	prb.y = zeros(Float64, length(prb.b))
	prb.S = S0
end

max_stars_repo_name: TomaszWasTaken/MySDPSolver

max_issues_repo_name: TomaszWasTaken/MySDPSolver

max_forks_repo_name: TomaszWasTaken/MySDPSolver


text:
# Autogenerated wrapper script for rocRAND_jll for x86_64-linux-gnu-cxx11
export libhiprand, librocrand

using hsa_rocr_jll
using ROCmCompilerSupport_jll
using ROCmOpenCLRuntime_jll
using HIP_jll
using rocminfo_jll
JLLWrappers.@generate_wrapper_header("rocRAND")
JLLWrappers.@declare_library_product(libhiprand, "libhiprand.so.1")
JLLWrappers.@declare_library_product(librocrand, "librocrand.so.1")
function __init__()
    JLLWrappers.@generate_init_header(hsa_rocr_jll, ROCmCompilerSupport_jll, ROCmOpenCLRuntime_jll, HIP_jll, rocminfo_jll)
    JLLWrappers.@init_library_product(
        libhiprand,
        "hiprand/lib/libhiprand.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

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

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
ImageCore.HasProperties(::Type{<:AbstractAxisIndices{T,N,P}}) where {T,N,P} = HasProperties(P)
ImageCore.HasProperties(::Type{<:NamedDimsArray{L,T,N,P}}) where {L,T,N,P} = HasProperties(P)
ImageCore.HasProperties(::Type{<:ImageMeta}) = HasProperties{true}()

ImageMetadata.properties(x::AbstractAxisIndices) = _properties(HasProperties(x), x)
ImageMetadata.properties(x::NamedDimsArray) = _properties(HasProperties(x), x)
_properties(::HasProperties{true}, x) = properties(parent(x))
_properties(::HasProperties{false}, x) = error("$(typeof(x)) does not have properties.")

function Base.hasproperty(x::AbstractAxisIndices{T,N,<:ImageMeta}, s::Symbol) where {T,N}
    return hasproperty(parent(x), s)
end
function Base.hasproperty(x::NamedDimsArray{L,T,N,<:AbstractAxisIndices{T,N,<:ImageMeta}}, s::Symbol) where {L,T,N}
    return hasproperty(parent(x), s)
end

function Base.getproperty(x::AbstractAxisIndices{T,N,<:ImageMeta}, s::Symbol) where {T,N}
    return getproperty(parent(x), s)
end
function Base.getproperty(x::NamedDimsArray{L,T,N,<:AbstractAxisIndices{T,N,<:ImageMeta}}, s::Symbol) where {L,T,N}
    return getproperty(parent(x), s)
end

function Base.setproperty!(x::AbstractAxisIndices{T,N,<:ImageMeta}, s::Symbol, val) where {T,N}
    return setproperty!(parent(x), s, val)
end
function Base.setproperty!(x::NamedDimsArray{L,T,N,<:AbstractAxisIndices{T,N,<:ImageMeta}}, s::Symbol, val) where {L,T,N}
    setproperty!(parent(x), s, val)
end


text:
struct MaximumProfile
    version::VersionNumber
    nglyphs::UInt16
    max_points::UInt16
    max_contours::UInt16
    max_composite_points::UInt16
    max_composite_contours::UInt16
    max_zones::UInt16
    max_twilight_points::UInt16
    max_storage::UInt16
    max_function_defs::UInt16
    max_instruction_defs::UInt16
    max_stack_elements::UInt16
    max_size_of_instructions::UInt16
    max_component_elements::UInt16
    max_component_depth::UInt16
end

function Base.parse(io::IO, ::Type{MaximumProfile})
    MaximumProfile(
        version_16_dot_16(read(io, UInt32)),
        (read(io, T) for T in fieldtypes(MaximumProfile)[2:end])...
    )
end


text:
# memory operations
function sync_gpu_to_cpu!(dev::MtlDevice, buf::MtlBuffer{T}) where T
    cmd = MTL.commit!(global_queue(dev)) do cmdbuf
        MtlBlitCommandEncoder(cmdbuf) do enc
            MTL.append_sync!(enc, buf)
        end
    end
    wait(cmd)
end

# GPU -> GPU
Base.unsafe_copyto!(dev::MtlDevice, dst::MtlBuffer{T}, src::MtlBuffer{T}, N::Integer) where T =
    unsafe_copyto!(dev, dst, 1, src, 1, n)
function Base.unsafe_copyto!(dev::MtlDevice, dst::MtlBuffer{T}, doff, src::MtlBuffer{T}, soff, N::Integer) where T
    cmd = MTL.commit!(global_queue(dev)) do cmdbuf
        MtlBlitCommandEncoder(cmdbuf) do enc
            MTL.append_copy!(enc, dst, doff, src, soff, N * sizeof(T))
        end
    end
    wait(cmd)
end

# GPU -> CPU
Base.unsafe_copyto!(dev::MtlDevice, dst::Ptr{T}, src::MtlBuffer{T}, N::Integer) where T =
    unsafe_copyto!(dev, dst, src, 1, n)
function Base.unsafe_copyto!(dev::MtlDevice, dst::Ptr{T}, src::MtlBuffer{T}, soff::Integer, N::Integer) where T
    storage_type = storage(src)
    if storage_type == Private
        tmp_buf = alloc(T, dev, N, storage=Shared)
        Base.unsafe_copyto!(dev, tmp_buf, 1, src, soff, N)
        Base.unsafe_copyto!(dst, content(tmp_buf), N)
        free(tmp_buf)
    elseif storage_type == Shared
        Base.unsafe_copyto!(dst, content(src), N)
    elseif storage_type == Managed
        sync_gpu_to_cpu!(dev, src)
        Base.unsafe_copyto!(dst, content(src), N)
    end
    return dst
end

# CPU -> GPU
Base.unsafe_copyto!(dev::MtlDevice, dst::MtlBuffer{T}, src::Ptr{T}, N::Integer) where T =
    unsafe_copyto!(dev, dst, 1, src, n)
function Base.unsafe_copyto!(dev::MtlDevice, dst::MtlBuffer{T}, doff::Integer,  src::Ptr{T}, N::Integer) where T
    storage_type = storage(dst)

    if storage_type == Private
        # Alloc a buffer containing a copy of src-ptr that is managed
        tmp_buf = alloc(T, dev, N, src, storage=Shared)
        # Copy from GPU to GPU buffer residing in CPU memory
        Base.unsafe_copyto!(dev, dst, doff, tmp_buf, 1, N)
        free(tmp_buf)
    elseif storage_type == Shared
        Base.unsafe_copyto!(dev, content(dst)+(doff-1)*sizeof(T), src, N)
    elseif storage_type == Managed
        Base.unsafe_copyto!(dev, content(dst)+(doff-1)*sizeof(T), src, N)
        MTL.DidModifyRange!(dst, 1:N)
    end
    return dst
end

function unsafe_fill!(dev::MtlDevice, ptr::MtlBuffer{T}, value::Union{UInt8,Int8}, N::Integer) where T
    cmd = MTL.commit!(global_queue(dev)) do cmdbuf
        MtlBlitCommandEncoder(cmdbuf) do enc
            MTL.append_fill!(enc, src, value, N * sizeof(T))
        end
    end
    wait(cmd)
end


text:
using StochasticDiffEq
using Bridge, StaticArrays
using BenchmarkTools, Test
#####

function lorenz!(t,u,du)
    @inbounds du[1] = 10.0(u[2]-u[1])
    @inbounds du[2] = u[1]*(28.0-u[3]) - u[2]
    @inbounds du[3] = u[1]*u[2] - (8/3)*u[3]
end

function σ_lorenz!(t,u,du)
    @inbounds du[1] = 3.0
    @inbounds du[2] = 3.0
    @inbounds du[3] = 3.0
end

n = 5000
T = (0.0, 10.0)
dt = (T[end]-T[1])/n
tt = collect(range(T[1], stop=T[2], length=n+1))
vx = [1.0,0.0,0.0]

prob_sde_lorenz = SDEProblem(lorenz!, σ_lorenz!, vx, T)
sol = StochasticDiffEq.solve(prob_sde_lorenz, StochasticDiffEq.EM(), dt=dt)

const SV = SVector{3, Float64}

mutable struct VLorenz <: ContinuousTimeProcess{Vector{Float64}}
end

Bridge.b(t, u, ::VLorenz) = [10.0(u[2]-u[1]), u[1]*(28.0-u[3]) - u[2], u[1]*u[2] - (8/3)*u[3]]
function Bridge.b!(t, u, du, ::VLorenz)
    @inbounds du[1] = 10.0(u[2]-u[1])
    @inbounds du[2] = u[1]*(28.0-u[3]) - u[2]
    @inbounds du[3] = u[1]*u[2] - (8/3)*u[3]
end
Bridge.σ(t, u, ::VLorenz) = 3.0I
function Bridge.σ!(t, u, dw, dm, ::VLorenz)
    @inbounds dm[1] = 3*dw[1]
    @inbounds dm[2] = 3*dw[2]
    @inbounds dm[3] = 3*dw[3]
end

W = sample(tt, Wiener{SV}())
VW = VSamplePath(tt, Bridge.mat(W.yy))
VX = VSamplePath(tt, copy(VW.yy))
VX2 = VSamplePath(tt, copy(VW.yy))
Bridge.solve!(Bridge.EulerMaruyama!(), VX, vx, VW, VLorenz())
Bridge.solve!(Bridge.EulerMaruyama(), VX2, vx, VW, VLorenz())

######

mutable struct Lorenz <: ContinuousTimeProcess{SV}
end

Bridge.b(s, u, ::Lorenz) = @SVector [10.0(u[2]-u[1]), u[1]*(28.0-u[3]) - u[2], u[1]*u[2] - (8/3)*u[3]           ]
Bridge.σ(t, u, ::Lorenz) = 3.0I

u = SVector{3}(vx)
X = Bridge.euler(u, W, Lorenz())
X2 = copy(X)
X3 = copy(X)

Bridge.solve!(Euler(), X, u, W, Lorenz())
Bridge.euler!(X2, u, W, Lorenz())
invoke(Bridge.solve!, Tuple{Bridge.EulerMaruyama,Any,SVector{3,Float64},Bridge.AbstractPath,Lorenz},
Bridge.EulerMaruyama(), X3, u, W, Lorenz())

@test maximum(abs.(VW.yy - Bridge.mat(W.yy))) < eps()

@test maximum(abs.(VX.yy - Bridge.mat(X.yy))) < eps()
@test maximum(abs.(VX.yy - VX2.yy)) < eps()
@test maximum(norm.(X.yy - X2.yy)) < eps()
@test maximum(norm.(X.yy - X3.yy)) < eps()

######

print("StochasticDiffEq.solve")
@btime StochasticDiffEq.solve(prob_sde_lorenz, StochasticDiffEq.EM(), dt=dt)
print("sample W")
@btime W = sample(tt, Wiener{SV}())
print("allocate X")
@btime X = copy(W)
print("euler! static")
@btime Bridge.euler!(X, u, W, Lorenz())
print("solve! static")
@btime Bridge.solve!(Bridge.EulerMaruyama(), X, u, W, Lorenz())
print("solve! static fallback")
@btime invoke(Bridge.solve!, Tuple{Bridge.EulerMaruyama,Any,SVector{3,Float64},Bridge.AbstractPath,Lorenz},
Bridge.EulerMaruyama(), X3, u, W, Lorenz())
print("solve! out of place fallback (don't do that)")
@btime Bridge.solve!(Bridge.EulerMaruyama(), VX, vx, VW, VLorenz())
print("solve! inplace")
@btime Bridge.solve!(Bridge.EulerMaruyama!(), VX, vx, VW, VLorenz())
;


text:
function abs_cs_safe(x::T) where {T<:Complex}
    return x*sign(real(x))
end

function abs_cs_safe(x)
    return abs(x)
end


text:
N = 10
filenames = readdir(joinpath("test", "data"))
for filename in filenames
	full_filename = Pkg.dir("DataFrames", "test", "data", filename)
	df = benchmark(() -> read_table(full_filename),
	                                "DataFrame I/O",
	                                filename,
	                                N)
	# TODO: Keep permanent record
	print_table(STDOUT, df, ',', '"', false)
end


text:
import Base: -, trace, inv, det, logdet, transpose, ctranspose, vecnorm, kron

############################# Unary sensitivities #############################
push!(ops, DiffOp(:(Base.:-), :(Tuple{DLA.AA}), [true]))
∇(::typeof(-), ::Arg1, p, Y::AA, Ȳ::AA, X::AA) = map(-, Ȳ)

push!(ops, DiffOp(:(Base.trace), :(Tuple{DLA.AM}), [true]))
∇(::typeof(trace), ::Arg1, p, Y::Real, Ȳ::Real, X::AM) = Diagonal(fill!(similar(X), Ȳ))

push!(ops, DiffOp(:(Base.inv), :(Tuple{DLA.AM}), [true]))
∇(::typeof(inv), ::Arg1, p, Y::AM, Ȳ::AM, X::AM) = -transpose(Y) * Ȳ * transpose(Y)

push!(ops, DiffOp(:(Base.det), :(Tuple{DLA.AM}), [true]))
∇(::typeof(det), ::Arg1, p, Y::Real, Ȳ::Real, X::AM) = Y * Ȳ * transpose(inv(X))

push!(ops, DiffOp(:(Base.logdet), :(Tuple{DLA.AM}), [true]))
∇(::typeof(logdet), ::Arg1, p, Y::Real, Ȳ::Real, X::AM) = Ȳ * transpose(inv(X))

push!(ops, DiffOp(:(Base.transpose), :(Tuple{DLA.AVM}), [true]))
∇(::typeof(transpose), ::Arg1, p, Y::AVM, Ȳ::AVM, X::AVM) = transpose(Ȳ)

push!(ops, DiffOp(:(Base.ctranspose), :(Tuple{DLA.AVM}), [true]))
∇(::typeof(ctranspose), ::Arg1, p, Y::AVM, Ȳ::AVM, X::AVM) = ctranspose(Ȳ)

push!(ops, DiffOp(:(Base.vecnorm), :(Tuple{DLA.AA}), [true]))
∇(::typeof(vecnorm), ::Arg1, p, Y::Real, Ȳ::Real, X::AA) = Ȳ ./ Y .* abs2.(X) ./ X

push!(ops, DiffOp(:(Base.vecnorm), :(Tuple{Real}), [true]))
∇(::typeof(vecnorm), ::Arg1, p, Y::Real, Ȳ::Real, X::Real) = Ȳ * sign(X)

############################# Binary sensitivities #############################
push!(ops, DiffOp(:(Base.:*), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(*), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = Ȳ * B'
∇(::typeof(*), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = A' * Ȳ

push!(ops, DiffOp(:(Base.At_mul_B), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(At_mul_B), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = B * Ȳ'
∇(::typeof(At_mul_B), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = A * Ȳ

push!(ops, DiffOp(:(Base.A_mul_Bt), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(A_mul_Bt), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = Ȳ * B
∇(::typeof(A_mul_Bt), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = Ȳ' * A

push!(ops, DiffOp(:(Base.At_mul_Bt), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(At_mul_Bt), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = B' * Ȳ'
∇(::typeof(At_mul_Bt), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = Ȳ' * A'

push!(ops, DiffOp(:(Base.Ac_mul_B), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(Ac_mul_B), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = B * Ȳ'
∇(::typeof(Ac_mul_B), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = A * Ȳ

push!(ops, DiffOp(:(Base.A_mul_Bc), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(A_mul_Bc), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = Ȳ * B
∇(::typeof(A_mul_Bc), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = Ȳ' * A

push!(ops, DiffOp(:(Base.Ac_mul_Bc), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(Ac_mul_Bc), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = B' * Ȳ'
∇(::typeof(Ac_mul_Bc), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = Ȳ' * A'

push!(ops, DiffOp(:(Base.:/), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(/), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = Ȳ / B'
∇(::typeof(/), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = -Y' * (Ȳ / B')

push!(ops, DiffOp(:(Base.At_rdiv_B), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(At_rdiv_B), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = B \ Ȳ'
∇(::typeof(At_rdiv_B), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = -Y' * (B \ Ȳ')'

push!(ops, DiffOp(:(Base.A_rdiv_Bt), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(A_rdiv_Bt), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = Ȳ / B
∇(::typeof(A_rdiv_Bt), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = -(B' \ Ȳ') * Y

push!(ops, DiffOp(:(Base.At_rdiv_Bt), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(At_rdiv_Bt), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = B' \ Ȳ'
∇(::typeof(At_rdiv_Bt), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = -(B' \ Ȳ') * Y

push!(ops, DiffOp(:(Base.Ac_rdiv_B), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(Ac_rdiv_B), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = B \ Ȳ'
∇(::typeof(Ac_rdiv_B), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = -Y' * (B \ Ȳ')'

push!(ops, DiffOp(:(Base.A_rdiv_Bc), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(A_rdiv_Bc), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = Ȳ / B
∇(::typeof(A_rdiv_Bc), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = -(B' \ Ȳ') * Y

push!(ops, DiffOp(:(Base.Ac_rdiv_Bc), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(Ac_rdiv_Bc), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = B' \ Ȳ'
∇(::typeof(Ac_rdiv_Bc), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = -(B' \ Ȳ') * Y

push!(ops, DiffOp(:(Base.:\), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(\), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = -(A' \ Ȳ) * Y'
∇(::typeof(\), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = A' \ Ȳ

push!(ops, DiffOp(:(Base.At_ldiv_B), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(At_ldiv_B), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = -Y *(A \ Ȳ)'
∇(::typeof(At_ldiv_B), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = A \ Ȳ

push!(ops, DiffOp(:(Base.A_ldiv_Bt), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(A_ldiv_Bt), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = -(Ȳ' / A)' * Y'
∇(::typeof(A_ldiv_Bt), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = Ȳ' / A

push!(ops, DiffOp(:(Base.At_ldiv_Bt), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(At_ldiv_Bt), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = -Y * (Ȳ' / A')
∇(::typeof(At_ldiv_Bt), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = Ȳ' / A'

push!(ops, DiffOp(:(Base.Ac_ldiv_B), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(Ac_ldiv_B), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = -Y * (A \ Ȳ)'
∇(::typeof(Ac_ldiv_B), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = A \ Ȳ

push!(ops, DiffOp(:(Base.A_ldiv_Bc), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(A_ldiv_Bc), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = -(Ȳ' / A)' * Y'
∇(::typeof(A_ldiv_Bc), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = Ȳ' / A

push!(ops, DiffOp(:(Base.Ac_ldiv_Bc), :(Tuple{DLA.AVM, DLA.AVM}), [true, true]))
∇(::typeof(Ac_ldiv_Bc), ::Arg1, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = -Y * (Ȳ' / A')
∇(::typeof(Ac_ldiv_Bc), ::Arg2, p, Y::ASVM, Ȳ::ASVM, A::AVM, B::AVM) = Ȳ' / A'

push!(ops, DiffOp(:(Base.vecnorm), :(Tuple{DLA.AA, Real}), [true, true]))
∇(::typeof(vecnorm), ::Arg1, p, Y::Real, Ȳ::Real, A::AA, B::Real) =
    Ȳ .* Y^(1 - B) .* abs.(A).^B ./ A
∇(::typeof(vecnorm), ::Arg2, p, Y::Real, Ȳ::Real, A::AA, B::Real) =
    Ȳ * (Y^(1 - B) * sum(abs.(A).^B .* log.(abs.(A))) - Y * log(Y)) / B

push!(ops, DiffOp(:(Base.vecnorm), :(Tuple{Real, Real}), [true, true]))
∇(::typeof(vecnorm), ::Arg1, p, Y::Real, Ȳ::Real, A::Real, B::Real) = Ȳ * sign(A)
∇(::typeof(vecnorm), ::Arg2, p, Y::Real, Ȳ::Real, A::Real, B::Real) = 0

push!(ops, DiffOp(:(Base.kron), :(Tuple{AM, AM}), [true, true]))
∇(::typeof(kron), ::Type{Val{1}}, p, Y::AM, Ȳ::AM, A::AM, B::AM) =
    ∇(zero(A), kron, Val{1}, p, Y, Ȳ, A, B)
∇(::typeof(kron), ::Type{Val{2}}, p, Y::AM, Ȳ::AM, A::AM, B::AM) =
    ∇(zero(B), kron, Val{2}, p, Y, Ȳ, A, B)
function ∇(Ā::AM, ::typeof(kron), ::Type{Val{1}}, p, Y::AM, Ȳ::AM, A::AM, B::AM)
    @assert size(Ā) == size(A)
    (I, J), (K, L), m = size(A), size(B), length(Y)
    @inbounds for j = reverse(1:J), l = reverse(1:L), i = reverse(1:I)
        āij = Ā[i, j]
        for k = reverse(1:K)
            āij += Ȳ[m] * B[k, l]
            m -= 1
        end
        Ā[i, j] = āij
    end
    return Ā
end
function ∇(B̄::AM, ::typeof(kron), ::Type{Val{2}}, p, Y::AM, Ȳ::AM, A::AM, B::AM)
    @assert size(B̄) == size(B)
    (I, J), (K, L), m = size(A), size(B), length(Y)
    @inbounds for j = reverse(1:J), l = reverse(1:L), i = reverse(1:I)
        aij = A[i, j]
        for k = reverse(1:K)
            B̄[k, l] += aij * Ȳ[m]
            m -= 1
        end
    end
    return B̄
end


text:
"""
    χ²(y::Vector, μ::Vector, σ::Vector)

Calculate reduced χ² of calculated values `y` with
uncertainty `σ` in comparision to expected values `μ`
"""
χ²(y::Vector, μ::Vector, σ::Vector) = sum((y .- μ).^2 ./ σ.^2)

"""
    compton(θ, p)

Compute compton energy at angle `θ` with initial energy `p[1]`
and electron rest mass `p[2]`.
"""
@. compton(θ, p) = p[1] * (1 / (1 + ((p[1]/p[2])*(1 - cos(θ * (π / 180))))))

"""
    scan(startE::Float64, stopE::Float64, startM::Float64,
        stopM::Float64, xData::Vector, yData::Vector, yError::Vector)

Compute chi squared values over a grid in parameter space defined by
`startE`, `stopE`, `startM`, and `stopM`. This returns the new parameters,
a chi squared value for these parameters, and a new, smaller space to cover.
"""
function scan(startE::Float64, stopE::Float64, startM::Float64, stopM::Float64, xData::Vector, yData::Vector, yError::Vector)
    #define box
    stepSizeE = (stopE - startE) / 4
    stepSizeM = (stopM - startM) / 4

    #initialize values
    chiSq = Inf64
    E = 0.0
    mₑ = 0.0
    newEstart = Inf64
    newEend = Inf64
    newMstart = Inf64
    newMend = Inf64

    # at each step, take the mid point and calculate the chi^2 value
    # Then, compare it to other values. If it's less, update E and m_e
    for i in 1:4
        for j in 1:4
            #get parameters from box
            p1 = (startE + stepSizeE * i) - 0.5 * stepSizeE
            p2 = (startM + stepSizeM * j) - 0.5 * stepSizeM
            params = [p1, p2]

            #calculate chi^2
            expected = compton(xData, params)
            chiSqTest = χ²(yData, expected, yError)

            #if claculated value is less than previous, update E and m
            if chiSqTest < chiSq
                chiSq = chiSqTest
                E = p1
                mₑ = p2
                newEstart = startE + stepSizeE * (i - 1)
                newEend = newEstart + stepSizeE
                newMstart = startM + stepSizeM * (j - 1)
                newMend = newMstart + stepSizeM
            end
        end
    end
    return E, mₑ, newEstart, newEend, newMstart, newMend, chiSq
end

"""
    scan_box(startE1::Float64, stopE1::Float64, startM1::Float64,
            stopM1::Float64, xData::Vector, yData::Vector, yError::Vector, tol::Float64)

Find best fit parameters and their errors and compute χ² given the parameters.

Require an upper limit of χ² of `tol::Float64`.
"""
function scan_box(startE1::Float64, stopE1::Float64, startM1::Float64, stopM1::Float64, xData::Vector, yData::Vector, yError::Vector, tol::Float64)

    #initialize values
    E, m, chiSq = Inf64, Inf64, Inf64
    dof = length(xData) - 2
    startE = startE1
    stopE = stopE1
    startM = startM1
    stopM = stopM1

    # keep computing until desired tolerance is reseach
    counter = 0
    while (chiSq / dof) > tol && counter < 10
        E, m, startE, stopE, startM, stopM, chiSq = scan(startE, stopE, startM, stopM, xData, yData, yError)
        counter += 1
    end

    #compute error
    σE = (stopE - startE) / 2
    σM = (stopM - startM) / 2

    return E, σE, m, σM, chiSq/dof
end


text:
@testset "Simulation Build Tests" begin
    problems = create_simulation_build_test_problems(get_template_basic_uc_simulation())
    sequence = SimulationSequence(
        problems = problems,
        feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
        intervals = Dict(
            "UC" => (Hour(24), Consecutive()),
            "ED" => (Hour(1), Consecutive()),
        ),
        feedforward = Dict(
            ("ED", :devices, :ThermalStandard) => SemiContinuousFF(
                binary_source_problem = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        ini_cond_chronology = InterProblemChronology(),
    )
    sim = Simulation(
        name = "test",
        steps = 1,
        problems = problems,
        sequence = sequence,
        simulation_folder = mktempdir(cleanup = true),
    )

    build_out = build!(sim)
    @test build_out == PSI.BuildStatus.BUILT

    @test isempty(values(sim.internal.simulation_cache))
    for field in fieldnames(SimulationSequence)
        if fieldtype(SimulationSequence, field) == Union{Dates.DateTime, Nothing}
            @test !isnothing(getfield(sim.sequence, field))
        end
    end
    @test isa(sim.sequence, SimulationSequence)

    @test length(findall(x -> x == 2, sequence.execution_order)) == 24
    @test length(findall(x -> x == 1, sequence.execution_order)) == 1
end

@testset "Simulation with provided initial time" begin
    problems = create_simulation_build_test_problems(get_template_basic_uc_simulation())
    sequence = SimulationSequence(
        problems = problems,
        feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
        intervals = Dict(
            "UC" => (Hour(24), Consecutive()),
            "ED" => (Hour(1), Consecutive()),
        ),
        feedforward = Dict(
            ("ED", :devices, :ThermalStandard) => SemiContinuousFF(
                binary_source_problem = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        ini_cond_chronology = InterProblemChronology(),
    )
    second_day = DateTime("1/1/2024  23:00:00", "d/m/y  H:M:S") + Hour(1)
    sim = Simulation(
        name = "test",
        steps = 1,
        problems = problems,
        sequence = sequence,
        simulation_folder = mktempdir(cleanup = true),
        initial_time = second_day,
    )
    build_out = build!(sim)
    @test build_out == PSI.BuildStatus.BUILT

    for (_, problem) in PSI.get_problems(sim)
        @test PSI.get_initial_time(problem) == second_day
    end
end

@testset "Negative Tests (Bad Parametrization)" begin
    problems = create_simulation_build_test_problems(get_template_basic_uc_simulation())
    sequence = SimulationSequence(
        problems = problems,
        feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
        intervals = Dict(
            "UC" => (Hour(24), Consecutive()),
            "ED" => (Hour(1), Consecutive()),
        ),
        feedforward = Dict(
            ("ED", :devices, :ThermalStandard) => SemiContinuousFF(
                binary_source_problem = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        ini_cond_chronology = InterProblemChronology(),
    )

    @test_throws UndefKeywordError sim = Simulation(name = "test", steps = 1)

    sim = Simulation(
        name = "test",
        steps = 1,
        problems = problems,
        sequence = sequence,
        simulation_folder = mktempdir(cleanup = true),
        initial_time = Dates.now(),
    )

    @test_throws IS.ConflictingInputsError build!(sim)

    sim = Simulation(
        name = "fake_path",
        steps = 1,
        problems = problems,
        sequence = sequence,
        simulation_folder = "fake_path",
    )

    @test_throws IS.ConflictingInputsError PSI._check_folder(sim)

    sequence.feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 30))
    sim = Simulation(
        name = "look_ahead",
        steps = 1,
        problems = problems,
        sequence = sequence,
        simulation_folder = mktempdir(cleanup = true),
    )
    @test_throws IS.ConflictingInputsError PSI._check_feedforward_chronologies(sim)

    sequence.feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24))
    sim = Simulation(
        name = "look_ahead",
        steps = 5,
        problems = problems,
        sequence = sequence,
        simulation_folder = mktempdir(cleanup = true),
    )
    @test_throws IS.ConflictingInputsError build!(sim)

    sequence.feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 30))
    sim = Simulation(
        name = "look_ahead",
        steps = 1,
        problems = problems,
        sequence = sequence,
        simulation_folder = mktempdir(cleanup = true),
    )
    @test_throws IS.ConflictingInputsError PSI._check_feedforward_chronologies(sim)

    @test_throws IS.ConflictingInputsError sim = Simulation(
        name = "disconnected problems",
        steps = 5,
        problems = create_simulation_build_test_problems(),
        sequence = sequence,
        simulation_folder = mktempdir(cleanup = true),
    )
end

# Pending tests to update

# @testset "Test Creation of Simulations with Cache" begin
#     stages_definition = create_stages(template_standard_uc, c_sys5_uc, c_sys5_ed)
#
#     # Cache is not defined all together
#     sequence_no_cache = SimulationSequence(
#         step_resolution = Hour(24),
#         order = Dict(1 => "UC", 2 => "ED"),
#         feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
#         horizons = Dict("UC" => 24, "ED" => 12),
#         intervals = Dict(
#             "UC" => (Hour(24), Consecutive()),
#             "ED" => (Hour(1), Consecutive()),
#         ),
#         feedforward = Dict(
#             ("ED", :devices, :Generators) => SemiContinuousFF(
#                 binary_source_stage = PSI.ON,
#                 affected_variables = [PSI.ACTIVE_POWER],
#             ),
#         ),
#         ini_cond_chronology = InterProblemChronology(),
#     )
#     sim = Simulation(
#         name = "cache",
#         steps = 1,
#         stages = stages_definition,
#         stages_sequence = sequence_no_cache,
#         simulation_folder = file_path,
#     )
#     build!(sim)
#
#     @test !isempty(sim.internal.simulation_cache)
#
#     stages_definition = create_stages(template_standard_uc, c_sys5_uc, c_sys5_ed)
#     sequence = SimulationSequence(
#         step_resolution = Hour(24),
#         order = Dict(1 => "UC", 2 => "ED"),
#         feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
#         horizons = Dict("UC" => 24, "ED" => 12),
#         intervals = Dict(
#             "UC" => (Hour(24), Consecutive()),
#             "ED" => (Hour(1), Consecutive()),
#         ),
#         feedforward = Dict(
#             ("ED", :devices, :Generators) => SemiContinuousFF(
#                 binary_source_stage = PSI.ON,
#                 affected_variables = [PSI.ACTIVE_POWER],
#             ),
#         ),
#         cache = Dict(("UC",) => TimeStatusChange(PSY.ThermalStandard, PSI.ON)),
#         ini_cond_chronology = InterProblemChronology(),
#     )
#     sim = Simulation(
#         name = "caches",
#         steps = 2,
#         stages = stages_definition,
#         stages_sequence = sequence,
#         simulation_folder = file_path,
#     )
#
#     build!(sim)
#
#     @test !isempty(sim.internal.simulation_cache)
#
#     stages_definition = create_stages(template_standard_uc, c_sys5_uc, c_sys5_ed)
#     # Uses IntraProblem but the cache is defined in the wrong stage
#     sequence_bad_cache = SimulationSequence(
#         step_resolution = Hour(24),
#         order = Dict(1 => "UC", 2 => "ED"),
#         feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
#         horizons = Dict("UC" => 24, "ED" => 12),
#         intervals = Dict(
#             "UC" => (Hour(24), Consecutive()),
#             "ED" => (Hour(1), Consecutive()),
#         ),
#         feedforward = Dict(
#             ("ED", :devices, :Generators) => SemiContinuousFF(
#                 binary_source_stage = PSI.ON,
#                 affected_variables = [PSI.ACTIVE_POWER],
#             ),
#         ),
#         cache = Dict(("ED",) => TimeStatusChange(PSY.ThermalStandard, PSI.ON)),
#         ini_cond_chronology = IntraProblemChronology(),
#     )
#
#     sim = Simulation(
#         name = "test",
#         steps = 1,
#         stages = stages_definition,
#         stages_sequence = sequence_bad_cache,
#         simulation_folder = file_path,
#     )
#     @test_throws IS.InvalidValue build!(sim)
# end


text:
import Base.getindex

struct SharedSparseMatrixCSC{Tv, Ti<:Integer} <: AbstractSparseMatrix{Tv, Ti}
    m::Int
    n::Int
    colptr::SharedVector{Ti}
    rowval::SharedVector{Ti}
    nzval::SharedVector{Tv}

    function SharedSparseMatrixCSC{Tv,Ti}(m::Integer, n::Integer, colptr::SharedVector{Ti}, rowval::SharedVector{Ti},
                                    nzval::SharedVector{Tv}) where {Tv,Ti<:Integer}
        m < 0 && throw(ArgumentError("number of rows (m) must be ≥ 0, got $m"))
        n < 0 && throw(ArgumentError("number of columns (n) must be ≥ 0, got $n"))
        new(Int(m), Int(n), colptr, rowval, nzval)
end
end

function SharedSparseMatrixCSC(m::Integer, n::Integer, colptr::SharedVector, rowval::SharedVector, nzval::SharedVector)
    Tv = eltype(nzval)
    Ti = promote_type(eltype(colptr), eltype(rowval))
    SharedSparseMatrixCSC{Tv,Ti}(m, n, colptr, rowval, nzval)
end

Base.size(S::SharedSparseMatrixCSC) = (S.m, S.n)

function A_mul_B!(alpha::Number, A::SharedSparseMatrixCSC, B::StridedVecOrMat, beta::Number, C::StridedVecOrMat)
    A.n == size(B, 1) || throw(DimensionMismatch(@sprintf("A.n=%d B.m=%d", A.n, size(B,1))))
    A.m == size(C, 1) || throw(DimensionMismatch(@sprintf("A.m=%d C.m=%d", A.m, size(C,1))))
    size(B, 2) == size(C, 2) || throw(DimensionMismatch(@sprintf("B.n=%d C.n=%d", size(B,2), size(C,2))))
    nzv = A.nzval
    rv = A.rowval
    if beta != 1
        beta != 0 ? scale!(C, beta) : fill!(C, zero(eltype(C)))
    end
    for k = 1:size(C, 2)
        for col = 1:A.n
            alphaxj = alpha*B[col,k]
            @inbounds for j = A.colptr[col]:(A.colptr[col + 1] - 1)
                C[rv[j], k] += nzv[j]*alphaxj
            end
        end
    end
    C
end

A_mul_B!(C::StridedVecOrMat, A::SharedSparseMatrixCSC, B::StridedVecOrMat) = A_mul_B!(one(eltype(B)), A, B, zero(eltype(C)), C)

convert(::Type{SharedSparseMatrixCSC}, M::Matrix) = sparse(M)
convert(::Type{SharedSparseMatrixCSC}, M::AbstractMatrix{Tv}) where {Tv} = convert(SharedSparseMatrixCSC{Tv,Int}, M)
convert(::Type{SharedSparseMatrixCSC{Tv}}, M::AbstractMatrix{Tv}) where {Tv} = convert(SharedSparseMatrixCSC{Tv,Int}, M)
function convert(::Type{SharedSparseMatrixCSC{Tv,Ti}}, M::AbstractMatrix) where {Tv,Ti}
    (I, J, V) = findnz(M)
    eltypeTiI = Base.convert(Vector{Ti}, I)
    eltypeTiJ = Base.convert(Vector{Ti}, J)
    eltypeTvV = Base.convert(Vector{Tv}, V)
    return shared_sparse_IJ_sorted(eltypeTiI, eltypeTiJ, eltypeTvV, size(M)...)
end

shared_sparse(A::AbstractMatrix{Tv}) where {Tv} = convert(SharedSparseMatrixCSC{Tv,Int}, A)

shared_sparse_IJ_sorted(I,J,V,m,n) = shared_sparse_IJ_sorted(I,J,V,m,n,+)
function shared_sparse_IJ_sorted(I::AbstractVector{Ti}, J::AbstractVector{Ti},
                                  V::AbstractVector{Tv},
                           m::Integer, n::Integer, combine::Function) where Ti<:Integer where Tv
    m = m < 0 ? 0 : m
    n = n < 0 ? 0 : n
    if isempty(V); return spzeros(eltype(V),Ti,m,n); end

    cols = zeros(Ti, n+1)
    cols[1] = 1  # For cumsum purposes
    cols[J[1] + 1] = 1

    lastdup = 1
    ndups = 0
    I_lastdup = I[1]
    J_lastdup = J[1]
    L = length(I)

    @inbounds for k=2:L
        if I[k] == I_lastdup && J[k] == J_lastdup
            V[lastdup] = combine(V[lastdup], V[k])
            ndups += 1
        else
            cols[J[k] + 1] += 1
            lastdup = k-ndups
            I_lastdup = I[k]
            J_lastdup = J[k]
            if ndups != 0
                I[lastdup] = I_lastdup
                V[lastdup] = V[k]
            end
        end
    end

    colptr = cumsum!(similar(cols), cols)

    # Allow up to 20% slack
    if ndups > 0.2*L
        numnz = L-ndups
        deleteat!(I, (numnz+1):L)
        deleteat!(V, (numnz+1):length(V))
    end

    return SharedSparseMatrixCSC(m, n, SharedVector{Ti}(colptr), SharedVector{Ti}(I), SharedVector{Tv}(V))
end

function getindex(A::SharedSparseMatrixCSC{T}, i0::Int, i1::Int) where T
    if !(1 <= i0 <= A.m && 1 <= i1 <= A.n); throw(BoundsError()); end
    first = A.colptr[i1]
    last = A.colptr[i1+1]-1
    while first <= last
        mid = (first + last) >> 1
        t = A.rowval[mid]
        if t == i0
            return A.nzval[mid]
        elseif t > i0
            last = mid - 1
        else
            first = mid + 1
        end
    end
    return zero(T)
end


text:
#

export Pose2AprilTag4Corners, PackedPose2AprilTag4Corners
export generateCostAprilTagsPreimageCalib

# suppressing but listing some internal functions
# export _defaultCameraCalib, _AprilTagToPose2

"""
    $TYPEDEF

Simplified constructor type to convert between 4 corner detection of AprilTags to a Pose2Pose2 factor for use in 2D

Notes
- Coordinate frames are:
  - assume robotics body frame is xyz <==> fwd-lft-up
  - assume AprilTags pose is xyz <==> rht-dwn-fwd
  - assume camera frame is xyz <==> rht-dwn-fwd
  - assume Images.jl frame is row-col <==> i-j  <==> dwn-rht
- Helper constructor uses `f_width, f_height, c_width, c_height,s` to build `K`,
  - setting `K` will overrule `f_width,f_height, c_width, c_height,s`.
- Finding preimage from deconv measurement sample `idx` in place of MvNormal mean:
  - see [`generateCostAprilTagsPreimageCalib`](@ref) for detauls.

Example

```julia
# bring in the packages
using AprilTags, Caesar, FileIO

# the size of the tag, as in the outer length of each side on of black square
taglength = 0.15

# load the image
img = load("photo.jpg")

# the image size
width, height = size(img)
# auto-guess `f_width=height, c_width=round(Int,width/2), c_height=round(Int, height/2)`

detector = AprilTagDetector()
tags = detector(img)

# new factor graph with Pose2 `:x0` and a Prior.
fg = generateCanonicalFG_ZeroPose(varType=Pose2)

# use a construction helper to add factors to all the tags
for tag in tags
  tagSym = Symbol("tag\$(tag.id)")
  exists(fg, tagSym) ? nothing : addVariable!(fg, tagSym, Pose2)
  pat = Pose2AprilTag4Corners(corners=tag.p, homography=tag.H, taglength=taglength)
  addFactor!(fg, [:x0; tagSym], pat)
end

# free AprilTags library memory
freeDetector!(detector)
```

DevNotes
- TODO IIF will get plumbing to combine many of preimage `obj` terms into single calibration search

Related

`AprilTags.detect`, `PackedPose2AprilTag4Corners`, [`generateCostAprilTagsPreimageCalib`](@ref)
"""
struct Pose2AprilTag4Corners{T <: SamplableBelief, F <: Function} <: IIF.AbstractManifoldMinimize
  # 4 corners as detected by AprilTags
  corners::NTuple{4,Tuple{Float64,Float64}}
  # homography matrix
  homography::Matrix{Float64}
  # camera calibration
  K::Matrix{Float64}
  # define the size of the tag in meters
  taglength::Float64
  # tag id
  id::Int
  # internally computed relative factor between binary variables-- from camera to tag in camera or body frame
  Zij::RoME.Pose2Pose2{T}

  # experimental for 'SLAM-aware' camera intrinsic calibration -- stores (lambda function, starting estimate), see [`generateCostAprilTagsPreimageCalib`](@ref).
  preimage::Tuple{F, Vector{Float64}}
end

getManifold(::IIF.InstanceType{<:Pose2AprilTag4Corners}) = getManifold(Pose2Pose2) # SE2E2_Manifold

"""
    $SIGNATURES

A default camera calibration matrix, useful for reference and quick tests.

Notes
- Hartley, Zisserman, Multiple View Geometry, p.157

DevNotes
- TODO likely to move to a common camera models Pkg
- TODO consolidate with similar functions in RoME?
"""
function _defaultCameraCalib(;f_width::Real = 480.0,
                              f_height::Real = 480.0,
                              c_width::Real = 320.0,
                              c_height::Real = 240.0,
                              s::Real  = 0.0)
  #
  K = [f_width s  c_width;
        0 f_height c_height;
        0 0  1.0]
  #
  return K
end

"""
    $SIGNATURES

Standardizing a function to convert a regular AprilTags.jl sighting into Pose2 format

Notes
- assume robotics body frame is xyz <==> fwd-lft-up
- assume AprilTags pose is xyz <==> rht-dwn-fwd
- assume camera frame is xyz <==> rht-dwn-fwd
- assume Images.jl frame is row-col <==> i-j  <==> dwn-rht
"""
function _AprilTagToPose2(corners,
                          homography::AbstractMatrix{<:Real},
                          f_width::Real,
                          f_height::Real,
                          c_width::Real,
                          c_height::Real,
                          taglength_::Real)
  #
  # pose <==> cTt
  pose, err1 = AprilTags.tagOrthogonalIteration(corners, homography, f_width, f_height, c_width, c_height, taglength=taglength_)
  cVt = Translation((pose[1:3,4])...)
  # bRc = bRy * yRc
  bRc = _Rot.UnitQuaternion(1/sqrt(2),0,1/sqrt(2),0) * _Rot.UnitQuaternion(1/sqrt(2),0,0,-1/sqrt(2))
  # for tag in body frame == bTt
  bTt = LinearMap(bRc) ∘ cVt

  # also extract the relative yaw between body and tag
  ld = LinearAlgebra.cross(bTt.linear*pose[1:3,1:3]*[0;0;1], [0;0;1])
  theta = TU.wrapRad(atan(ld[2],ld[1]) + pi/2)

  M = getManifold(Pose2Pose2)
  e0 = identity_element(M)

  # create a tangent vector of the measurement
  hat( M, e0, [bTt.translation[1:2,];theta] )
end

const _CornerVecTuple = Union{NTuple{4,Tuple{Float64,Float64}}, <:AbstractVector{Tuple{Float64,Float64}}, <:AbstractVector{<:AbstractVector{<:Real}}}

function Pose2AprilTag4Corners(;corners::_CornerVecTuple=((0.0,0.0),(1.0,0.0),(0.0,1.0),(1.0,1.0)),
                                homography::AbstractMatrix{<:Real}=diagm(ones(3)),
                                f_width::Real=480.0,
                                f_height::Real=f_width,
                                c_width::Real=320.0,
                                c_height::Real=240.0,
                                s::Real =0.0,
                                K::AbstractMatrix{<:Real}=_defaultCameraCalib(f_width=f_width,
                                                                              f_height=f_height,
                                                                              c_width=c_width,
                                                                              c_height=c_height,
                                                                              s=s),
                                id::Int=-1,
                                taglength::Real=0.25,
                                covariance::Matrix{<:Real}=diagm([0.1;0.1;0.1].^2) )
  #
  # make standard tuple
  corners_ = if corners isa Tuple
    corners
  else
    c = corners
    ((c[1][1],c[1][2]),(c[2][1],c[2][2]),(c[3][1],c[3][2]),(c[4][1],c[4][2]))
  end

  # calculate the transform
  # f_width, f_height, c_width, c_height = K[1,1], K[2,2], K[1,3], K[2,3]
  x0 = [K[1,1], K[2,2], K[1,3], K[2,3], taglength]
  Dtag = _AprilTagToPose2(corners, homography, x0...)
  # println("$(tag.id), $(ld[3]), $(round.(Dtag, digits=3))")
  M = SpecialEuclidean(2)
  e0 = identity_element(M)
  p2l2 = Pose2Pose2(MvNormal(vee(M,e0,Dtag), covariance))

  # preimage optimize function
  ## TODO not the greatest use of Manifolds.jl, can do better
  preImgFnc = (y, x) -> sum( vee(M, e0, y - _AprilTagToPose2(corners_, homography, x...)).^2 )
  # preImgFnc = (y, x) -> sum( (y - _AprilTagToPose2(corners_, homography, x...)).^2 )

  #
  return Pose2AprilTag4Corners(corners_, homography, K, taglength, id, p2l2, (preImgFnc, x0))
end

function Pose2AprilTag4Corners( tag::AprilTag;
                                corners::_CornerVecTuple=tag.p,
                                homography::AbstractMatrix{<:Real}=tag.H,
                                f_width::Real=480.0,
                                f_height::Real=f_width,
                                c_width::Real=320.0,
                                c_height::Real=240.0,
                                s::Real =0.0,
                                K::AbstractMatrix{<:Real}=_defaultCameraCalib(f_width=f_width,
                                                                              f_height=f_height,
                                                                              c_height=c_height,
                                                                              c_width=c_width,
                                                                              s=s),
                                id::Int=tag.id,
                                taglength::Real=0.25,
                                covariance::Matrix{<:Real}=diagm([0.1;0.1;0.1].^2) )
  #
  # this is a lazy construction helper for direct use in Julia

  Pose2AprilTag4Corners(;corners=corners,
                        homography=homography,
                        f_width=f_width,
                        f_height=f_height,
                        c_width=c_width,
                        c_height=c_height,
                        s=s,
                        K=K,
                        id=id,
                        taglength=taglength,
                        covariance=covariance )
end

# just pass through sample and factor evaluations

function getSample( pat4c::CalcFactor{<:Pose2AprilTag4Corners} )
  #
  M = getManifold(pat4c.factor.Zij)
  ϵ = getPointIdentity(Pose2)

  X = sampleTangent(M, pat4c.factor.Zij.z, ϵ)
  return X
end

function (pat4c::CalcFactor{<:Pose2AprilTag4Corners})(X,
                                                      p,
                                                      q)
  #

  @assert X isa ProductRepr "Pose2AprilTag4Corners expects measurement sample X to be a Manifolds tangent vector, not coordinate or point representation.  Got X=$X"
  M = getManifold(pat4c.factor.Zij)
  q̂ = Manifolds.compose(M, p, exp(M, identity_element(M, p), X)) #for groups
  #TODO allocalte for vee! see Manifolds #412, fix for AD
  Xc = zeros(3)
  vee!(M, Xc, q, log(M, q, q̂))
  return Xc
  # wXjhat = SE2(wxi)*SE2(z)
  # jXjhat = SE2(wxj) \ wXjhat
  # se2vee(jXjhat)
end

## serialization

struct PackedPose2AprilTag4Corners <: AbstractPackedFactor
  # format of serialized data
  mimeType::String
  # corners, as detected by AprilTags library
  corners::Vector{Float64}
  # homography matrix
  homography::Vector{Float64}
  # camera calibration
  K::Vector{Float64}
  # define the size of the tag in meters
  taglength::Float64
  # tag id
  id::Int
end

function convert( ::Type{<:PackedInferenceType},
                  obj::Pose2AprilTag4Corners)
  #
  corVec = zeros(8)
  corVec[1] = obj.corners[1][1]
  corVec[2] = obj.corners[1][2]
  corVec[3] = obj.corners[2][1]
  corVec[4] = obj.corners[2][2]
  corVec[5] = obj.corners[3][1]
  corVec[6] = obj.corners[3][2]
  corVec[7] = obj.corners[4][1]
  corVec[8] = obj.corners[4][2]

  return PackedPose2AprilTag4Corners( "/application/JuliaLang/PackedPose2AprilTag4Corners",
                                      corVec,
                                      obj.homography[:],
                                      obj.K[:],
                                      obj.taglength,
                                      obj.id )
end

function convert( ::Type{<:DFG.AbstractRelative},
                  obj::PackedPose2AprilTag4Corners)
  #
  cv = obj.corners
  corners = ((cv[1],cv[2]),(cv[3],cv[4]),(cv[5],cv[6]),(cv[7],cv[8]))

  return Pose2AprilTag4Corners( corners=corners,
                                homography=reshape(obj.homography,:,3),
                                K=reshape(obj.K,3,3),
                                taglength=obj.taglength,
                                id=obj.id  )
end

## calibrate via preimage

"""
    $SIGNATURES

Helper function to generate and calculate the aggregate cost indicative of the discrepancy between
the deconvolution prediction and prescribed measurement of mutliple `Pose2AprilTag4Corners` factors
in the factor graph.

The idea is that a bad calibration will give some kind of SLAM result, and if there is enough information
then the SLAM result can be used to bootstrap better and better calibration estimates of the camera that
was used to capture the AprilTag sightings.  This function is meant to help do that secondary parameter
search inside a factor graph objection, after a regular solution has been found.

Notes
- `pred, _ = approxDeconv(dfg, fct)`
- `fct.preimage[1](pred[:,idx], [f_width, f_height, c_width, c_height, taglength])`
  - `fct.preimage[1]` is a function to find the preimage.
- `obj = (fc_wh) -> fct.preimage[1](pred[:,idx], fc_wh)`
  - `fc_wh = [f_width, f_height, c_width, c_height, 0.172]`
- `obj2 = (fcwh) -> obj([fcwh[1]; fcwh[1]; fcwh[2]; c_height; taglength])`
- `result = Optim.optimize(obj, fct.preimage[2], BFGS(), Optim.options(x_tol=1e-8))`
  - A stored starting estimate for optimization `fct.preimage[2]`
"""
function generateCostAprilTagsPreimageCalib(dfg::AbstractDFG,
                                            fsyms::Vector{Symbol}=lsf(dfg, Pose2AprilTag4Corners );
                                            idx::Int = 1, # the sample number
                                            f_width::Real=480.0,
                                            f_height::Real=f_width,
                                            c_width::Real=320.0,
                                            c_height::Real=240.0,
                                            taglength::Real=0.172,
                                            args::Function = (f_cw)->[f_cw[1]; f_cw[1]; f_cw[2]; c_height; taglength]  )
  #

  # temporary memory containers for all the different lambda functions
  fcts = []
  objs = []
  obj2s = []
  preds = []

  for i in 1:length(fsyms)
    fsym = fsyms[i]
    pred, _ = approxDeconv(dfg, fsym)
    push!(preds, pred)
    fct = getFactorType(dfg, fsym)
    push!(fcts, fct)
    # fct.preimage[1](pred[:,idx], f_cw)
    obj = (fc_wh) -> fcts[i].preimage[1](preds[i][:,idx], fc_wh)
    push!(objs, deepcopy(obj))
    # use keyword args mapping on which parameter should be optimized
    obj2 = (f_cw) -> objs[i](args(f_cw))
    push!(obj2s, deepcopy(obj2))
  end

  # set up the cummulative cost over all functions in `obj2s`
  cost = x->((f->f(x)).(obj2s) |> sum)

  # test the function is working
  # @show cost([f_width; c_width])

  return cost
end

##

#


text:
# using Bridge, StaticArrays, Distributions
# using Bridge:logpdfnormal
# using Test, Statistics, Random, LinearAlgebra
# using Bridge.Models
# using DelimitedFiles, DataFrames, RCall
# using CSV
# # install.packages(ggforce)
# using Base.Iterators, SparseArrays, LowRankApprox, Trajectories
# using ForwardDiff
# using Plots,  PyPlot #using Makie
# using RecursiveArrayTools
# using NPZ # for reading python datafiles
# using Random

using BridgeLandmarks

workdir = @__DIR__
println(workdir)
cd(workdir)

using RCall
using Plots
using Random
using Distributions
using NPZ # for reading python datafiles
using DataFrames
using DelimitedFiles
using CSV

pyplot()

include("generatedata.jl")
include("plotting.jl")

Random.seed!(3)

################################# start settings #################################
n = 6  # nr of landmarks
models = [:ms, :ahs]
model = models[1]
println("model: ",model)

ITER = 5
subsamples = 0:1:ITER

showplotσq = false # only for ahs model

sampler =[:sgd, :mcmc][2]
println("sampler: ",sampler)

#------------------------------------------------------------------
# prior on θ = (a, c, γ)
prior_a = Exponential(5.0)
prior_c = Exponential(5.0)
prior_γ = Exponential(5.0)
#------------------------------------------------------------------

datasets =["forwardsimulated", "shifted","shiftedextreme",
            "bear","heart","peach",
            "generatedstefan","forwardsimulated_multiple", "cardiac"]
dataset = datasets[7]
println("dataset: ",dataset)

fixinitmomentato0 = true#false
#------------------------------------------------------------------
# for sgd (FIX LATER)
ϵ = 0.01  # sgd step size
ϵstep(i) = 1/(1+i)^(0.7)

#------------------------------------------------------------------
### MCMC tuning pars
ρ = 0.8              # pcN-step

# proposal for θ = (a, c, γ)
σ_a = 0.2  # update a to aᵒ as aᵒ = a * exp(σ_a * rnorm())
σ_c = 0.2  # update c to cᵒ as cᵒ = c * exp(σ_c * rnorm())
σ_γ = 0.2  # update γ to γᵒ as γᵒ = γ * exp(σ_γ * rnorm())

η(n) = min(0.1, 10/sqrt(n))  # adaptation rate for adjusting tuning pars

#------------------------------------------------------------------
σobs = 0.01   # noise on observations

#------------------------------------------------------------------
# set time grids
dt = 0.01
T = 1.0; t = 0.0:dt:T; tt_ =  BridgeLandmarks.tc(t,T)

outdir = "./figs/"
if false # to use later on, when plotting is transferred from R to Julia
    outdir = "./figs/"* string(model) * "_" * string(sampler) *"_" * string(dataset) * "/"
    if !isdir(outdir) mkdir(outdir) end
end
################################# end settings #################################

### Specify landmarks models
a = 2.0     # Hamiltonian kernel parameter (the larger, the stronger landmarks behave similarly)
c = 0.1     # multiplicative factor in kernel
γ = 1.0     # Noise level

# a = 0.2
# c = 0.1
# γ = 0.002

if model == :ms
    λ = 0.0;    # Mean reversion par in MS-model = not the lambda of noise fields  =#
    nfs = 0 # needs to have value for plotting purposes
    Ptrue = MarslandShardlow(a, c, γ, λ, n)
else
    db = 2.0 # domainbound
    nfstd = 1.0 # tau , width of noisefields
    nfs = construct_nfs(db, nfstd, γ) # 3rd argument gives average noise of positions (with superposition)
    Ptrue = Landmarks(a, c, n, db, nfstd, nfs)
end

if (model == :ahs) & showplotσq
    plotσq(db, nfs)
end

x0, xobs0, xobsT, Xf, Ptrue, pb, obs_atzero  = generatedata(dataset,Ptrue,t,σobs)

# step-size on mala steps for initial state
if obs_atzero
    δ = [0.0, 0.01] # in this case first comp is not used
else
    δ = [0.0005, 0.005]
end

(ainit, cinit, γinit) = (1.0, 1.0, 1.0) #(0.1, 0.1, 0.1) # (a,c,γ)
if model == :ms
    P = MarslandShardlow(ainit, cinit, γinit, Ptrue.λ, Ptrue.n)
elseif model == :ahs
    nfsinit = construct_nfs(Ptrue.db, Ptrue.nfstd, γinit)
    P = Landmarks(ainit, cinit, Ptrue.n, Ptrue.db, Ptrue.nfstd, nfsinit)
end

mT = zeros(PointF,P.n)   # vector of momenta at time T used for constructing guiding term
#mT = randn(PointF,P.n)

start = time() # to compute elapsed time

if obs_atzero     # only update mom
    xobsT = [xobsT]  # should be a vector
    xinit = State(xobs0, zeros(PointF,P.n))
    initstate_updatetypes = [:mala_mom]
elseif !obs_atzero & !fixinitmomentato0 # multiple shapes: update both pos and mom
    θ = π/6
    rot =  SMatrix{2,2}(cos(θ), sin(θ), -sin(θ), cos(θ))
    xinit = 1.3*State([rot * xobsT[1][i] for i in 1:P.n], zeros(PointF,P.n))
    #xinit = 1.3*State([rot * xobsTvec[1][i] for i in 1:n], randn(PointF,P.n))
    #xinit = State(xobsT[1], zeros(PointF, P.n))
    initstate_updatetypes = [:mala_mom, :rmmala_pos]
elseif !obs_atzero & fixinitmomentato0   # multiple shapes only update pos, so initialised momenta (equal toz zero) are maintained throughout the iterations
    xinit = 1.2*State(xobsT[1], zeros(PointF,P.n))
    θ = π/6
    rot =  SMatrix{2,2}(cos(θ), sin(θ), -sin(θ), cos(θ))
    xinit = 1.3*State([rot * xobsT[1][i] for i in 1:P.n], zeros(PointF,P.n))
    #xinit = State(xobsT[1], zeros(PointF, P.n))
    initstate_updatetypes = [:rmmala_pos]
end

anim, Xsave, parsave, objvals, acc_pcn, accinfo = lm_mcmc(tt_, (xobs0,xobsT), σobs, mT, P,
         sampler, obs_atzero, fixinitmomentato0,
         xinit, ITER, subsamples,
        (ρ, δ, prior_a, prior_c, prior_γ, σ_a, σ_c, σ_γ, η), initstate_updatetypes, adaptskip,
        outdir,  pb; updatepars = true, makefig=true, showmomenta=false)

elapsed = time() - start

println("Elapsed    time: ",round(elapsed/60;digits=2), " minutes")
println("Acceptance percentage pCN step: ",acc_pcn)

include("./postprocessing.jl")


text:
using TensorContractionsXD, Test

@testset "Optimal contraction order" begin
    _,cost = @optimalcontractiontree A[-1,1,2,3]*B[2,4,5,6]*C[1,5,7,-3]*D[3,8,4,9]*
                                    E[6,9,7,10]*F[-2,8,11,12]*G[10,11,12,-4]
    @test cost == TensorContractionsXD.Poly{:χ,Int}([0,0,0,0,0,0,2,2,2])

    _,cost = @optimalcontractiontree A[1,4,5,6,7,8]*B[2,9,10,11,12,13]*C[3,14,15,16,17,18]*
                                    D[-6,9,23,24,25,26]*E[-5,5,19,20,21,22]*F[8,14,27,28,29,30]*
                                    G[12,15,31,32,33,34]*H[22,25,28,31,35,36,37,38]*I[-4,20,23,35,39,40,41,42]*
                                    J[42,36,37,38,43,44,45,46]*K[41,24,44,26,47,48]*L[19,40,21,43,49,50]*
                                    M[27,45,29,30,51,52]*N[46,32,33,34,53,54]*O[-2,-3,39,49,47,55]*
                                    P[4,50,6,7,51,56]*Q[48,10,11,53,13,57]*R[52,54,16,17,18,58]*S[55,56,57,58,-1,1,2,3]
    @test cost == TensorContractionsXD.Poly{:χ,Int}([0,0,0,0,0,0,0,0,0,5,5,0,1,1,3,0,3])

    _,cost = @optimalcontractiontree a[64,72,73,19,22]*b[65,74,75,23,76]*c[66,77,20,79,28]*d[67,21,24,29,34]*
                                    e[68,25,78,35,80]*f[69,81,30,83,84]*g[70,31,36,85,86]*h[71,37,82,87,88]*
                                    i[-5,19,20,21,1,2,4,5]*j[22,23,24,25,3,26,6,27]*k[28,29,30,31,7,8,32,33]*
                                    l[34,35,36,37,9,38,39,40]*m[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18]*
                                    n[10,11,13,14,41,42,43,44]*o[12,26,15,27,47,48,49,50]*
                                    p[16,17,32,33,53,54,55,56]*q[18,38,39,40,60,61,62,63]*r[-2,-3,-4,41,89]*
                                    s[72,73,42,47,90]*t[74,75,48,76,45]*u[77,43,79,53,46]*v[44,49,54,60,51]*
                                    w[50,78,61,80,52]*x[81,55,83,84,57]*y[56,62,85,86,58]*z[63,82,87,88,59]*
                                    aa[89,90,45,46,51,52,57,58,59,-1,64,65,66,67,68,69,70,71]
    @test cost == TensorContractionsXD.Poly{:χ,Int}([0,0,0,0,0,0,0,4,4,0,0,0,1,1,1,0,1,0,0,0,3,0,3,2,0,2,4])
end


text:
using Documenter
using PointerArithmetic

makedocs(
    sitename = "PointerArithmetic.jl",
    format = Documenter.HTML(),
    modules = [PointerArithmetic]
)

deploydocs(
    repo = "github.com/nickolasrm/PointerArithmetic.jl.git"
)


text:
include("../src/SymFinel.jl")

using .SymFinelTests

include("ExprTests.jl")
include("InterpTests.jl")
include("DomainTests.jl")
include("FinelTests.jl")
include("ResidualTest.jl")

include("../examples/AxialBeamElement.jl")
include("../examples/READMEexample.jl")