text:
@testset "SparseRadialMapComponent (Nx,p) = (1,-1)" begin
    x = randn()

    C = SparseRadialMapComponent(1,-1)

    @test C.Nx == 1
    @test C.p == [-1]
    @test C.activedim == Int64[]
    @test size(C.ξ, 1)==1
    @test size(C.ξ[1], 1)==0

    @test size(C.σ, 1)==1
    @test size(C.σ[1], 1)==0

    @test size(C.coeff,1)==1
    @test size(C.coeff[1],1)==0

    @test isapprox(C(x), x, atol = 1e-10)
end

@testset "SparseRadialMapComponent (Nx,p) = (1,0)" begin

    C = SparseRadialMapComponent(1,0)
    x = randn()

    @test C.Nx == 1
    @test C.p == [0]
    @test C.activedim == Int64[1]
    @test size(C.ξ, 1)==1
    @test size(C.ξ[1], 1)==0

    @test size(C.σ, 1)==1
    @test size(C.σ[1], 1)==0

    @test size(C.coeff,1)==1
    @test size(C.coeff[1],1)==2

    C.coeff[1]=[2.0;3.0]
    @test isapprox(C(x), 2.0 + 3.0*x, atol = 1e-10)
end

@testset "SparseRadialMapComponent (Nx,p) = (1,1)" begin

    C = SparseRadialMapComponent(1,1)
    x = randn()

    @test C.Nx == 1
    @test C.p == [1]
    @test C.activedim == Int64[1]

    @test size(C.ξ,1) ==1
    @test size(C.ξ[1],1) ==1+2

    @test size(C.σ,1) ==1
    @test size(C.σ[1],1) ==1+2

    @test size(C.ξ,1) ==1
    @test size(C.ξ[1],1) ==1+2

    @test size(C.coeff,1) ==1
    @test size(C.coeff[1],1) ==1+3

    C.coeff[1] .= [1.0; 2.0; 3.0; -1.0]

    @test isapprox(C(x), 1.0 + 2.0*ψ₀(0.0,1.0)(x) +
                     3.0*ψj(0.0,1.0)(x) - 1.0*ψpp1(0.0,1.0)(x), atol = 1e-10)

    @test isapprox(C([x]), 1.0 + 2.0*ψ₀(0.0,1.0)(x) +
                      3.0*ψj(0.0,1.0)(x) - 1.0*ψpp1(0.0,1.0)(x), atol = 1e-10)
end

@testset "SparseRadialMapComponent (Nx,p) = (Nx,-1) with Nx>1" begin

    C = SparseRadialMapComponent(3, -1)
    x = randn()
    y = randn(3)

    @test C.Nx == 3
    @test C.p == [-1; -1 ; -1]
    @test C.activedim == Int64[]

    @test size(C.ξ,1)==3
    @test size(C.σ,1)==3
    @test size(C.coeff,1)==3

    @test size(C.ξ[1],1)==0
    @test size(C.σ[1],1)==0
    @test size(C.coeff[1],1)==0

    @test size(C.ξ[2],1)==0
    @test size(C.σ[2],1)==0
    @test size(C.coeff[2],1)==0

    @test size(C.ξ[3],1)==0
    @test size(C.σ[3],1)==0
    @test size(C.coeff[3],1)==0

    @test isapprox(C(x), x, atol = 1e-10)
    @test isapprox(C(y), y[3], atol = 1e-10)
end

@testset "SparseRadialMapComponent (Nx,p) = (Nx,0) with Nx>1" begin

    C = SparseRadialMapComponent(3, 0)
    x = randn()
    y = randn(3)

    @test C.Nx == 3
    @test C.p == [0 ; 0 ; 0]
    @test C.activedim == collect(1:3)

    @test size(C.ξ,1)==3
    @test size(C.σ,1)==3
    @test size(C.coeff,1)==3

    @test size(C.ξ[1],1)==0
    @test size(C.σ[1],1)==0
    @test size(C.coeff[1],1)==1

    @test size(C.ξ[2],1)==0
    @test size(C.σ[2],1)==0
    @test size(C.coeff[2],1)==1

    @test size(C.ξ[3],1)==0
    @test size(C.σ[3],1)==0
    @test size(C.coeff[3],1)==2

    @test isapprox(C(x), 0.0, atol = 1e-10)

    C.coeff[1][1] = 2.0

    @test isapprox(C([x;0.0;0.0]), 2.0*x, atol = 1e-10)

    C.coeff[2][1] = 3.0
    @test isapprox(C([y[1];y[2];0.0]), 2*y[1] + 3*y[2], atol = 1e-10)

    C.coeff[3][1] = 4.0
    @test isapprox(C(y), 2.0*y[1] + 3.0*y[2] + 4.0, atol = 1e-10)
    @test isapprox(C(x), 2.0*x+3.0*x+ 4.0, atol = 1e-10)

    C.coeff[3][2] = -1.0
    @test isapprox(C(y), 2.0*y[1] + 3.0*y[2]+ 4.0 - 1.0*y[3], atol = 1e-10)
    @test isapprox(C(x), 2.0*x + 3.0*x + 4.0 + -1.0*x, atol = 1e-10)

    clearcoeff!(C)
    @test norm(extractcoeff(C))<1e-10

end

@testset "SparseRadialMapComponent (Nx,p) = (Nx,p) with Nx>1 and p>1" begin

    C = SparseRadialMapComponent(3, 2)
    x = randn()
    y = randn(3)

    @test C.Nx == 3
    @test C.p == [2;2;2]
    @test C.activedim == collect(1:3)

    @test size(C.ξ,1)==3
    @test size(C.σ,1)==3
    @test size(C.coeff,1)==3

    @test size(C.ξ[1],1)==2
    @test size(C.σ[1],1)==2
    @test size(C.coeff[1],1)==2+1

    @test size(C.ξ[2],1)==2
    @test size(C.σ[2],1)==2
    @test size(C.coeff[2],1)==2+1

    @test size(C.ξ[3],1)==2+2
    @test size(C.σ[3],1)==2+2
    @test size(C.coeff[3],1)==2+3

    @test C.ξ[1]==zeros(2)
    @test C.ξ[2]==zeros(2)
    @test C.ξ[3]==zeros(2+2)

    @test C.σ[1]==ones(2)
    @test C.σ[2]==ones(2)
    @test C.σ[3]==ones(2+2)

    @test C.coeff[1]==zeros(2+1)
    @test C.coeff[2]==zeros(2+1)
    @test C.coeff[3]==zeros(2+3)

    C.coeff[1] .= ones(3)

    @test isapprox(C(x), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(x), atol = 1e-10)

    @test isapprox(C([x; 0.0; 0.0]), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(x), atol = 1e-10)

    C.coeff[2] .= 2*ones(3)

    @test isapprox(C([y[1]; y[2];  0.0]), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(y[1]) + ui(2,C.ξ[2], C.σ[2], C.coeff[2])(y[2]), atol = 1e-10)

    C.coeff[3] .= 3*ones(5)

    @test isapprox(C(y), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(y[1]) +
                                 ui(2,C.ξ[2], C.σ[2], C.coeff[2])(y[2]) +
                                 uk(2,C.ξ[3], C.σ[3], C.coeff[3])(y[3]), atol = 1e-10)

    @test isapprox(C(x), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(x) +
                                 ui(2,C.ξ[2], C.σ[2], C.coeff[2])(x) +
                                 uk(2,C.ξ[3], C.σ[3], C.coeff[3])(x), atol = 1e-10)

     clearcoeff!(C)
     @test norm(extractcoeff(C))<1e-10
end

@testset "SparseRadialMapComponent (Nx,p) = (Nx,p) with Nx>1 and p =[-1  0 -1]" begin

    C = SparseRadialMapComponent(3, [-1; 0; -1])
    x = randn()
    y = randn(3)

    @test C.Nx == 3
    @test C.p == [-1; 0; -1]
    @test C.activedim == Int64[2]

    @test size(C.ξ,1)==3
    @test size(C.σ,1)==3
    @test size(C.coeff,1)==3

    @test size(C.ξ[1],1)==0
    @test size(C.σ[1],1)==0
    @test size(C.coeff[1],1)==0

    @test size(C.ξ[2],1)==0
    @test size(C.σ[2],1)==0
    @test size(C.coeff[2],1)==1

    @test size(C.ξ[3],1)==0
    @test size(C.σ[3],1)==0
    @test size(C.coeff[3],1)==0

    @test C.ξ[1]==zeros(0)
    @test C.ξ[2]==zeros(0)
    @test C.ξ[3]==zeros(0)

    @test C.σ[1]==ones(0)
    @test C.σ[2]==ones(0)
    @test C.σ[3]==ones(0)

    @test C.coeff[1]==zeros(0)
    @test C.coeff[2]==zeros(1)
    @test C.coeff[3]==zeros(0)

    C.coeff[2] .= [randn()]

    @test isapprox(C(x), ui(0,C.ξ[2], C.σ[2], C.coeff[2])(x)+uk(-1,C.ξ[3], C.σ[3], C.coeff[3])(x), atol = 1e-10)
    @test isapprox(C(y), ui(-1,C.ξ[1], C.σ[1], C.coeff[1])(y[1])+
                         ui(0,C.ξ[2], C.σ[2], C.coeff[2])(y[2])+
                         uk(-1,C.ξ[3], C.σ[3], C.coeff[3])(y[3]), atol = 1e-10)

     clearcoeff!(C)
     @test norm(extractcoeff(C))<1e-10

end

@testset "SparseRadialMapComponent (Nx,p) = (Nx,p) with Nx>1 and p =[1  0 -1]" begin

    C = SparseRadialMapComponent(3, [1; 0; -1])
    x = randn()
    y = randn(3)

    @test C.Nx == 3
    @test C.p == [1; 0; -1]
    @test C.activedim == [1; 2]

    @test size(C.ξ,1)==3
    @test size(C.σ,1)==3
    @test size(C.coeff,1)==3

    @test size(C.ξ[1],1)==1
    @test size(C.σ[1],1)==1
    @test size(C.coeff[1],1)==2

    @test size(C.ξ[2],1)==0
    @test size(C.σ[2],1)==0
    @test size(C.coeff[2],1)==1

    @test size(C.ξ[3],1)==0
    @test size(C.σ[3],1)==0
    @test size(C.coeff[3],1)==0

    @test C.ξ[1]==zeros(1)
    @test C.ξ[2]==zeros(0)
    @test C.ξ[3]==zeros(0)

    @test C.σ[1]==ones(1)
    @test C.σ[2]==ones(0)
    @test C.σ[3]==ones(0)

    @test C.coeff[1]==zeros(2)
    @test C.coeff[2]==zeros(1)
    @test C.coeff[3]==zeros(0)

    C.coeff[1] .= randn(2)
    @test isapprox(C(x), ui(1,C.ξ[1], C.σ[1], C.coeff[1])(x)+
                         ui(0,C.ξ[2], C.σ[2], C.coeff[2])(x)+
                         uk(-1,C.ξ[3], C.σ[3], C.coeff[3])(x), atol = 1e-10)
    C.coeff[2] .= [randn()]
    @test isapprox(C(x), ui(1,C.ξ[1], C.σ[1], C.coeff[1])(x) +
                         ui(0,C.ξ[2], C.σ[2], C.coeff[2])(x) + x, atol = 1e-10)

    @test isapprox(C(y), ui(1,C.ξ[1], C.σ[1], C.coeff[1])(y[1]) +
                         ui(0,C.ξ[2], C.σ[2], C.coeff[2])(y[2]) + y[3], atol = 1e-10)

     clearcoeff!(C)
     @test norm(extractcoeff(C))<1e-10

end

@testset "SparseRadialMapComponent (Nx,p) = (Nx,p) with Nx>1 and p =[-1  0 0]" begin

    C = SparseRadialMapComponent(3, [-1; 0; 0])
    x = randn()
    y = randn(3)

    @test C.Nx == 3
    @test C.p == [-1; 0; 0]
    @test C.activedim == [2; 3]

    @test size(C.ξ,1)==3
    @test size(C.σ,1)==3
    @test size(C.coeff,1)==3

    @test size(C.ξ[1],1)==0
    @test size(C.σ[1],1)==0
    @test size(C.coeff[1],1)==0

    @test size(C.ξ[2],1)==0
    @test size(C.σ[2],1)==0
    @test size(C.coeff[2],1)==1

    @test size(C.ξ[3],1)==0
    @test size(C.σ[3],1)==0
    @test size(C.coeff[3],1)==2

    @test C.ξ[1]==zeros(0)
    @test C.ξ[2]==zeros(0)
    @test C.ξ[3]==zeros(0)

    @test C.σ[1]==ones(0)
    @test C.σ[2]==ones(0)
    @test C.σ[3]==ones(0)

    @test C.coeff[1]==zeros(0)
    @test C.coeff[2]==zeros(1)
    @test C.coeff[3]==zeros(2)

    C.coeff[2] .= [randn()]
    @test isapprox(C(x), ui(0,C.ξ[2], C.σ[2], C.coeff[2])(x), atol = 1e-10)
    @test isapprox(C(y), ui(0,C.ξ[2], C.σ[2], C.coeff[2])(y[2]), atol = 1e-10)

    C.coeff[3] .= randn(2)
    @test isapprox(C(x), ui(0,C.ξ[2], C.σ[2], C.coeff[2])(x) +
                         uk(0,C.ξ[3], C.σ[3], C.coeff[3])(x), atol = 1e-10)

    @test isapprox(C(y), ui(0,C.ξ[2], C.σ[2], C.coeff[2])(y[2]) +
                         uk(0,C.ξ[3], C.σ[3], C.coeff[3])(y[3]), atol = 1e-10)

     clearcoeff!(C)
     @test norm(extractcoeff(C))<1e-10
end

@testset "SparseRadialMapComponent (Nx,p) = (Nx,p) with Nx>1 and p =[2  -1  1]" begin

    C = SparseRadialMapComponent(3, [2; -1; 1])
    x = randn()
    y = randn(3)

    @test C.Nx == 3
    @test C.p == [2; -1; 1]
    @test C.activedim == [1; 3]

    @test size(C.ξ,1)==3
    @test size(C.σ,1)==3
    @test size(C.coeff,1)==3

    @test size(C.ξ[1],1)==2
    @test size(C.σ[1],1)==2
    @test size(C.coeff[1],1)==3

    @test size(C.ξ[2],1)==0
    @test size(C.σ[2],1)==0
    @test size(C.coeff[2],1)==0

    @test size(C.ξ[3],1)==3
    @test size(C.σ[3],1)==3
    @test size(C.coeff[3],1)==4

    @test C.ξ[1]==zeros(2)
    @test C.ξ[2]==zeros(0)
    @test C.ξ[3]==zeros(3)

    @test C.σ[1]==ones(2)
    @test C.σ[2]==ones(0)
    @test C.σ[3]==ones(3)

    @test C.coeff[1]==zeros(3)
    @test C.coeff[2]==zeros(0)
    @test C.coeff[3]==zeros(4)

    @test isapprox(C(x), 0.0, atol =  1e-10)
    @test isapprox(C(y), 0.0, atol = 1e-10)

    C.coeff[1] .= randn(3)
    @test isapprox(C(x), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(x), atol  = 1e-10)
    @test isapprox(C(y), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(y[1]), atol = 1e-10)

    C.coeff[3] .= randn(4)
    @test isapprox(C(x), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(x)+uk(1,C.ξ[3], C.σ[3], C.coeff[3])(x), atol = 1e-10)

    @test isapprox(C(y), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(y[1]) + uk(1,C.ξ[3], C.σ[3], C.coeff[3])(y[3]), atol = 1e-10)

     clearcoeff!(C)
     @test norm(extractcoeff(C))<1e-10

end

@testset "SparseRadialMapComponent (Nx,p) = (Nx,p) with Nx>1 and p =[2  -1  2]" begin

    C = SparseRadialMapComponent(3, [2; -1; 2])
    x = randn()
    y = randn(3)

    @test C.Nx == 3
    @test C.p == [2; -1; 2]
    @test C.activedim == [1; 3]

    @test size(C.ξ,1)==3
    @test size(C.σ,1)==3
    @test size(C.coeff,1)==3

    @test size(C.ξ[1],1)==2
    @test size(C.σ[1],1)==2
    @test size(C.coeff[1],1)==3

    @test size(C.ξ[2],1)==0
    @test size(C.σ[2],1)==0
    @test size(C.coeff[2],1)==0

    @test size(C.ξ[3],1)==4
    @test size(C.σ[3],1)==4
    @test size(C.coeff[3],1)==5

    @test C.ξ[1]==zeros(2)
    @test C.ξ[2]==zeros(0)
    @test C.ξ[3]==zeros(4)

    @test C.σ[1]==ones(2)
    @test C.σ[2]==ones(0)
    @test C.σ[3]==ones(4)

    @test C.coeff[1]==zeros(3)
    @test C.coeff[2]==zeros(0)
    @test C.coeff[3]==zeros(5)

    @test isapprox(C(x), 0.0, atol = 1e-10)
    @test isapprox(C(y), 0.0, atol = 1e-10)

    C.coeff[1] .= randn(3)
    @test isapprox(C(x), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(x), atol = 1e-10)
    @test isapprox(C(y), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(y[1]), atol = 1e-10)

    C.coeff[3] .= randn(5)
    @test isapprox(C(x), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(x) + uk(2,C.ξ[3], C.σ[3], C.coeff[3])(x), atol = 1e-10)
    @test isapprox(C(y), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(y[1]) + uk(2,C.ξ[3], C.σ[3], C.coeff[3])(y[3]), atol = 1e-10)

     clearcoeff!(C)
     @test norm(extractcoeff(C))<1e-10
end

@testset "SparseRadialMapComponent (Nx,p) = (Nx,p) with Nx>1 and p =[2  2  -1]" begin

    C = SparseRadialMapComponent(3, [2; 2; -1])
    x = randn()
    y = randn(3)

    @test C.Nx == 3
    @test C.p == [2; 2; -1]
    @test C.activedim == [1; 2]

    @test size(C.ξ,1)==3
    @test size(C.σ,1)==3
    @test size(C.coeff,1)==3

    @test size(C.ξ[1],1)==2
    @test size(C.σ[1],1)==2
    @test size(C.coeff[1],1)==3

    @test size(C.ξ[2],1)==2
    @test size(C.σ[2],1)==2
    @test size(C.coeff[2],1)==3

    @test size(C.ξ[3],1)==0
    @test size(C.σ[3],1)==0
    @test size(C.coeff[3],1)==0

    @test C.ξ[1]==zeros(2)
    @test C.ξ[2]==zeros(2)
    @test C.ξ[3]==zeros(0)

    @test C.σ[1]==ones(2)
    @test C.σ[2]==ones(2)
    @test C.σ[3]==ones(0)

    @test C.coeff[1]==zeros(3)
    @test C.coeff[2]==zeros(3)
    @test C.coeff[3]==zeros(0)

    @test isapprox(C(x), x, atol = 1e-10)
    @test isapprox(C(y), y[3], atol = 1e-10)

    C.coeff[1] .= randn(3)
    @test isapprox(C(x), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(x) + x, atol = 1e-10)
    @test isapprox(C(y), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(y[1])+ y[3], atol = 1e-10)

    C.coeff[2] .= randn(3)
    @test isapprox(C(x), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(x) +
                         ui(2,C.ξ[2], C.σ[2], C.coeff[2])(x) +
                         uk(-1,C.ξ[3], C.σ[3], C.coeff[3])(x), atol = 1e-10)

    @test isapprox(C(y), ui(2,C.ξ[1], C.σ[1], C.coeff[1])(y[1])+ ui(2,C.ξ[2], C.σ[2], C.coeff[2])(y[2]) +
                         y[3], atol = 1e-10)

    clearcoeff!(C)
    @test norm(extractcoeff(C))<1e-10
end

@testset "set_id function" begin
    C = SparseRadialMapComponent(1, 1)
    set_id(C)
    @test C.Nx==1
    @test C.ξ[1]== Float64[]
    @test C.σ[1]== Float64[]
    @test C.coeff[1]== Float64[]
    @test C.p == [-1]

    C = SparseRadialMapComponent(4, [1;2;-1;-1])
    set_id(C)
    @test C.Nx ==4
    for i=1:4
        @test C.ξ[i]== Float64[]
        @test C.σ[i]== Float64[]
        @test C.coeff[i]== Float64[]
        @test C.p[i]== -1
    end
end

@testset "Component function for SparseRadialMapComponent type" begin

    C = SparseRadialMapComponent(1, 2)

    @test C.Nx == 1
    @test C.p == [2]

    @test  TransportBasedInference.component(C, 1).ξk == zeros(2+2)
    @test  TransportBasedInference.component(C, 1).σk == ones(2+2)
    @test  TransportBasedInference.component(C, 1).coeffk == zeros(2+3)

    C = SparseRadialMapComponent(5, 2)

    @test C.Nx == 5
    @test C.p == [2; 2; 2; 2; 2]

    @test  TransportBasedInference.component(C, 5).ξk == zeros(2+2)
    @test  TransportBasedInference.component(C, 5).σk == ones(2+2)
    @test  TransportBasedInference.component(C, 5).coeffk == zeros(2+3)

    @test  TransportBasedInference.component(C, 2).ξi == zeros(2)
    @test  TransportBasedInference.component(C, 2).σi == ones(2)
    @test  TransportBasedInference.component(C, 2).coeffi == zeros(2+1)

    C = SparseRadialMapComponent(5, [1; 0; -1; 2; 1])

    solξ = Int64[1; 0 ; 0; 2; 3]
    solσ = Int64[1; 0 ; 0; 2; 3]
    sola = Int64[2; 1 ; 0; 3; 4]

    for i=1:4
        @test TransportBasedInference.component(C, i).p == C.p[i]
        @test typeof(TransportBasedInference.component(C, i))<:ui
        @test  TransportBasedInference.component(C, i).ξi == zeros(solξ[i])
        @test  TransportBasedInference.component(C, i).σi == ones(solσ[i])
        @test  TransportBasedInference.component(C, i).coeffi == zeros(sola[i])
    end

    for i=5:5
        @test TransportBasedInference.component(C, i).p == C.p[i]
        @test typeof(TransportBasedInference.component(C, i))<:uk
        @test  TransportBasedInference.component(C, i).ξk == zeros(solξ[i])
        @test  TransportBasedInference.component(C, i).σk == ones(solσ[i])
        @test  TransportBasedInference.component(C, i).coeffk == zeros(sola[i])
    end

end

@testset "Verify off_diagonal function" begin

    # Nx=1  & p=0
    C = SparseRadialMapComponent(1,0)
    C.coeff[1] .= randn(2)

    @test off_diagonal(C, randn()) == 0.0

    # Nx=1 & p= 3
    C = SparseRadialMapComponent(1, 3)
    C.coeff[1] .= rand(6)

    @test off_diagonal(C, randn()) == 0.0

    # Nx=3 & p = 0
    C = SparseRadialMapComponent(3, 0)
    a1 =randn()
    a2 = randn()
    a3 = rand(2)
    C.coeff[1] .= a1
    C.coeff[2] .= a2
    C.coeff[3] .= a3

    z = randn(3)
    @test norm(off_diagonal(C, z) - (C(z) - TransportBasedInference.component(C,3)(z[3])))<1e-10

    # Nx=3 & p = 3
    C = SparseRadialMapComponent(3, 3)
    for i=1:2
    C.ξ[i] .= randn(3)
    C.σ[i] .= rand(3)
    end

    C.ξ[3] .= randn(5)
    C.σ[3] .= rand(5)

    a1 =randn(4)
    a2 = randn(4)
    a3 = rand(6)
    C.coeff[1] .= a1
    C.coeff[2] .= a2
    C.coeff[3] .= a3

    z = randn(3)
    @test norm(off_diagonal(C, z) - (C(z) - TransportBasedInference.component(C,3)(z[3])))<1e-10

    # Nx=1  & p=-1
    C = SparseRadialMapComponent(1,-1)

    @test off_diagonal(C, 2.0)==0.0

    # Nx=3  & p=[-1 -1 -1]
    C = SparseRadialMapComponent(3,-1)

    @test off_diagonal(C, 2.0) == 0.0

    # Nx=3  & p=[0 -1 -1]
    C = SparseRadialMapComponent(3, [0; -1; -1])
    a1 = randn(1)
    C.coeff[1] .= a1

    z = randn(3)
    @test norm(off_diagonal(C, z) - (C(z) - TransportBasedInference.component(C,3)(z[3])))<1e-10

    # Nx=3  & p=[-1  2 -1]
    C = SparseRadialMapComponent(3, [-1; 2; -1])
    ξ2 = randn(2)
    σ2 = rand(2)
    a2 = randn(3)

    C.ξ[2] .= ξ2
    C.σ[2] .= σ2
    C.coeff[2] .= a2

    z = randn(3)
    @test norm(off_diagonal(C, z) - (C(z) - TransportBasedInference.component(C,3)(z[3])))<1e-10
end

@testset "extract and modify coefficients of RadialMapComponent" begin
    # Nx=1 and p=0
    C = SparseRadialMapComponent(1, 0)

    @test C.Nx == 1
    @test C.p == [0]

    modifycoeff!(C, [1.0;2.0])
    @test extractcoeff(C) == [1.0; 2.0]

    C = RadialMapComponent(1, 3)
    modifycoeff!(C, [1.0;2.0; 3.0; 4.0; 5.0; 6.0])
    @test extractcoeff(C) == [1.0;2.0; 3.0; 4.0; 5.0; 6.0]

    # Nx =3 and p=2
    A = collect(1.0:1.0:11.0)
    C = SparseRadialMapComponent(3, 2)

    modifycoeff!(C, A)
    @test  extractcoeff(C) == A

    C = RadialMapComponent(3, 2)
    C.coeff[1] = [1.0; 2.0; 3.0]
    C.coeff[2] = [4.0; 5.0; 6.0]
    C.coeff[3] = [7.0; 8.0; 9.0; 10.0; 11.0]

    @test extractcoeff(C) ==A

    # Nx = 3 and p = 0

    C = SparseRadialMapComponent(3, 0)

    @test C.Nx == 3
    @test C.p == [0; 0; 0]

    C.coeff[1] .= [1.0]
    C.coeff[2] .= [2.0]
    C.coeff[3] .= [3.0; 4.0]

    @test extractcoeff(C)==collect(1.0:4.0)

    modifycoeff!(C, collect(5.0:8.0))

    @test extractcoeff(C)==collect(5.0:8.0)

    # Nx = 5 and p = [-1 2 -1 0 -1]
    C = SparseRadialMapComponent(5, [-1; 2; -1; 0; -1])
    A = randn(4)
    modifycoeff!(C, A)
    @test extractcoeff(C)==A

    # Nx = 11 and p = [2 -1 0 -1 0 2 -1 0 -1 2  1]
    C = SparseRadialMapComponent(11, [2; -1; 0; -1; 0; 2; -1; 0; -1; 2; 2])
    A = randn(17)
    modifycoeff!(C, A)
    @test C.coeff[1]  == A[1:3]
    @test C.coeff[2]  == Float64[]
    @test C.coeff[3]  == [A[4]]
    @test C.coeff[4]  == Float64[]
    @test C.coeff[5]  == [A[5]]
    @test C.coeff[6]  == A[6:8]
    @test C.coeff[7]  == Float64[]
    @test C.coeff[8]  == [A[9]]
    @test C.coeff[9]  == Float64[]
    @test C.coeff[10] == A[10:12]
    @test C.coeff[11] == A[13:17]

    @test extractcoeff(C)==A

    # There is an extra verification that we have made exactly as much assignment
    # that there are component of C.coeff
end


text:
using Test
using KnetMetrics
using Random

@testset "regression-metrics" begin

    Random.seed!(42)
    u = rand(100)
    Random.seed!(123)
    v = rand(100)

    @testset "max-error" begin
        @test isapprox(max_error(u,v), 0.8624100780860133)
    end

    @testset "max-error" begin
        @test isapprox(max_error(u,v), 0.8624100780860133)
    end

    @testset "mean-absolute-error" begin
        @test isapprox(mean_absolute_error(u,v), 0.3666602951971807)
    end

    @testset "mean-squared-error" begin
        @test isapprox(mean_squared_error(u,v), 0.3666602951971807)
    end

    @testset "mean-squared-log-error" begin
        @test isapprox(mean_squared_log_error(u,v), 0.24971305530593674)
    end

    @testset "median-absolute-error" begin
        @test isapprox(median_absolute_error(u,v), 0.3183007014762824)
    end

    @testset "mean-absolute-percentage-error" begin
        @test isapprox(mean_absolute_percentage_error(u,v), 5.329695963484449)
    end

end

ext: jl


text:
using Combinatorics

@testset "Testing Pure Nash Equilibrium Routines" begin

    # Pure Action Nash equilibrium
    @testset "2x2 Game with 1 Pure Action Nash equilibrium" begin
        A = [9.0 1.0
             10.0 3.0]

        nfg = NormalFormGame(A)
        ne = pure_nash(nfg)

        @test ne == [(2, 2)]
    end

    @testset "Matching Pennies Game with 0 Pure Action Nash equilibrium" begin
        MP = [1.0 -1.0
              -1.0 1.0]
        p1 = Player(MP)
        p2 = Player(-MP)

        g_MP = NormalFormGame(p1, p2)
        ne = pure_nash(g_MP)

        @test ne == []
    end

    @testset "Coordination game with 2 Pure Action Nash equilibria" begin
        Coo = [4.0 0.0
               3.0 2.0]

        g_Coo = NormalFormGame(Coo)
        ne = pure_nash(g_Coo)

        @test sort(ne) == sort([(1,1); (2,2)])
    end

    @testset "Coordination game with 2 Pure Action Nash equilibria but only find 1" begin
        Coo = [4.0 0.0
               3.0 2.0]

        g_Coo = NormalFormGame(Coo)
        ne = pure_nash(g_Coo; ntofind=1)

        @test length(ne) == 1
    end

    @testset "Unanimity Game with more than two players" begin
        N = 4
        a, b = 1, 2
        g_Unanimity = NormalFormGame(tuple(fill(2, N)...))
        g_Unanimity[fill(1, N)...] = fill(a, N)
        g_Unanimity[fill(2, N)...] = fill(b, N)

        Unanimity_NE = [tuple(fill(1, N)...)]
        for k in 2:N-2
            for ind in combinations(1:N, k)
                a = fill(1, N)
                a[ind] .= 2
                push!(Unanimity_NE, tuple(a...))
            end
        end
        push!(Unanimity_NE, tuple(fill(2, N)...))

        ne = pure_nash(g_Unanimity)

        @test sort(ne) == sort(Unanimity_NE)
    end

    @testset "Tolerance" begin
        epsilon = 1e-08

        g = NormalFormGame((2, 2))
        g[1, 1] = [1, 1]
        g[1, 2] = [-2, 1 + epsilon]
        g[2, 1] = [1 + epsilon, -2]
        g[2, 2] = [0, 0];

        NEs = [[(2, 2)]]
        epsilon_NEs = [[(1, 1); (2, 2)]]

        for (tol, answer) in zip([0 epsilon], [NEs epsilon_NEs])
            @test sort(pure_nash(g, tol=tol)) ==  sort(answer)
        end
    end

    @testset "Trivial game with 1 player" begin
        n = 3
        g1 = NormalFormGame(Player(collect(1:n)))
        @test pure_nash(g1) == [(n,)]
        @test sort(pure_nash(g1, tol=1.)) == [(n-1,), (n,)]
    end

end


text:
function _ladderise!(p::AbstractNode, q::AbstractNode, rev::Bool)
	istip(q) && return 1
	stsizes = zeros(Int, length(q.links))
	@inbounds for (i, link) in enumerate(q.links)
		link.to == p && continue
		stsizes[i] = _ladderise!(q, link.to, rev)
	end
	permute!(q.links, sortperm(stsizes; rev=rev))

	return sum(stsizes) + 1
end

"""
	ladderise!(tree [, sense = :right])

Rearrange the links in each node of `tree` so that bigger subtrees are always to the right (`sense = :right`) or the left (`sense = :left`).

Left and right orientations are relative to the origin node of the tree.
"""
function ladderise!(tree::AbstractTree; sense=:right)
	if sense == :right
		_ladderise!(tree.origin, tree.origin, false)
	elseif sense == :left
		_ladderise!(tree.origin, tree.origin, true)
	else
		throw(
			ErrorException("invalid ladderisation sense " * ":" * string(sense))
			)
	end

	return nothing
end


text:
#!/usr/bin/env julia

using FITSIO
using Celeste
using JLD
using DataFrames

function optimized_sources_to_fits(inpath, outdir)
    optimized_sources = load(inpath)["results"]
    sources_df = Celeste.AccuracyBenchmark.celeste_to_df(optimized_sources)

    # stores the optimized sources in an "attribute-major" format
    # suitable for writing to a fits table
    sources_dict = Dict{String, Vector}()
    sources_dict["thingid"] = [os.thingid for os in optimized_sources]
    for col_symbol in names(sources_df)
        if !any(isna, sources_df[:, col_symbol]) && col_symbol != :objid
            sources_dict[string(col_symbol)] = sources_df[col_symbol]
        end
    end

    outpath = joinpath(outdir, "$(basename(inpath)[1:end-4]).fits")
    outfits = FITS(outpath, "w")
    write(outfits, sources_dict)
    println("wrote $(length(optimized_sources)) to $outpath")
end

if length(ARGS) != 2 || ARGS[1][end-3:end] != ".jld"
    println("usage: optimized_sources_to_fits.jl <celeste-*.jld> <outdir>")
else
    optimized_sources_to_fits(ARGS[1], ARGS[2])
end


text:
"""
    function f_TCRDD(args...; kwargs...)

Calculates the Transmission Constrained Residual Demand Derivative (TCRDD) of the slack
bus of an optimised power system (optimised with OPF minimising total generation cost).

Detailed description: The TCRDD is calculated using the PTDF matrix of the system. This
matrix is divided in accordance with the generator segments (specified below) and load
buses as well as the binding lines of the solved system [1]. To determine if a generator
or a branch are binding, the dual variables from the OPF are compared against a tolerance.
The tolerance for binding lines and generators can be specified individually.

Generator Segments:
    1) Slack Generator
    2) Generators with binding active power
    3) Generators with non zero slopes in their cost function
    4) Generators with only constant price in their cost function

[1] L. Xu, R. Baldick and Y. Sutjandra, "Bidding Into Electricity Markets:
    A Transmission-Constrained Residual Demand Derivative Approach," in IEEE
    Transactions on Power Systems, vol. 26, no. 3, pp. 1380-1388, Aug. 2011,
    doi: 10.1109/TPWRS.2010.2083702.

#Arguments
- `sys::System`:                            Power system in p.u. (from PowerSystems.jl)
- `res::PowerSimulations.
    OperationsProblemResults`:              Results of the solved OPF for the system
                                             (from PowerSimulations.jl)

#Keywords
- `dual_lines_tol::Float64 = 1e-1`:         Tolerance to identify any binding lines
- `dual_gen_tol::Float64 = 1e-1`:           Tolerance to identify any binding generators

"""
function f_TCRDD(
    sys::System,
    res::PowerSimulations.OperationsProblemResults;
    dual_lines_tol::Float64 = 1e-1,
    dual_gen_tol::Float64 = 1e-1
    )
    # PTDF Matrices
    # Calculate PTDF matrix of the system
    PTDF_matrix = PTDF(sys)

    # PTDF for binding Lines
    PTDF_binding_lines = get_PTDF_bindingLines(sys, res, PTDF_matrix; dual_lines_tol)
    nbind_lines = length(PTDF_binding_lines[:,1]) #number binding lines

    # PTDF Generators separation
    # PTDF Slack is zeros
    gens_thermal_slack_busnumber, PTDFb_slack = get_PTDF_thermal_slack(
        sys,
        PTDF_binding_lines
    )
    ngens_thermal_slack = length(gens_thermal_slack_busnumber)
    # PTDF Pg Binding Gens (Checked with Duals)
    gens_thermal_bindPg_busnumber, PTDFb_bindingPg = get_PTDF_thermal_bindingPg(
        sys,
        res,
        PTDF_binding_lines;
        dual_gen_tol
    )
    ngens_thermal_bindPg = length(gens_thermal_bindPg_busnumber)
    # PTDF Pmax Non Zero Slope
    nonzeroslope_nb_PTDF = get_PTDF_thermal_nonzeroslope(
        sys,
        res,
        PTDF_binding_lines;
        dual_gen_tol
    )
    gens_thermal_nonzeroslope_busnumber, PTDFb_nonzeroslope = nonzeroslope_nb_PTDF

    ngens_thermal_nonzeroslope = length(gens_thermal_nonzeroslope_busnumber)
    # PTDF Pmax Non Constant Price
    gens_thermal_constprice_busnumber,PTDFb_constprice = get_PTDF_thermal_constprice(
        sys,
        res,
        PTDF_binding_lines;
        dual_gen_tol
    )
    ngens_thermal_constprice = length(gens_thermal_constprice_busnumber)

    # Load PTDF separation
    PQ_buses, PTDF_PQLoad = get_PTDF_load(sys, PTDF_binding_lines)

    # Obtaining vector of ones of size Slack, Pmax Pg, nonzeroslopes, constant price
    ones_slack = ones(ngens_thermal_slack)
    ones_bindPg = ones(ngens_thermal_bindPg)
    ones_nonzeroslope = ones(ngens_thermal_nonzeroslope)
    ones_constprice = ones(ngens_thermal_constprice)

    # Build zeros matrices just for Ax=b
    zeros_nonzeroslope_constprice = zeros(
        ngens_thermal_nonzeroslope,
        ngens_thermal_constprice
    )
    zeros_constprice_nonzeroslope = zeros_nonzeroslope_constprice'
    zeros_constprice_constprice = zeros(ngens_thermal_constprice,ngens_thermal_constprice)
    zeros_lines_lines = zeros(nbind_lines,nbind_lines)
    zeros_lines = zeros(nbind_lines,1)

    # Calculation of the second derivative of the objective function w.r.t
    # to the active power of the non zero slope generators
    d2f_PgPgnonzeroslope, d2f_PgPgnonzeroslope_MatrixInfo = d2f_PgPg_nonzeroslope(
        sys,
        res;
        dual_gen_tol
    )

    # Build the system of equations depending on the case
    if nbind_lines ≠ 0  && ngens_thermal_constprice ≠ 0
        #build A x = b
        A = [d2f_PgPgnonzeroslope zeros_nonzeroslope_constprice (-(PTDFb_nonzeroslope'));
            zeros_constprice_nonzeroslope zeros_constprice_constprice -(PTDFb_constprice');
            PTDFb_nonzeroslope PTDFb_constprice zeros_lines_lines]
        b = [ones_nonzeroslope; ones_constprice; zeros_lines]
        #x = [dPgnonzeroslope_dlam dPgconstprice_dlam dmulines_dlam]

        #Solve the system using JuMP
        model_a = Model(GLPK.Optimizer)
        @variable(model_a, dPgnonzeroslope_dlam_a[1:ngens_thermal_nonzeroslope])
        @variable(model_a, dPgconstprice_dlam_a[1:ngens_thermal_constprice])
        @variable(model_a, dmubind_dlam_a[1:nbind_lines])

        @objective(model_a, Min, 5) #Minimise any objective

        @constraint(model_a,const_a,
            A * [dPgnonzeroslope_dlam_a; dPgconstprice_dlam_a; dmubind_dlam_a] .== b)

        print(model_a)
        optimize!(model_a)

        println("Termination status : ", termination_status(model_a))
        println("Primal status      : ", primal_status(model_a))

        #obj_value = objective_value(model)
        dPgnonzeroslope_dlam_a_value = value.(dPgnonzeroslope_dlam_a)
        dPgconstprice_dlam_a_value = value.(dPgconstprice_dlam_a)
        dmubind_dlam_a_value = value.(dmubind_dlam_a)

        #Calculate the TCRDD
        tcrdd_slack = -((ones_nonzeroslope') * dPgnonzeroslope_dlam_a_value) -
            ((ones_constprice') * dPgconstprice_dlam_a_value)

    elseif nbind_lines ≠ 0  && ngens_thermal_constprice == 0
        #build A x = b
        A = [d2f_PgPgnonzeroslope (-(PTDFb_nonzeroslope'));
            PTDFb_nonzeroslope zeros_lines_lines]
        b = [ones_nonzeroslope; zeros_lines]
        #x = [dPgnonzeroslope_dlam dmulines_dlam]

        #Solve the system using JuMP
        model_b = Model(GLPK.Optimizer)
        @variable(model_b, dPgnonzeroslope_dlam_b[1:ngens_thermal_nonzeroslope])
        @variable(model_b, dmubind_dlam_b[1:nbind_lines])

        @objective(model_b, Min, 5) #Minimise any objective

        @constraint(model_b, const_b, A * [dPgnonzeroslope_dlam_b ; dmubind_dlam_b] .== b)

        print(model_b)
        optimize!(model_b)

        println("Termination status : ", termination_status(model_b))
        println("Primal status      : ", primal_status(model_b))

        #obj_value = objective_value(model)
        dPgnonzeroslope_dlam_b_value = value.(dPgnonzeroslope_dlam_b)
        dmubind_dlam_b_value = value.(dmubind_dlam_b)

        #Calculate the TCRDD
        tcrdd_slack = -(ones_nonzeroslope') * dPgnonzeroslope_dlam_b_value

    elseif nbind_lines == 0 #&& ngens_thermal_constprice ≠ 0
        A = d2f_PgPgnonzeroslope
        b = ones_nonzeroslope
        #x = [dPgnonzeroslope_dlam];

        #Solve the system using JuMP
        model_d = Model(GLPK.Optimizer)
        @variable(model_d, dPgnonzeroslope_dlam_d[1:ngens_thermal_nonzeroslope])

        @objective(model_d, Min, 5) #Minimise any objective

        @constraint(model_d, const_d, A * dPgnonzeroslope_dlam_d .== b)

        print(model_d)
        optimize!(model_d)

        println("Termination status : ", termination_status(model_d))
        println("Primal status      : ", primal_status(model_d))

        #obj_value = objective_value(model)
        dPgnonzeroslope_dlam_d_value = value.(dPgnonzeroslope_dlam_d)

        #Calculate the TCRDD
        tcrdd_slack = -(ones_nonzeroslope') * dPgnonzeroslope_dlam_d_value

    end
    return tcrdd_slack
end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
module NonconvexMultistart

export  HyperoptAlg,
        HyperoptOptions,
        RandomSampler,
        Hyperband,
        LHSampler,
        CLHSampler,
        GPSampler,
        DeflatedAlg,
        DeflatedOptions

import Hyperopt
using Reexport, Parameters, Setfield, ChainRulesCore, Sobol
@reexport using NonconvexCore
using NonconvexCore: reset!, isfeasible, @params, VecModel, AbstractOptimizer, AbstractResult
import NonconvexCore: optimize!, Workspace

include("hyperopt.jl")
include("deflation.jl")

end


text:
module TestJuliaManager

using Test

if lowercase(get(ENV, "CI", "false")) == "true"
    @testset begin include("destructive_tests.jl") end
end

end  # module


text:
function train_layer_1_sigmoid_fast(neurons,vars,obs,n_per_part,inds_all,x,y)

    a_all=Vector{Vector{Float64}}()
    ooops=0
    for i=1:neurons
        # if i in 1:10:neurons println(Dates.format(now(), "HH:MM:SS")," Calculating weights w for ",i," neuron, of ",neurons," Total.") end
        try
            atmp=[ones(n_per_part[i+1]-n_per_part[i]) x[inds_all[n_per_part[i]+1:n_per_part[i+1]],:]]\isigm1.(y[inds_all[n_per_part[i]+1:n_per_part[i+1]]])
            obs1=length(inds_all[n_per_part[i]+1:n_per_part[i+1]])
            if sum(isnan.(atmp))==0 && sum(isinf.(atmp))==0 && sum(abs.(atmp[2:end]))>0
                push!(a_all,atmp)
            else
                ooops+=1
            end
        catch ex
            ooops+=1
            println("problem neuron=",i," total ooops=",ooops," ex=",ex)
            push!(a_all,zeros(vars+1))
        end
    end
    layer1=sigm1.([ones(obs) x]*reduce(hcat, a_all))
    mat1=inv([layer1 ones(obs)]'*[layer1 ones(obs)])
    # println(Dates.format(now(), "HH:MM:SS")," Solving system (",obs," by ",size(layer1,2)+1,") for weights v in the output Layer")
    a_layer1=mat1*[layer1 ones(obs)]'*y
    # println(Dates.format(now(), "HH:MM:SS")," All weights computed.")
    return a_all,a_layer1,layer1,mat1

end


text:
# Common

compresscoo(a::AbstractSparseMatrix,args...) = compresscoo(typeof(a),args...)

struct SubSparseMatrix{T,A,B,C} <: AbstractMatrix{T}
  parent::A
  indices::B
  inv_indices::C
  flag::Tuple{Int,Int}
  function SubSparseMatrix(
    parent::AbstractSparseMatrix{T},
    indices::Tuple,
    inv_indices::Tuple,
    flag::Tuple{Int,Int}) where T

    A = typeof(parent)
    B = typeof(indices)
    C = typeof(inv_indices)
    new{T,A,B,C}(parent,indices,inv_indices,flag)
  end
end

Base.size(a::SubSparseMatrix) = map(length,a.indices)
Base.IndexStyle(::Type{<:SubSparseMatrix}) = IndexCartesian()
function Base.getindex(a::SubSparseMatrix,i::Integer,j::Integer)
  I = a.indices[1][i]
  J = a.indices[2][j]
  a.parent[I,J]
end

struct EltypeVector{T,V} <: AbstractVector{T}
  parent::V
  function EltypeVector(::Type{T},parent::V) where {T,V<:AbstractVector}
    new{T,V}(parent)
  end
end
EltypeVector(::Type{T},parent::AbstractVector{T}) where T = parent
Base.size(v::EltypeVector) = size(v.parent)
Base.axes(v::EltypeVector) = axes(v.parent)
Base.@propagate_inbounds Base.getindex(v::EltypeVector{T},i::Integer) where T = convert(T,v.parent[i])
Base.@propagate_inbounds Base.setindex!(v::EltypeVector,w,i::Integer) = (v.parent[i] = w)
Base.IndexStyle(::Type{<:EltypeVector{T,V}}) where {T,V} = IndexStyle(V)

# interface methods

indextype(::Type{<:AbstractSparseMatrix}) = @notimplemented
indextype(a::AbstractSparseMatrix) = indextype(typeof(a))

function LinearAlgebra.mul!(
  C::AbstractVector,
  A::SubSparseMatrix,
  B::AbstractVector,
  α::Number,
  β::Number)

  @notimplemented
end

function compresscoo(
  ::Type{<:AbstractSparseMatrix},
  I::AbstractVector,
  J::AbstractVector,
  V::AbstractVector,
  n::Integer,
  m::Integer,
  combine=+)
  @notimplemented
end

function nzindex(A::AbstractSparseMatrix, i0::Integer, i1::Integer)
  @notimplemented
end

nziterator(a::AbstractSparseMatrix) = @notimplemented

# Implementation of SparseMatrixCSC

indextype(a::Type{SparseMatrixCSC{Tv,Ti}}) where {Tv,Ti} = Ti

function compresscoo(
  ::Type{SparseMatrixCSC{Tv,Ti}},
  I::AbstractVector,
  J::AbstractVector,
  V::AbstractVector,
  m::Integer,
  n::Integer,
  combine=+) where {Tv,Ti}

  sparse(
    EltypeVector(Ti,I),
    EltypeVector(Ti,J),
    EltypeVector(Tv,V),
    m,n,combine)
end

function nzindex(A::SparseArrays.AbstractSparseMatrixCSC, i0::Integer, i1::Integer)
    if !(1 <= i0 <= size(A, 1) && 1 <= i1 <= size(A, 2)); throw(BoundsError()); end
    ptrs = SparseArrays.getcolptr(A)
    r1 = Int(ptrs[i1])
    r2 = Int(ptrs[i1+1]-1)
    (r1 > r2) && return -1
    r1 = searchsortedfirst(rowvals(A), i0, r1, r2, Base.Order.Forward)
    ((r1 > r2) || (rowvals(A)[r1] != i0)) ? -1 : r1
end

nziterator(a::SparseArrays.AbstractSparseMatrixCSC) = NZIteratorCSC(a)

struct NZIteratorCSC{A}
  matrix::A
end

Base.length(a::NZIteratorCSC) = nnz(a.matrix)
Base.eltype(::Type{<:NZIteratorCSC{A}}) where A = Tuple{Int,Int,eltype(A)}
Base.eltype(::T) where T <: NZIteratorCSC = eltype(T)
Base.IteratorSize(::Type{<:NZIteratorCSC}) = Base.HasLength()
Base.IteratorEltype(::Type{<:NZIteratorCSC}) = Base.HasEltype()

@inline function Base.iterate(a::NZIteratorCSC)
  if nnz(a.matrix) == 0
    return nothing
  end
  col = 0
  knext = nothing
  while knext === nothing
    col += 1
    ks = nzrange(a.matrix,col)
    knext = iterate(ks)
  end
  k, kstate = knext
  i = Int(rowvals(a.matrix)[k])
  j = col
  v = nonzeros(a.matrix)[k]
  (i,j,v), (col,kstate)
end

@inline function Base.iterate(a::NZIteratorCSC,state)
  col, kstate = state
  ks = nzrange(a.matrix,col)
  knext = iterate(ks,kstate)
  if knext === nothing
    while knext === nothing
      if col == size(a.matrix,2)
        return nothing
      end
      col += 1
      ks = nzrange(a.matrix,col)
      knext = iterate(ks)
    end
  end
  k, kstate = knext
  i = Int(rowvals(a.matrix)[k])
  j = col
  v = nonzeros(a.matrix)[k]
  (i,j,v), (col,kstate)
end

function LinearAlgebra.mul!(
  C::AbstractVector,
  A::SubSparseMatrix{T,<:SparseArrays.AbstractSparseMatrixCSC} where T,
  B::AbstractVector,
  α::Number,
  β::Number)

  size(A, 2) == size(B, 1) || throw(DimensionMismatch())
  size(A, 1) == size(C, 1) || throw(DimensionMismatch())
  size(B, 2) == size(C, 2) || throw(DimensionMismatch())
  if β != 1
      β != 0 ? rmul!(C, β) : fill!(C, zero(eltype(C)))
  end
  rows, cols = A.indices
  invrows, invcols = A.inv_indices
  rflag, cflag = A.flag
  Ap = A.parent
  nzv = nonzeros(Ap)
  rv = rowvals(Ap)
  for (j,J) in enumerate(cols)
      αxj = B[j] * α
      for p in nzrange(Ap,J)
        I = rv[p]
        i = invrows[I]*rflag
        if i>0
          C[i] += nzv[p]*αxj
        end
      end
  end
  C
end

# Implementation of SparseMatrixCSR

indextype(a::Type{SparseMatrixCSR{Bi,Tv,Ti}}) where {Bi,Tv,Ti} = Ti

function compresscoo(
  ::Type{SparseMatrixCSR{Bi,Tv,Ti}},
  I::AbstractVector,
  J::AbstractVector,
  V::AbstractVector,
  m::Integer,
  n::Integer,
  combine=+) where {Bi,Tv,Ti}

  sparsecsr(
    Val(Bi),
    EltypeVector(Ti,I),
    EltypeVector(Ti,J),
    EltypeVector(Tv,V),
    m,n,combine)
end

function nzindex(A::SparseMatrixCSR, i0::Integer, i1::Integer)
  if !(1 <= i0 <= size(A, 1) && 1 <= i1 <= size(A, 2)); throw(BoundsError()); end
  o = getoffset(A)
  Bi = getBi(A)
  r1 = Int(A.rowptr[i0]+o)
  r2 = Int(A.rowptr[i0+1]-Bi)
  (r1 > r2) && return -1
  i1o = i1-o
  k = searchsortedfirst(colvals(A), i1o, r1, r2, Base.Order.Forward)
  ((k > r2) || (colvals(A)[k] != i1o)) ? -1 : k
end

function LinearAlgebra.mul!(
  C::AbstractVector,
  A::SubSparseMatrix{T,<:SparseMatrixCSR} where T,
  B::AbstractVector,
  α::Number,
  β::Number)

  size(A, 2) == size(B, 1) || throw(DimensionMismatch())
  size(A, 1) == size(C, 1) || throw(DimensionMismatch())
  size(B, 2) == size(C, 2) || throw(DimensionMismatch())
  if β != 1
    β != 0 ? rmul!(C, β) : fill!(C, zero(eltype(C)))
  end
  rows, cols = A.indices
  invrows, invcols = A.inv_indices
  rflag, cflag = A.flag
  Ap = A.parent
  nzv = nonzeros(Ap)
  cv = colvals(Ap)
  o = getoffset(Ap)
  for (i,I) in enumerate(rows)
    @inbounds for p in nzrange(Ap,I)
      J = cv[p]+o
      j = invcols[J]*cflag
      if j>0
        C[i] += nzv[p]*B[j]*α
      end
    end
  end
  C
end

nziterator(a::SparseMatrixCSR) = NZIteratorCSR(a)

struct NZIteratorCSR{A}
  matrix::A
end

Base.length(a::NZIteratorCSR) = nnz(a.matrix)
Base.eltype(::Type{<:NZIteratorCSR{A}}) where A = Tuple{Int,Int,eltype(A)}
Base.eltype(::T) where T <: NZIteratorCSR = eltype(T)
Base.IteratorSize(::Type{<:NZIteratorCSR}) = Base.HasLength()
Base.IteratorEltype(::Type{<:NZIteratorCSR}) = Base.HasEltype()

@inline function Base.iterate(a::NZIteratorCSR)
  if nnz(a.matrix) == 0
    return nothing
  end
  row = 0
  ptrs = a.matrix.rowptr
  knext = nothing
  while knext === nothing
    row += 1
    ks = nzrange(a.matrix,row)
    knext = iterate(ks)
  end
  k, kstate = knext
  i = row
  j = Int(colvals(a.matrix)[k]+getoffset(a.matrix))
  v = nonzeros(a.matrix)[k]
  (i,j,v), (row,kstate)
end

@inline function Base.iterate(a::NZIteratorCSR,state)
  row, kstate = state
  ks = nzrange(a.matrix,row)
  knext = iterate(ks,kstate)
  if knext === nothing
    while knext === nothing
      if row == size(a.matrix,1)
        return nothing
      end
      row += 1
      ks = nzrange(a.matrix,row)
      knext = iterate(ks)
    end
  end
  k, kstate = knext
  i = row
  j = Int(colvals(a.matrix)[k]+getoffset(a.matrix))
  v = nonzeros(a.matrix)[k]
  (i,j,v), (row,kstate)
end


text:
#=
This module is based on the design procedure described in Designing Planar Magnetics.
The introduction approximatly follows the examples in this paper.

Load the module.
=#
using PlanerTransformer

#=
Define the operating frequency and maximum specific power loss for the transformer
=#
frequency = 0.5e6 # Hz
spl_max = 150e3 # W/m^3

#=
Standard core geometries are available in `core_geometry_dict`.
  An ER25 core / plate set will be used.
=#
core = core_geometry_dict["er25/5.5_plt"]

#=
The properities of some ferrites are available in `ferrite_dict`.  The functions
`specific_power_loss` and `flux_density` interpolate the specific power loss curves
from the data sheet.

Use `flux_density` to calculate the maximum flux density for given specific power loss
for two ferrites.
=#
flux_density_max_pp_3f35 = flux_density(ferrite_dict["3f35"], spl_max, frequency) * 2.0
flux_density_max_pp_r = flux_density(ferrite_dict["r"], spl_max, frequency) * 2.0

#=

Calculate volt seconds per turn.
=#
vspt_3f35 = volt_seconds_per_turn(core, flux_density_max_pp_3f35)
vspt_r = volt_seconds_per_turn(core, flux_density_max_pp_r)

#=

Maximum output voltage for forward converter with 1T secondary @ 50% duty cycle
=#
max_vo_3f35 = vspt_3f35 * frequency
max_vo_r = vspt_r * frequency
#=
Define the windings.

`PlanerTrnasformer` only allows transformers with two windings.  First the PCB
to be used for the windings is defined with `PCB_Specification`.  `windings` then
creats a `Windings` object.  All primary layers must have the same number of turns.
All secondary layers must have the same number of turns.  Primary and secondary layers
may be in series or parallel.  This transformer will be built with a double sided PCB
with 4oz copper.
=#
trace_edge_gap = 0.5e-3
trace_trace_gap = 0.4e-3
copper_thickness = copper_weight_to_meters(4.0) # 4.0oz copper
dielectric_thickness = 0.125e-3
stackup = Stackup((copper, fr408, copper),(copper_thickness, dielectric_thickness, copper_thickness))
pcb = PCB_Specification(trace_edge_gap,
                        trace_trace_gap,
                        stackup)
primary_turns_per_layer = 6
secondary_turns_per_layer = 1
is_layer_primary = (true, false) # which layers belong to primary and secondary
are_primary_layers_in_series = false
are_secondary_layers_in_series = false
my_windings = windings(pcb, core,
                    primary_turns_per_layer,
                    secondary_turns_per_layer,
                    is_layer_primary,
                    are_primary_layers_in_series,
                    are_secondary_layers_in_series)
#=
Maximum voltage at the output of the converter (not the peak voltage at the transformer output)
=#
max_volts = frequency.*volt_seconds(my_windings,flux_density_max_pp_3f35)
#=
Create a transformer by combining a ferrite with windings.
=#
transformer_3f35_2layer = Transformer(ferrite_dict["3f35"], my_windings)
#=
Verify 6:1 turns ratio.
=#
turns(transformer_3f35_2layer) # (primary, secondary)
#=
Laekage inductance is referenced to primary.
=#
leakage_inductance(transformer_3f35_2layer)
#=
Define parameters for the switching power supply using the transformer.  Power must flow
from primary to secondary (v_in positive, i_out negative).  Voltages and currents are
peak+ at the transformer, not the input and output of the switching power supply.
=#
v_in = 48.0
i_out = -7.0
duty = 0.50
my_converter = Forward(v_in, i_out, frequency, duty) # or PushPull
#=
Calculate power dissipation, and other stuff,  of a transformer in a
 switching power supply.    Voltages and currents are peak at the transformer,
not the input and output of the switching power supply.  Powers are averge.
=#
tpa = transformer_power_analysis(transformer_3f35_2layer, my_converter);
#=
Many parameters can be displayed
=#
voltage(tpa) # (primary, secondary)
current(tpa) # (primary, secondary)
input_power(tpa)
output_power(tpa)
ac_winding_resistance(tpa) # (primary, secondary)
dc_winding_resistance(tpa) # (primary, secondary)
equilivent_resistance(tpa) # (primary, secondary)
core_specific_power(tpa) # W/m^3
core_total_power(tpa)
winding_power(tpa) # (primary, secondary)
total_power(tpa)
#=
Total power dissipation with 20A output
=#
total_power(transformer_power_analysis(transformer_3f35_2layer, Forward(v_in,-20.0,frequency,duty)))
#=
The next three sections Double and interleave the Windings,
Copper strip secondary, and Primary in series are a series of improvments on
this transformer.
=#
#=
Double and interleave the Windings
=#

b_stage_thickness = 0.0762
double_stackup = Stackup((copper,           fr408,                copper,           fr408,             copper,           fr408,                copper),
                         (copper_thickness, dielectric_thickness, copper_thickness, b_stage_thickness, copper_thickness, dielectric_thickness, copper_thickness))
double_pcb = PCB_Specification(trace_edge_gap,
                               trace_trace_gap,
                               double_stackup)
double_windings = windings(double_pcb, core,
                           primary_turns_per_layer,
                           secondary_turns_per_layer,
                           prisec"S-P-P-S", # use prisec string to specify which layers are primary and secondary.
                           are_primary_layers_in_series,
                           are_secondary_layers_in_series)

double_transformer_3f35 = Transformer(ferrite_dict["3f35"], double_windings)
turns(double_transformer_3f35)
leakage_inductance(double_transformer_3f35)
double_tpa = transformer_power_analysis(double_transformer_3f35, my_converter)
equilivent_resistance(double_tpa) # (primary, secondary)
core_total_power(double_tpa)
winding_power(double_tpa) # (primary, secondary)
total_power(double_tpa)
total_power(transformer_power_analysis(double_transformer_3f35, Forward(v_in,-20.0,frequency,duty)))
#=
Copper strip secondary
=#
strip_stackup = Stackup((copper,           fr408,                copper,           fr408,             copper,           fr408,                copper),
                         (copper_thickness+1e-3, dielectric_thickness, copper_thickness, b_stage_thickness, copper_thickness, dielectric_thickness, copper_thickness+1e-3))
strip_pcb = PCB_Specification(trace_edge_gap,
                              trace_trace_gap,
                              strip_stackup)
strip_windings = windings(strip_pcb, core,
                          primary_turns_per_layer,
                          secondary_turns_per_layer,
                          prisec"S-P-P-S",
                          are_primary_layers_in_series,
                          are_secondary_layers_in_series)

strip_transformer_3f35 = Transformer(ferrite_dict["3f35"], strip_windings)
turns(strip_transformer_3f35)
leakage_inductance(strip_transformer_3f35)
strip_tpa = transformer_power_analysis(strip_transformer_3f35, my_converter)
equilivent_resistance(strip_tpa) # (primary, secondary)
core_total_power(strip_tpa)
winding_power(strip_tpa) # (primary, secondary)
total_power(strip_tpa)
total_power(transformer_power_analysis(strip_transformer_3f35, Forward(v_in,-20.0,frequency,duty)))
#=
Primary in series
=#
series_windings = windings(strip_pcb, core,
                          3, # turns per layer, primary
                          secondary_turns_per_layer,
                          prisec"S-P-P-S",
                          true, # primary layers are in series
                          are_secondary_layers_in_series)

series_transformer_3f35 = Transformer(ferrite_dict["3f35"], series_windings)
series_tpa = transformer_power_analysis(series_transformer_3f35, my_converter)
equilivent_resistance(series_tpa) # (primary, secondary)
core_total_power(series_tpa)
winding_power(series_tpa) # (primary, secondary)
total_power(series_tpa)
total_power(transformer_power_analysis(series_transformer_3f35, Forward(v_in,-20.0,frequency,duty)))
#=
LTspice model

Parameters for a cicuit simulation can be extracted.
=#
turns(series_transformer_3f35)
leakage_inductance(series_transformer_3f35)
chan_inductor(series_transformer_3f35)
equilivent_resistance(series_tpa)


text:
using DelimitedFiles

file = "input04.txt"
blocks = DelimitedFiles.readdlm(file, skipstart=2)
num_blocks = Int(size(blocks, 1) / size(blocks, 2))
block_len = size(blocks, 2)

f = open(file, "r")
    l = readlines(f)
close(f)
rng = split(l[1], ",")

function iswinner(b)
    for row in eachrow(b)
        if count(x->x==1, row) == length(row)
            return true
        end
    end
    for row in eachrow(transpose(b))
        if count(x->x==1, row) == length(row)
            return true
        end
    end
    return false
end

markers = zeros(Int, size(blocks, 1), size(blocks, 2))

function play()
for draw in rng
    for (i, x) in enumerate(blocks)
        row = mod(i, size(blocks, 1))
        if row == 0
            row = size(blocks, 1)
        end
        column = trunc(Int, ((i - 1) / size(blocks, 1))) +1
        #println("row $row")
        #println("column $column")
        #println(x)
        if x == parse(Int, draw)
            markers[row, column] = 1
        end
    end
    for j in range(1, length=num_blocks, step=block_len)
        bl = markers[j:j+4,:]
        if iswinner(bl)
            global winner = j
            return draw
            #@goto escape_label
        end
    end
end
end

#@label escape_label
last_draw = play()
winning_mask = markers[winner:winner+4,:]
winning_block = blocks[winner:winner+4,:]
println(winning_mask)
println(winning_block)
println("Last number drawn: $last_draw")

winning_sum = 0
for (i, j) in zip(winning_mask, winning_block)
    if i == 0
        global winning_sum += j
    end
end
println("Winning sum: $winning_sum")
result1 = winning_sum * parse(Int, last_draw)
println("Result part 1: $result1")

markers = zeros(Int, size(blocks, 1), size(blocks, 2))

function play_to_loose()
bingo = []
for draw in rng
    for (i, x) in enumerate(blocks)
        row = mod(i, size(blocks, 1))
        if row == 0
            row = size(blocks, 1)
        end
        column = trunc(Int, ((i - 1) / size(blocks, 1))) +1
        if x == parse(Int, draw)
            markers[row, column] = 1
        end
    end
    for j in range(1, length=num_blocks, step=block_len)
        bl = markers[j:j+4,:]
        if j in bingo
            continue
        elseif iswinner(bl)
            append!(bingo, j)
            if length(bingo) == num_blocks
                global looser = j
                return draw
            end
        end
    end
end
end

final_draw = play_to_loose()
loosing_mask = markers[looser:looser+4,:]
loosing_block = blocks[looser:looser+4,:]
println("Last number drawn: $final_draw")

loosing_sum = 0
for (i, j) in zip(loosing_mask, loosing_block)
    if i == 0
        global loosing_sum += j
    end
end
println("Loosing: $loosing_sum")
result2 = loosing_sum * parse(Int, final_draw)
println("Result part 2: $result2")

max_stars_repo_name: jung-benjamin/advent-of-code

max_issues_repo_name: jung-benjamin/advent-of-code

max_forks_repo_name: jung-benjamin/advent-of-code


text:
module NeurDS

using ModelingToolkit, OrdinaryDiffEq, Plots

# helper functions for ODE system
r(v)= 127v/105 +8265/105
m∞(v) = 1/(exp((-22-v)/8)+1)
H(x) = x>0. ? 1. : 0.

function PlantCell(;name)
    @parameters μh μn μx μCa gna gn gx gCa xs Cas gleak
    sts = @variables t v(t) n(t) h(t) x(t) Ca(t) Iext(t)
    D = Differential(t)

    ODESystem(Equation[
        D(v) ~ gna*^(m∞(v),3.0)*h*(30-v)+gn*^(n,4.0)*(-75 -v)+gx*x*(30-v)+gCa*Ca/(.5+Ca)*(-75-v)+gleak*(-40-v) + Iext,
        D(h) ~ μh*((1-h)*(.07*exp((25-r(v))/20))-h*(1/(1+exp((55-r(v))/10)))),
        D(n) ~ μn*((1-n)*(.01*(55-r(v)))/(exp((55-r(v))/10)-1)-n*0.125*exp((45 - r(v))/80)),
        D(x) ~ μx*((1/(exp(.15*(-50 + xs -v))+1))-x),
        D(Ca) ~ μCa*(.0085*x*(140 - v +Cas)-Ca),
    ]; name, defaults = Dict([
        #Initial Conditions
        v => 0.,
        n => 0.,
        h => 0.,
        x => 0.9,
        Ca => .7,
        #parameters
        μh => .08,
        μn => .08,
        gna => 4.,
        gn => 3.,
        xs => -4.,
        Cas => -60.,
        gleak => .0025,
        μx => .012,
        μCa => .00025
    ]))
end

@named c1 = PlantCell()
@named c2 = PlantCell()

struct AlphaSynapse
    sys
    pre
    post
end

function AlphaSynapse(pre, post; name)
    @parameters α β threshold k E g
    sts = @variables t s(t) v(t) vpre(t)
    D = Differential(t)
    sys = ODESystem(Equation[
        D(s) ~ α*(1-s)/(1+exp(threshold-vpre)/k) - β*s,
        v ~ g*s*(E-v)
    ]; name, defaults = Dict([
        #initial conditions
        s => 0.,
        #parameters
        α => .1,
        β => .1,
        threshold => 0.,
        k => 0.,
        E => 30.,
        g => .1
    ]))
    AlphaSynapse(sys, pre, post)
end

@named s = AlphaSynapse(c1, c2)

function network(cells, synapses; name)
    external_currents = []
    for cell in cells
        iexts = [syn.sys.v for syn in synapses if syn.post.name == cell.name]
        push!(external_currents, cell.Iext ~ isempty(iexts) ? 0 : sum(iexts))
    end
    presynaptic_drives = [syn.sys.vpre ~ syn.pre.v for syn in synapses]
    connections = vcat(external_currents, presynaptic_drives)

    compose(ODESystem(connections), cells..., [syn.sys for syn in synapses]...)
end

@named net = network([c1,c2], [s])

structural_simplify(net)


text:
# Sort for dicts
import Base: sort, sort!

function sort!(d::OrderedDict; byvalue::Bool=false, args...)
    if d.ndel > 0
        rehash!(d)
    end

    if byvalue
        p = sortperm(d.vals; args...)
    else
        p = sortperm(d.keys; args...)
    end

    for (i,key) in enumerate(d.keys)
        idx = ht_keyindex(d, key, false)
        d.slots[idx] = p[i]
    end

    d.keys = d.keys[p]
    d.vals = d.vals[p]

    return d
end

sort(d::OrderedDict; args...) = sort!(copy(d); args...)
sort(d::Dict; args...) = sort!(OrderedDict(d); args...)
## Uncomment these after #224 is merged
# sort!(d::DefaultOrderedDict; args...) = (sort!(d.d.d; args...); d)
# sort(d::DefaultDict; args...) = DefaultOrderedDict(d.d.default, sort(d.d.d; args...))
# sort(d::DefaultOrderedDict; args...) = DefaultOrderedDict(d.d.default, sort!(copy(d.d.d); args...))


text:
immutable LocationIterator{w,T}
    range::UnitRange{Int}
    index::FMIndex{w,T}
end

Base.length(iter::LocationIterator) = length(iter.range)
@inline Base.start(iter::LocationIterator) = 1
@inline Base.done(iter::LocationIterator, i) = i > length(iter)
@inline Base.next(iter::LocationIterator, i) = sa_value(iter.range[i], iter.index) + 1, i + 1


text:
#=
    Given three points (in their natural sequence) and a contained or nearby fourth abcissa,
    interpolate or extrapolate the fourth ordinate(s) using a polynomial or a rational approach.

        quadratic(pt₁, pt₂, pt₃, x)

        thiele3(pt₁, pt₂, pt₃, x)
=#

function quadratic(pt₁, pt₂, pt₃, x)
    result = quadratic(pt₁[1], pt₂[1], pt₃[1],
                       pt₁[2:end], pt₂[2:end], pt₃[2:end], x)
    # check for degenerate case
    if !allfinite(result)
	result = linear(pt₁, pt₂, pt₃, x)
    end
	return result
end

function quadratic(pt₁x, pt₂x, pt₃x, pt₁y, pt₂y, pt₃y, x)
    t1 = @. pt₂x - pt₁x
    t2 = @. pt₁x - pt₃x
    t3 = @. pt₃x - pt₂x
    t4 = @. t2   * pt₂y
    t5 = @. pt₃x * pt₃x
    t6 = @. pt₂x * pt₂x
    t7 = @. pt₁x * pt₁x

    s = @. -(inv(t1) * inv(t2) * inv(t3))

    a = @. (pt₁y * t3 + pt₃y * t1 + t4) * x
    b = @. t5 * (pt₁y - pt₂y)
    c = @. t6 * (pt₃y - pt₁y)
    d = @. t7 * (pt₂y - pt₃y)
    q = @. t6 * (pt₁y * pt₃x - pt₃y * pt₁x)
    r = @. t4 *  pt₁x * pt₃x
    n = @. (-pt₁y * t5 + pt₃y * t7) * pt₂x
    aa = @. a - b - c - d
    bb = @. r - n - q

    result = @. (s * (aa * x + bb))
    return result # may contain !finite values
end

function thiele3(pt₁, pt₂, pt₃, x)
    result = thiele3(pt₁[1], pt₂[1], pt₃[1],
                     pt₁[2:end], pt₂[2:end], pt₃[2:end], x)
    # check for degenerate case
    if !allfinite(result)
        result = quadratic(pt₁, pt₂, pt₃, x)
    end
    return result
end

function thiele3(pt₁x, pt₂x, pt₃x, pt₁y, pt₂y, pt₃y, x)
    t1 = @. pt₁y - pt₂y
    t2 = @. pt₂y - pt₃y
    t1 = @. inv(t1)
    t2 = @. inv(t2)
    t1 = @. (pt₁x - pt₂x) * t1
    t2 = @. -(pt₂x - pt₃x) * t2 + t1
    t2 = @. inv(t2)
    t2 = @. (pt₁x - pt₃x) * t2 - pt₁y + pt₂y
    t2 = @. inv(t2)
    t1 = @. -(pt₂x - x) * t2 + t1
    t1 = @. inv(t1)

    result = @. (-(pt₁x - x) * t1 + pt₁y)
    if !allfinite(result)
        result = (pt₁x, pt₂x, pt₃x, pt₁y, pt₂y, pt₃y, x)
    end

    return result
end


text:
@testset "twix input svs_lcosy" begin
    twix = load_twix("twix/sub-SiemensBrainPhantom_seq-svslcosy_inc-1.twix")
    @test length(twix.data) == 1
    acq = twix.data[1]
    @test size(acq.data) == (2080,34)
    @test eltype(acq.data) == ComplexF32
    @test acq.cutoff_pre  == 0x0002
    @test acq.cutoff_post == 0x001e

    @test sprint(show, twix) == """
        MRExperiment metadata:
          Protocol           = svs_lcosy
          Sequence File Name = %CustomerSeq%\\svs_lcosy
          Software Version   = N4_VE11C_LATEST_20160120
          Device Id          = ("Siemens", "Prisma", "166042")
          Reference Date     = 2018-11-27T12:22:52
          Frequency          = 123255189 Hz
          Coils              = ["HeadNeck_64"]
        Acquisition summary:
          Number   = 1
          Duration = 0.0 s
        """
end

@testset "twix input press" begin
    # Avg = 1
    twix = load_twix("twix/sub-SiemensBrainPhantom_seq-svsse_ref-1_avg-1.twix")
    @test length(twix.data) == 2
    acq = twix.data[2]
    @test size(acq.data) == (2080,34)
    @test eltype(acq.data) == ComplexF32
    @test acq.cutoff_pre  == 0x0000
    @test acq.cutoff_post == 0x0020

    @test sprint(show, twix) == """
        MRExperiment metadata:
          Protocol           = svs_se-ref-1_avg-1
          Sequence File Name = %SiemensSeq%\\svs_se
          Software Version   = N4_VE11C_LATEST_20160120
          Device Id          = ("Siemens", "Prisma", "166042")
          Reference Date     = 2019-08-22T15:00:28
          Frequency          = 123254849 Hz
          Coils              = ["HeadNeck_64"]
        Acquisition summary:
          Number   = 2
          Duration = 2.0 s
          Loop index 6 in [0,1]   (phase)
        """
    expt = mr_load(twix)
    @test length(expt.press_scans) == 1
    @test length(expt.navigator) == 0
    @test length(expt.ref_scans) == 1

    # Avg = 2
    twix = load_twix("twix/sub-SiemensBrainPhantom_seq-svsse_ref-1_avg-2.twix")
    @test length(twix.data) == 3
    all_acq = twix.data[2:end]
    for acq in all_acq
        @test size(acq.data) == (2080,34)
        @test eltype(acq.data) == ComplexF32
        @test acq.cutoff_pre  == 0x000c
        @test acq.cutoff_post == 0x0014
    end

    @test sprint(show, twix) == """
        MRExperiment metadata:
          Protocol           = svs_se-ref-1_avg-2
          Sequence File Name = %SiemensSeq%\\svs_se
          Software Version   = N4_VE11C_LATEST_20160120
          Device Id          = ("Siemens", "Prisma", "166042")
          Reference Date     = 2019-08-22T15:00:28
          Frequency          = 123254849 Hz
          Coils              = ["HeadNeck_64"]
        Acquisition summary:
          Number   = 3
          Duration = 4.0 s
          Loop index 2 in [0,1]   (acquisition)
          Loop index 6 in [0,1]   (phase)
        """
    expt = mr_load(twix)
    @test length(expt.press_scans) == 2
    @test length(expt.navigator) == 0
    @test length(expt.ref_scans) == 1
end

@testset "twix quality control" begin
    twix = @test_logs (Logging.Warn,r"Unexpected empty meas packet") #=
               =# (Logging.Warn,r"Different values found in MeasYaps") #=
               =# load_twix("twix/sub-SiemensBrainPhantom_seq-svslcosy_incomplete.twix")
    @test MagneticResonanceSignals.AcquisitionsIncomplete in twix.quality_control

    valid_twix_bytes = read("twix/sub-SiemensBrainPhantom_seq-svslcosy_inc-1.twix")

    partially_zeroed_twix = copy(valid_twix_bytes)
    partially_zeroed_twix[0x69:end] .= 0
    twix = @test_logs (Logging.Warn,r"Unexpected empty measurement header sections") match_mode=:any #=
        =# load_twix(IOBuffer(partially_zeroed_twix))
    @test MagneticResonanceSignals.MeasHeaderEmpty in twix.quality_control

    truncated_twix = copy(valid_twix_bytes)
    truncated_twix = truncated_twix[1:900_000]
    twix = @test_logs (Logging.Warn,r"Twix acquisition truncated at position 900000") #=
        =# (Logging.Warn,r"Different values found in MeasYaps") #=
        =# load_twix(IOBuffer(truncated_twix))
    @test MagneticResonanceSignals.AcquisitionsIncomplete in twix.quality_control
end

@testset "metadata parsing" begin
    twix = load_twix("twix/sub-SiemensBrainPhantom_seq-svslcosy_inc-1.twix")
    @test twix.metadata["tProtocolName"] == "svs_lcosy"
    @test twix.metadata["sTXSPEC.asNucleusInfo[0].lFrequency"] == 123255189
    @test twix.metadata["sTXSPEC.asNucleusInfo[0].lFrequency"] == 123255189

    @test twix.metadata["Dicom.ManufacturersModelName"] == "Prisma"
    @test twix.metadata["Dicom.DeviceSerialNumber"] == "166042"

    @test twix.metadata["Meas.tReferenceImage0"] == "1.3.12.2.1107.5.2.43.166042.2018112712225215359633498"
    @test twix.metadata["Meas.tReferenceImage1"] == "1.3.12.2.1107.5.2.43.166042.2018112712230494362433518"
    @test twix.metadata["Meas.tReferenceImage2"] == "1.3.12.2.1107.5.2.43.166042.2018112712231775153533538"
    @test twix.metadata["Meas.tFrameOfReference"] == "1.3.12.2.1107.5.2.43.166042.1.20181127122016414.0.0.4988"
end

@testset "data extraction" begin
    twix = load_twix("twix/sub-SiemensBrainPhantom_seq-svslcosy_inc-1.twix")
    samples = sampledata(twix, 1) # downsample==1
    @test size(samples) == (2048,34)
    @test first(getaxis(samples, :time)) == 0u"μs"
    @test step(getaxis(samples, :time)) == 500u"μs"
    # This test data used a 64-channel head and neck coil, with a subset of
    # coil elements activated
    @test getaxis(samples, :channel) == [
        :H42,:H43,:H31,:H30,:H63,:H62,:H48,:H49,:H35,:H34,:H40,:H41,:H50,:H70,:H45,
        :H44,:H67,:H66,:H52,:H51,:H46,:H47,:H54,:H53,:H65,:H64,:H71,:H55,:H33,:H32,
        :H61,:H60,:H69,:H68
    ]
    # Test values of first and last time samples of channel data
    @test samples[1,1]     == -1.3986137f-6 + 1.442153f-6im
    @test samples[1,end]   ==  2.022367f-6  - 5.987473f-6im
    @test samples[end,1]   == -1.3408717f-6 - 3.1664968f-7im
    @test samples[end,end] == -3.306195f-8  - 7.7229924f-7im

    # With downsampling
    samples = sampledata(twix, 1, downsample=2)
    @test size(samples) == (1024,34)
    @test first(getaxis(samples, :time)) == 0u"μs"
    @test step(getaxis(samples, :time)) == 1000u"μs"
end

@testset "coil select parsing" begin
    yaps_dict = MagneticResonanceSignals.parse_header_yaps(String(read("twix/PRESS_TE135_Breast_Coil_Headers/MeasYaps")))
    @test length(@test_logs((:warn, "Could not find some metadata for coil 7"),
                            MagneticResonanceSignals.parse_yaps_rx_coil_selection(yaps_dict))) == 7
end

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


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

mutable struct ConnectionMonitorOutput <: SwaggerModel
    type::Any # spec type: Union{ Nothing, String } # spec name: type
    workspaceSettings::Any # spec type: Union{ Nothing, ConnectionMonitorWorkspaceSettings } # spec name: workspaceSettings

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

const _property_map_ConnectionMonitorOutput = Dict{Symbol,Symbol}(Symbol("type")=>Symbol("type"), Symbol("workspaceSettings")=>Symbol("workspaceSettings"))
const _property_types_ConnectionMonitorOutput = Dict{Symbol,String}(Symbol("type")=>"String", Symbol("workspaceSettings")=>"ConnectionMonitorWorkspaceSettings")
Base.propertynames(::Type{ ConnectionMonitorOutput }) = collect(keys(_property_map_ConnectionMonitorOutput))
Swagger.property_type(::Type{ ConnectionMonitorOutput }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ConnectionMonitorOutput[name]))}
Swagger.field_name(::Type{ ConnectionMonitorOutput }, property_name::Symbol) =  _property_map_ConnectionMonitorOutput[property_name]

const _allowed_ConnectionMonitorOutput_type = ["Workspace"]

function check_required(o::ConnectionMonitorOutput)
    true
end

function validate_property(::Type{ ConnectionMonitorOutput }, name::Symbol, val)
    if name === Symbol("type")
        Swagger.validate_param(name, "ConnectionMonitorOutput", :enum, val, _allowed_ConnectionMonitorOutput_type)
    end
end


text:
module GeodesicRendering

import Base.Threads: @threads

using Parameters

import SciMLBase

using GeodesicBase
using GeodesicTracer

include("utility.jl")
include("render-cache.jl")
include("point-functions.jl")
include("render.jl")

function rendergeodesics(
    m::AbstractMetricParams{T},
    init_pos,
    max_time;
    vf = ConstValueFunctions.shadow,
    kwargs...,
) where {T}
    __rendergeodesics(
        m,
        init_pos;
        image_width = 350,
        image_height = 250,
        fov_factor = 3.0,
        max_time = convert(T, max_time),
        vf = vf,
        kwargs...,
    )
end

function prerendergeodesics(
    m::AbstractMetricParams{T},
    init_pos,
    max_time;
    cache = DefaultCache(),
    kwargs...,
) where {T}
    __prerendergeodesics(
        m,
        init_pos,
        cache;
        image_width = 350,
        image_height = 250,
        fov_factor = 3.0,
        max_time = convert(T, max_time),
        kwargs...,
    )
end

export rendergeodesics, prerendergeodesics

end # module


text:
function saveCSV(com::Community, file::String)

    #Save local and identity
    d = DataFrame()
    for (i,j) in enumerate(com.declaredSymbols_["Local"])
        d[!,j] = com.local_[:,i]
    end
    for (i,j) in enumerate(com.declaredSymbols_["Identity"])
        d[!,j] = com.identity_[:,i]
    end

    CSV.write(string(file,"_local.csv"),d)

    #Save global and globalArray
    d = DataFrame()
    d[!,:t] = [com.t]
    d[!,:N] = [com.N]
    for (i,j) in enumerate(com.declaredSymbols_["Global"])
        d[!,j] = com.global_[:,i]
    end
    for (i,j) in enumerate(com.declaredSymbols_["GlobalArray"])
        dims = size(com.globalArray_[i])
        values = vcat(com.globalArray_[i]...)
        positions = vcat(CartesianIndices(com.globalArray_[i])...)
        name = string(j,"_")
        for (pos,p) in enumerate(positions)
            name2 = name
            for c in 1:length(dims)
                name2 = string(name2,"_",p[c])
            end
            d[!,Meta.parse(name2)] = [values[pos]]
        end
    end
    CSV.write(string(file,"_global.csv"),d)

    return
end

function saveCSV(com::CommunityInTime, file::String)

    for (i,j) in enumerate(com.com)
        save(j,string(file,"_",i))
    end

    return
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
push!(LOAD_PATH,"../src/")

using Documenter, MOSInterface

makedocs(
    sitename="MOS Interface",
    modules = [MOSInterface],
    pages=[
        "Getting Started" => "index.md",
        "API Reference" => [
            "interface.md",
            "model.md"
        ],
    ],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    )
)

deploydocs(
    repo = "github.com/Fuinn/mos-interface-jl.git",
)

max_stars_repo_name: Fuinn/mos-interface-jl

max_issues_repo_name: Fuinn/mos-interface-jl

max_forks_repo_name: Fuinn/mos-interface-jl


text:
# This file was generated, do not modify it. # hide
xgb = fitted_params(mtm).best_model
@show xgb.max_depth
@show xgb.min_child_weight


text:
using Plots
plotly()

x = 0:0.01:2pi
y_sin = sin.(x)
y_cos = cos.(x)

plot(x, [y_sin y_cos], labels=["sin" "cos"])

_round(number) = round(number, sigdigits=4)
_point_label(x, y, fun) = "$fun($(_round(x))) = $(_round(y))"
tooltips_sin = [_point_label(x, y, "sin") for (x, y) in zip(x, y_sin)]
tooltips_cos = [_point_label(x, y, "cos") for (x, y) in zip(x, y_cos)]

plot(x, [y_sin y_cos], labels=["sin" "cos"], hover=[tooltips_sin tooltips_cos])


text:
@testset "Test How to State" begin
###############################################################################
#
# The data used through the algorithmic process in the Stopping framework
# are stored in a State.
# We illustrate here the GenericState and its features
#
###############################################################################
#using Test, Stopping

###############################################################################
#The GenericState contains only two entries:
# a vector x, and a Float current_time
state1 = GenericState(ones(2)) #takes a Vector as a mandatory input
state2 = GenericState(ones(2), current_time = 1.0)

#By default if a non-mandatory entry is not specified it is void:
@test isnan(state1.current_time)
@test state2.current_time == 1.0

###############################################################################
#The GenericState has two functions: update! and reinit!
#update! is used to update entries of the State:
update!(state1, current_time = 1.0)
@test state1.current_time == 1.0
#Note that the update select the relevant entries
update!(state1, fx = 1.0) #does nothing as there are no fx entry
@test state1.current_time == 1.0 && state1.x == ones(2)

#The update! can be done only if the new entry is void or has the same type
#as the existing one.
update!(state1, current_time = 2) #does nothing as it is the wrong type
@test state1.current_time == 1.0
#An advanced user can force the update even if the type is not the same by
#turning the keyword convert as true (it is false by default).
#update!(state1, convert = true, current_time = 2) NON!!!
#@test state1.current_time == 2
#Non-required entry in the State can always be set as void without convert
update!(state1, current_time = NaN)
@test isnan(state1.current_time)

#A shorter way to empty the State is to use the reinit! function.
#This function is particularly useful, when there are many entries.
reinit!(state2)
@test state2.x == ones(2) && isnan(state2.current_time)
#If we want to reinit! with a different value of the mandatory entry:
reinit!(state2, zeros(2))
@test state2.x == zeros(2) && isnan(state2.current_time)
#After reinitializing the State reinit! can update entries passed as keywords.
#either in the default call:
reinit!(state2, current_time = 1.0)
@test state2.x == zeros(2) && state2.current_time == 1.0
#or in the one changing x:
reinit!(state2, ones(2), current_time = 1.0)
@test state2.x == ones(2) && state2.current_time == 1.0

###############################################################################
#The State has also a private function guaranteeing there are no NaN
OK = Stopping._domain_check(state1) #function returns a boolean
@test OK == false #no NaN

@test Stopping._domain_check(state2) == false
update!(state2, x=[NaN, 0.0])
@test Stopping._domain_check(state2) == true

end


text:
@static if VERSION < v"0.7-"
    struct Adjoint{T,A} <: AbstractMatrix{T}
        parent::A
    end
else
    import LinearAlgebra: Adjoint
end

abstract type AbstractRecursiveMatrix{T,M,N,L} <: StaticArrays.StaticArray{Tuple{M,N}, T, 2} end
abstract type MutableRecursiveMatrix{T,M,N,L} <: AbstractRecursiveMatrix{T,M,N,L} end

struct StaticRecursiveMatrix{T,M,N,L} <: AbstractRecursiveMatrix{T,M,N,L}
    data::NTuple{L,T}
    StaticRecursiveMatrix{T,M,N,L}(data::NTuple{L,T}) where {T,M,N,L} = new(data)
end
mutable struct RecursiveMatrix{T,M,N,L} <: MutableRecursiveMatrix{T,M,N,L}
    data::NTuple{L,T}
    RecursiveMatrix{T,M,N,L}(data::NTuple{L,T}) where {T,M,N,L} = new(data)
    RecursiveMatrix{T,M,N,L}() where {T,M,N,L} = new()
    @generated RecursiveMatrix{T,M,N}() where {T,M,N} = :(RecursiveMatrix{$T,$M,$N,$(M*N)}())
end

function RecursiveMatrix(data::NTuple{L,T}, ::Val{M}, ::Val{N}) where {M,N,L,T}
    RecursiveMatrix{T,M,N,L}(data)
end
function RecursiveMatrix(::Type{T}, ::Val{M}, ::Val{N}, ::Val{L}) where {M,N,L,T}
    RecursiveMatrix{T,M,N,L}()
end
function RecursiveMatrix(::Type{T}, ::Val{M}, ::Val{N}) where {T,M,N}
    RecursiveMatrix(T,Val{M}(),Val{N}(),ValProd(Val{M}(),Val{L}()))
end
function RecursiveMatrix(data::NTuple{L,T}, ::Val{M}) where {M,L,T}
    RecursiveMatrix(data, Val{M}(), ValDiv(Val{L}(),Val{M}()))
end

function randmat(::Val{M},::Val{N},::Val{L}) where {M,N,L}
    out = RecursiveMatrix{Float64,M,N,L}()
    @inbounds for i = 1:L
        out[i] = randn()
    end
    out
end
randmat(m,n) = randmat(Val(m),Val(n),Val(m*n))
randmat(n) = randmat(n,n)
srandmat(m,n) = StaticRecursiveMatrix{Float64,m,n,m*n}(ntuple(i -> randn(), Val(m*n)))
srandmat(n) = srandmat(n,n)

"""
Creating a recursive matrix without compile time size information is not type stable.
"""
function RecursiveMatrix(data::AbstractMatrix)
    M, N = size(data)
    RecursiveMatrix(data, Val{M}(), Val{N}(), Val{M*N}())
end
function RecursiveMatrix(data::AbstractMatrix, ::Val{M}, ::Val{N}) where {M,N}
    RecursiveMatrix(data, Val{M}(), Val{N}(), ValProd(Val{M}(), Val{N}()))
end
function RecursiveMatrix(data::AbstractMatrix{T}, ::Val{M}, ::Val{N}, ::Val{L}) where {T,M,N,L}
    out = RecursiveMatrix(T, Val{M}(), Val{N}() ,Val{L}())
    # Not an efficient order of indexing into the recursive matrix.
    # Optimize later by doing it in branch-less chunks.
    for i ∈ 1:N, j ∈ 1:M
        out[j,i] = data[j,i]
    end
    out
end

"""
The methods in this library assume that operations on the pointer matrices do not alias.
It is discouraged to use these externally at all. If you do, make sure you don't violate these assumptions.

For this reason, few methods are actually implemented on this type. It is used internally via direct llvmcalls.
"""
struct PointerRecursiveMatrix{T,M,N,L} <: MutableRecursiveMatrix{T,M,N,L}
    data::Ptr{StaticRecursiveMatrix{T,cutoff,cutoff,cutoff2}}
    # PointerRecursiveMatrix(ptr::Ptr{StaticRecursiveMatrix{T,cutoff,cutoff,cutoff2}}, ::Val{M}, ::Val{N}, ::Val{L}) where {T,M,N,L} = new{T,M,N,L}(ptr)
end
function PointerRecursiveMatrix(parent::RecursiveMatrix{T,MP,NP,LP}, ::Val{M}, ::Val{N}, offset=0) where {T,M,N,MP,NP,LP}
    PointerRecursiveMatrix(
        Base.unsafe_convert(Ptr{StaticRecursiveMatrix{T,cutoff,cutoff,cutoff2}}, pointer_from_objref(parent) + offset * sizeof(T)),
        Val{M}(), Val{N}(), ValProd(Val{M}(), Val{N}())
    )
end

# Needs to be made faster.
function Base.Matrix(A::RecursiveMatrix{T,M,N,L}) where {T,M,N,L}
    out = Matrix{T}(undef, M, N)
    @inbounds for i = 1:N, j = 1:M
        out[j,i] = A[j,i]
    end
    out
end

struct Block{m,n} end
struct HalfBlock{m,n,h} end
Base.@pure Block(m,n) = Block{m,n}()
Base.@pure HalfBlock(m,n,h) = HalfBlock{m,n,h}()
struct BlockIndex{T,M,N,L}
    i::UInt
end

const RecursiveMatrixOrTranpose{T,M,N,L} = Union{
    RecursiveMatrix{T,M,N,L},
    Adjoint{T,RecursiveMatrix{T,N,M,L}}
}
const RecursivePointerMatrixOrTranpose{T,M,N,L} = Union{
    PointerRecursiveMatrix{T,M,N,L},
    Adjoint{T,PointerRecursiveMatrix{T,N,M,L}}
}
const RecurseOrTranpose{T,M,N,L} = Union{
    RecursiveMatrix{T,M,N,L},
    Adjoint{T,RecursiveMatrix{T,N,M,L}},
    StaticRecursiveMatrix{T,M,N,L},
    Adjoint{T,StaticRecursiveMatrix{T,M,N,L}},
    PointerRecursiveMatrix{T,M,N,L},
    Adjoint{T,PointerRecursiveMatrix{T,N,M,L}}
}
const MutableRecurseOrTranpose{T,M,N,L} = Union{
    RecursiveMatrix{T,M,N,L},
    Adjoint{T,RecursiveMatrix{T,N,M,L}},
    PointerRecursiveMatrix{T,M,N,L},
    Adjoint{T,PointerRecursiveMatrix{T,N,M,L}}
}
function dereftype(::Type{PointerRecursiveMatrix{T,M,N,L}}) where {T,M,N,L}
    RecursiveMatrix{T,M,N,L}
end
function dereftype(::Type{Adjoint{T,PointerRecursiveMatrix{T,N,M,L}}}) where {T,M,N,L}
    Adjoint{T,RecursiveMatrix{T,N,M,L}}
end

# abstract type transposition end
# struct t <: transposition end
# struct n <: transposition end
# tranpose(::t) = n()
# tranpose(::n) = t()
istransposed(::AbstractRecursiveMatrix) = false #n()
istransposed(::Type{<:AbstractRecursiveMatrix}) = false #n()
istransposed(::Adjoint{T,<:AbstractRecursiveMatrix{T}}) where T = true #t()
istransposed(::Type{<:Adjoint{T,<:AbstractRecursiveMatrix{T}}}) where T = true #t()
vistransposed(::AbstractRecursiveMatrix) = Val{false}() #n()
vistransposed(::Type{<:AbstractRecursiveMatrix}) = Val{false}() #n()
vistransposed(::Adjoint{T,<:AbstractRecursiveMatrix{T}}) where T = Val{true}() #t()
vistransposed(::Type{<:Adjoint{T,<:AbstractRecursiveMatrix{T}}}) where T = Val{true}() #t()

sub2ind(t, dims, i, j) = t ? j + (i-1)*dims[2] : i + (j-1)*dims[1]

Base.size(::RecurseOrTranpose{T,M,N}) where {T,M,N} = (M,N)
# @static if VERSION < v"0.7-"
@inline Base.getindex(A::Adjoint{T,<:MutableRecursiveMatrix{T}}, i) where T = A.parent[i]
@inline Base.setindex!(A::Adjoint{T,<:MutableRecursiveMatrix{T}}, v, i) where T = (A.parent[i] = v)
# end

# function forloopmul!(C, A, B)
#     @inbounds for i = 1:size(B,2)
#         for j = 1:size(A,1)
#             t = A[j,1] * B[1,j]
#             for k = 1:size(A,2)
#                 t += A[j,k] * B[k,j]
#             end
#             C[j,i] = t
#         end
#     end
# end

@inline float_point(A::RecursiveMatrix{T}) where T = Base.unsafe_convert(Ptr{T}, pointer_from_objref(A))
@inline float_point(A::PointerRecursiveMatrix{T}) where T = Base.unsafe_convert(Ptr{T}, A.data)
@inline point(A::RecursiveMatrix) = pointer_from_objref(A)
@inline point(A::PointerRecursiveMatrix) = A.data
@inline convert_point(A::RecursiveMatrix, ::Type{T}) where T = Base.unsafe_convert(Ptr{T}, pointer_from_objref(A))
@inline convert_point(A::RecursiveMatrix{T}) where T = Base.unsafe_convert(Ptr{StaticRecursiveMatrix{T,cutoff,cutoff,cutoff2}}, pointer_from_objref(A))
@inline convert_point(A::PointerRecursiveMatrix, ::Type{T}) where T = Base.unsafe_convert(Ptr{T}, A.data)
@inline convert_point(A::PointerRecursiveMatrix) = A.data
@inline point(A::Adjoint{T,<:MutableRecurseOrTranpose{T}}) where T = point(A.parent)

@generated function Base.getindex(A::RecursiveMatrixOrTranpose{T,M,N,L}, i::Int) where {T,M,N,L}
    bounds_error_string = "($M, $N) array at index (\$i,\$j)."

    quote #Add recursion.
        Base.@_inline_meta
        @boundscheck begin
            $L < i && throw(BoundsError($bounds_error_string))
        end
        A.data[i]
        # A.data[][i]
    end
end
# @generated function Base.getindex(A::RecursiveMatrixOrTranpose{T,M,N,L}, i::Int, ::Type{K}) where {T,M,N,L,K}
#     bounds_error_string = "($M, $N) array at index (\$i,\$j)."
#     quote #Add recursion.
#         Base.@_inline_meta
#         # @boundscheck begin
#         #     $L < i && throw(BoundsError($bounds_error_string))
#         # end
#         unsafe_load(Base.unsafe_convert(Ptr{$K}, pointer_from_data(A.data) + i*$(sizeof(T)) ), 1 )
#         # A.data[][i]
#     end
# end
@inline function Base.getindex(A::StaticRecursiveMatrix{T,M,N,L}, i::Int) where {T,M,N,L}
    @boundscheck i > L && throw(BoundsError())
    @inbounds out = A.data[i]
    out
end
@inline function Base.getindex(A::StaticRecursiveMatrix{T,M,N,L}, i::Int, j::Int) where {T,M,N,L}
    @boundscheck ( i > M || j > N ) && throw(BoundsError())
    @inbounds out = A.data[ (j-1)*M + i ]
    out
end

@inline function Base.setindex!(S::Symmetric{T,RecursiveMatrix{T,N,N,L}}, val, i::Integer, j::Integer) where {T,N,L}
    S.data[i,j] = val
end
@inline function Base.setindex!(S::Symmetric{T,RecursiveMatrix{T,N,N,L}}, val, i::Integer) where {T,N,L}
    S.data[i] = val
end

# Base.getindex(A::RecursiveMatrixOrTranpose{T,M,N}, i::Int, j::Int) where {T,M,N} = A.data[dense_sub2ind((M,N),i,j)]
@generated function Base.getindex(A::RecursiveMatrixOrTranpose{T,M,N,L}, i::Int, j::Int) where {T,M,N,L}
    bounds_error_string = "($M, $N) array at index (\$i,\$j)."
    if istransposed(A)
        linear_ind_expr = :(dense_sub2ind(Val{$M}(), Val{$N}(), j, i))
    else
        linear_ind_expr = :(dense_sub2ind(Val{$M}(), Val{$N}(), i, j))
    end
    quote
        # Base.@_inline_meta
        @boundscheck ($M < i || $N < j) && throw(BoundsError($bounds_error_string))
        A.data[$linear_ind_expr]
    end
end
# function getindex_block_quote(::Type{T}, full_m_blocks, mdim, m_r, full_n_blocks, ndim, n_r, m, n) where T
#     if m <= full_m_blocks && n <= full_n_blocks
#         Lout = mdim*ndim
#         if mdim == ndim == cutoff
#             q = quote
#                     Base.@_inline_meta
#                     Base.unsafe_load(convert_point(A) + $(sizeof(T)*( (m-1)*Lout + (n-1)*M*ndim) )   )
#                 end
#         else
#             q = quote
#                     Base.@_inline_meta
#                     Base.unsafe_load(convert_point(A,
#                         Ptr{StaticRecursiveMatrix{$T,$mdim,$ndim,$Lout}}) +
#                         $(sizeof(T)*( (m-1)*Lout + (n-1)*M*ndim) )   )
#             end
#         end

#     elseif m <= full_m_blocks # we're on the far right
#         Lout = mdim * n_r
#         q = quote
#                 Base.@_inline_meta
#                 Base.unsafe_load(convert_point( A,
#                     Ptr{StaticRecursiveMatrix{$T,$mdim,$n_r,$Lout}}) +
#                     $(sizeof(T)*( (m-1)*Lout + (n-1)*M*ndim) )  )
#             end

#     elseif n <= full_n_blocks
#         Lout = m_r * ndim
#         q = quote
#                 Base.@_inline_meta
#                 Base.unsafe_load(convert_point( A,
#                     Ptr{StaticRecursiveMatrix{$T,$m_r,$ndim,$Lout}}) +
#                     $(sizeof(T)*( (m-1)*ndim*mdim + (n-1)*M*ndim) ) )
#             end

#     else
#         Lout = m_r*n_r
#         q = quote
#                 Base.@_inline_meta
#                 Base.unsafe_load(convert_point( A,
#                     Ptr{StaticRecursiveMatrix{$T,$m_r,$n_r,$Lout}}) +
#                     $(sizeof(T)*( L - Lout ))  )
#             end

#     end
#     q
# end

# @generated function Base.getindex(A::MutableRecursiveMatrix{T,M,N,L}, ::Block{m,n}) where {T,M,N,L,m,n}
#     if max(M,N) <= 3cutoff # We are returning a StaticRecursiveMatrix
#         full_m_blocks, incomplete_m_block = divrem(M, cutoff)
#         full_n_blocks, incomplete_n_block = divrem(N, cutoff)

#         q = getindex_block_quote(T, full_m_blocks, cutoff, incomplete_m_block, full_n_blocks, cutoff, incomplete_n_block, m, n)

#     else # We are returning a pointer matrix.

#         if reasonably_square(minmax(M,N)...) # split both dims
#             M1, M2 = split_dim(M)
#             N1, N2 = split_dim(N)

#             Mout = (M1, M2)[m]
#             Nout = (N1, N2)[n]
#             if n == 1
#                 offset = sizeof(T) * (m-1)*M1*N1
#             else # n == 2
#                 offset = sizeof(T) * ( (m-1)*M1*N2 + M*N1 )
#             end
#             q = :( PointerRecursiveMatrix{$T,$Mout,$Nout,$(Mout*Nout)}( convert_point(A) + $offset ) )

#         elseif M > N
#             M1, M2 = split_dim(M)

#             Mout = (M1, M2)[m]
#             offset = sizeof(T) * (m-1)*M1*N
#             q = :( PointerRecursiveMatrix{$T,$Mout,$N,$(Mout*N)}( convert_point(A) + $offset ) )

#         else # M < N
#             N1, N2 = split_dim(N)

#             Nout = (N1, N2)[n]
#             offset = sizeof(T) * (n-1)*M*N1
#             q = :( PointerRecursiveMatrix{$T,$M,$Nout,$(M*Nout)}( convert_point(A) + $offset ) )

#         end

#     end
#     q
# end

# @generated function Base.getindex(A::MutableRecursiveMatrix{T,M,N,L}, ::HalfBlock{m,n}) where {T,M,N,L,m,n}
#     full_m_blocks, incomplete_m_block = divrem(M, cutoff)
#     full_n_blocks, incomplete_n_block = divrem(N, halfcutoff)
#     getindex_block_quote(T, full_m_blocks, cutoff, incomplete_m_block, full_n_blocks, halfcutoff, incomplete_n_block, m, n)
# end

@generated function Base.setindex!(A::RecursiveMatrixOrTranpose{T,M,N,L}, val::T, i::Int) where {T,M,N,L}
    bounds_error_string = "($M, $N) array at index (\$i,\$j)."
    quote #Add recursion.
        Base.@_inline_meta
        @boundscheck $L < i && throw(BoundsError($bounds_error_string))
        unsafe_store!(float_point(A), val, i )
    end
end
@generated function Base.setindex!(A::RecursiveMatrixOrTranpose{T,M,N,L}, val::I, i::Int) where {T,M,N,L,I<:Integer}
    bounds_error_string = "($M, $N) array at index (\$i,\$j)."
    quote #Add recursion.
        # Base.@_inline_meta
        @boundscheck begin
            $L < i && throw(BoundsError($bounds_error_string))
        end
        unsafe_store!(float_point(A), convert(T,val), i )
    end
end
@generated function Base.setindex!(A::RecursiveMatrixOrTranpose{T,M,N,L}, val::V, i::Int) where {T,M,N,L,V}
    bounds_error_string = "($M, $N) array at index (\$i,\$j)."
    quote #Add recursion.
        Base.@_inline_meta
        @boundscheck begin
            $L < i && throw(BoundsError($bounds_error_string))
        end
        unsafe_store!(Base.unsafe_convert(Ptr{V}, pointer_from_objref(A) + sizeof(T)*(i-1) ), val)#, 1 )
    end
end

@generated function Base.setindex!(A::RecursiveMatrixOrTranpose{T,M,N,L}, val::T, i::Int, j::Int) where {T,M,N,L}
    bounds_error_string = "($M, $N) array at index (\$i,\$j)."
    if istransposed(A)
        linear_ind_expr = :(dense_sub2ind(Val{$M}(), Val{$N}(), j, i))
    else
        linear_ind_expr = :(dense_sub2ind(Val{$M}(), Val{$N}(), i, j))
    end
    quote #Add recursion.
        # Base.@_inline_meta
        @boundscheck begin
            ($M < i || $N < j) && throw(BoundsError($bounds_error_string))
        end
        # unsafe_store!(point(A), convert(T, val), sub2ind(istransposed(A), ($M,$N), i, j) )

        unsafe_store!(float_point(A), val, $linear_ind_expr )
    end
end

@generated function Base.setindex!(A::RecursiveMatrixOrTranpose{T,M,N,L}, val::V, i::Int, j::Int) where {T,M,N,L,V}
    bounds_error_string = "($M, $N) array at index (\$i,\$j)."
    if istransposed(A)
        linear_ind_expr = :(dense_sub2ind(Val{$M}(), Val{$N}(), j, i))
    else
        linear_ind_expr = :(dense_sub2ind(Val{$M}(), Val{$N}(), i, j))
    end
    quote #Add recursion.
        # Base.@_inline_meta
        @boundscheck begin
            ($M < i || $N < j) && throw(BoundsError($bounds_error_string))
        end
        # unsafe_store!(point(A), convert(T, val), sub2ind(istransposed(A), ($M,$N), i, j) )
        unsafe_store!(Base.unsafe_convert(Ptr{V}, pointer_from_objref(A) + sizeof(T)*(($linear_ind_expr)-1) ), val)#, 1 )
    end
end

# """
# The second block indicates that the number of dimensions we are splitting things into.
# """
# @generated function BlockIndex(::Type{MutableRecursiveMatrix{T,M,N,L}}, ::Block{m,n}, ::Block{MC,NC}) where {T,M,N,L,m,n,MC,NC}
#     full_m_blocks, m_r = divrem(M, cutoff)
#     full_n_blocks, n_r = divrem(N, cutoff)

#     if m <= full_m_blocks && n <= full_n_blocks
#         return BlockIndex{T,cutoff,cutoff,cutoff2}(sizeof(T)*( (m-1)*cutoff2 + (n-1)*M*cutoff) )
#     elseif m <= full_m_blocks # we're on the far right
#         return BlockIndex{T,cutoff,n_r,cutoff*n_r}(sizeof(T)*( (m-1)*cutoff*n_r + (n-1)*M*cutoff))
#     elseif n <= full_n_blocks
#         return BlockIndex{T,m_r,cutoff,m_r*cutoff}(sizeof(T)*( (m-1)*cutoff2 + (n-1)*M*cutoff ) )
#     else
#         return BlockIndex{T,m_r,n_r,m_r*n_r}(sizeof(T)*( L - m_r*n_r ) )
#     end
# end
#
# This function shouldn't really be typed on m,n, or...
#
@generated function BlockIndex(::Type{RecursiveMatrix{T,M,N,L}}, ::Block{m,n}) where {T,M,N,L,m,n}
    full_n_blocks, n_r = divrem(N, cutoff)
    full_m_blocks, m_r = divrem(M, cutoff)

    offset = sizeof(T) * ( dense_sub2ind((M,N), (m-1)*cutoff + 1, (n-1)*cutoff + 1 ) - 1 )

    if m <= full_m_blocks && n <= full_n_blocks
        return BlockIndex{T,cutoff,cutoff,cutoff2}( offset )
    elseif m <= full_m_blocks # we're on the far right
        return BlockIndex{T,cutoff,n_r,cutoff*n_r}( offset )
    elseif n <= full_n_blocks
        return BlockIndex{T,m_r,cutoff,m_r*cutoff}( offset )
    else
        return BlockIndex{T,m_r,n_r,m_r*n_r}( offset )
    end
end

@generated function BlockIndex(::Type{RecursiveMatrix{T,M,N,L}}, ::HalfBlock{m,n,h}) where {T,M,N,L,m,n,h}
    full_n_blocks, n_r = divrem(N, cutoff)
    full_m_blocks, m_r = divrem(M, cutoff)

    offset = sizeof(T) * ( dense_sub2ind((M,N), (m-1)*cutoff + 1, (n-1)*cutoff + (h-1)*halfcutoff + 1 ) - 1 )

    if m <= full_m_blocks && n <= full_n_blocks
        return BlockIndex{T,cutoff,halfcutoff,cutoff*halfcutoff}( offset )
    elseif m <= full_m_blocks # we're on the far right
        ndim = h == 1 ? min(n_r, halfcutoff) : n_r - halfcutoff
        return BlockIndex{T,cutoff,ndim,cutoff*ndim}( offset )
    elseif n <= full_n_blocks # we're on the bottom
        return BlockIndex{T,m_r,halfcutoff,m_r*halfcutoff}( offset )
    else
        ndim = h == 1 ? min(n_r, halfcutoff) : n_r - halfcutoff
        return BlockIndex{T,m_r,ndim,m_r*ndim}( offset )
    end
end

@generated function Base.getindex(A::MutableRecursiveMatrix{T,M,N,L}, i::BlockIndex{T,m,n,mn}) where {T,M,N,L,m,n,mn}
    # size_check sizeof(T)*(L-mn)
    # SM = StaticRecursiveMatrix{T,m,n,mn}
    quote
        Base.@_inline_meta
        # @boundscheck i.i > $size_check && throw(BoundsError())
        # Base.unsafe_load(convert_point(A, $SM ) + i.i)
        @boundscheck i.i > $(sizeof(T)*(L-mn))  && throw(BoundsError())
        Base.unsafe_load(convert_point(A, StaticRecursiveMatrix{$T,$m,$n,$mn} ) + i.i)
    end
end
@inline function Base.getindex(A::MutableRecursiveMatrix{T,M,N,L}, i::BlockIndex{T,cutoff,cutoff,cutoff2}) where {T,M,N,L}
    @boundscheck i.i > (L-cutoff2)*sizeof(T) && throw(BoundsError())
    Base.unsafe_load(convert_point(A) + i.i)
end

"""
The setindex methods don't work too well...
"""
@generated function Base.setindex!(A::MutableRecursiveMatrix{T,M,N,L},
                                    val::StaticRecursiveMatrix{T,m,n,mn},
                                    i::BlockIndex{T,m,n,mn}) where {T,M,N,L,m,n,mn}
    quote
        Base.@_inline_meta
        @boundscheck i.i > $(sizeof(T)*(L-mn))  && throw(BoundsError())
        Base.unsafe_store!(convert_point(A, StaticRecursiveMatrix{$T,$m,$n,$mn} ) + i.i, val); nothing
    end
end

@inline function Base.setindex!(A::MutableRecursiveMatrix{T,M,N,L},
                                    val::StaticRecursiveMatrix{T,cutoff,cutoff,cutoff2},
                                    i::BlockIndex{T,cutoff,cutoff,cutoff2}) where {T,M,N,L}
    @boundscheck i.i > (L-cutoff2)*sizeof(T) && throw(BoundsError())
    Base.unsafe_store!(convert_point(A) + i.i, val); nothing
end

@generated function point(A::MutableRecursiveMatrix{T}, i::BlockIndex{T,m,n,mn}) where {T,m,n,mn}
    :(convert_point(A, StaticRecursiveMatrix{$T,$m,$n,$mn} ) + i.i)
end
function point(A::MutableRecursiveMatrix{T}, i::BlockIndex{T,cutoff,cutoff,cutoff2}) where T
    convert_point(A) + i.i
end

# function setindex_block_quote(::Type{T}, ::Type{V}, full_m_blocks, mdim, m_r, full_n_blocks, ndim, n_r, m, n) where {T,V}
#     if m <= full_m_blocks && n <= full_n_blocks
#         Lout = mdim*ndim
#         if mdim == ndim == cutoff
#             q = quote
#                     Base.@_inline_meta
#                     Base.unsafe_store!(convert_point(A) + $(sizeof(T)*( (m-1)*Lout + (n-1)*M*ndim)), val )
#                 end
#         else
#             q = quote
#                     Base.@_inline_meta
#                     Base.unsafe_store!(convert_point(A, Ptr{V}) + $(sizeof(T)*( (m-1)*Lout + (n-1)*M*ndim) ), val )
#             end
#         end

#     elseif m <= full_m_blocks # we're on the far right
#         Lout = mdim * n_r
#         q = quote
#                 Base.@_inline_meta
#                 Base.unsafe_store!(convert_point( A, Ptr{V}) + $(sizeof(T)*( (m-1)*Lout + (n-1)*M*ndim) ). val )
#             end

#     elseif n <= full_n_blocks
#         Lout = m_r * ndim
#         q = quote
#                 Base.@_inline_meta
#                 Base.unsafe_store!(convert_point( A, Ptr{V}) + $(sizeof(T)*( (m-1)*ndim*mdim + (n-1)*M*ndim) ), val )
#             end

#     else
#         Lout = m_r*n_r
#         q = quote
#                 Base.@_inline_meta
#                 Base.unsafe_store!(convert_point( A, Ptr{V}) + $(sizeof(T)*( L - Lout )), val )
#             end

#     end
#     q
# end

# @generated function Base.setindex!(A::MutableRecursiveMatrix{T,M,N,L}, val::V, ::HalfBlock{m,n}) where {T,M,N,L,V,m,n}
#     full_m_blocks, incomplete_m_block = divrem(M, cutoff)
#     full_n_blocks, incomplete_n_block = divrem(N, halfcutoff)
#     setindex_block_quote(T, V, full_m_blocks, cutoff, incomplete_m_block, full_n_blocks, halfcutoff, incomplete_n_block, m, n)
# end

# @generated function Base.setindex!(A::MutableRecursiveMatrix{T,M,N,L}, val::V, ::Block{m,n}) where {T,M,N,L,V,m,n}
#     full_m_blocks, incomplete_m_block = divrem(M, cutoff)
#     full_n_blocks, incomplete_n_block = divrem(N, cutoff)
#     setindex_block_quote(T, V, full_m_blocks, cutoff, incomplete_m_block, full_n_blocks, cutoff, incomplete_n_block, m, n)
# end

# rmatA8 = RecursiveMatrix{Float64,8,8,abs2(8)}(Ref(ntuple(i -> randn(), Val(abs2(8)))));
# rmatB8 = RecursiveMatrix{Float64,8,8,abs2(8)}(Ref(ntuple(i -> randn(), Val(abs2(8)))));
# rmatC8 = RecursiveMatrix{Float64,8,8,abs2(8)}(Ref{NTuple{64,Float64}}());
# rmatA12 = RecursiveMatrix{Float64,12,12,abs2(12)}(Ref(ntuple(i -> randn(), Val(abs2(12)))));
# rmatB12 = RecursiveMatrix{Float64,12,12,abs2(12)}(Ref(ntuple(i -> randn(), Val(abs2(12)))));
# rmatC12 = RecursiveMatrix{Float64,12,12,abs2(12)}(Ref{NTuple{144,Float64}}());
# rmatA16 = RecursiveMatrix{Float64,16,16,abs2(16)}(Ref(ntuple(i -> randn(), Val(abs2(16)))));
# rmatB16 = RecursiveMatrix{Float64,16,16,abs2(16)}(Ref(ntuple(i -> randn(), Val(abs2(16)))));
# rmatC16 = RecursiveMatrix{Float64,16,16,abs2(16)}(Ref{NTuple{256,Float64}}());