text:
using ArrayInterface, Test
using Base: setindex
import ArrayInterface: has_sparsestruct, findstructralnz, fast_scalar_indexing, lu_instance, device, contiguous_axis, contiguous_batch_size, stride_rank, dense_dims, StaticInt
@test ArrayInterface.ismutable(rand(3))

using Aqua
Aqua.test_all(ArrayInterface)

using StaticArrays
x = @SVector [1,2,3]
@test ArrayInterface.ismutable(x) == false
@test ArrayInterface.ismutable(view(x, 1:2)) == false
x = @MVector [1,2,3]
@test ArrayInterface.ismutable(x) == true
@test ArrayInterface.ismutable(view(x, 1:2)) == true
@test ArrayInterface.ismutable(1:10) == false
@test ArrayInterface.ismutable((0.1,1.0)) == false
@test ArrayInterface.ismutable(Base.ImmutableDict{Symbol,Int64}) == false
@test ArrayInterface.ismutable((;x=1)) == false

@test isone(ArrayInterface.known_first(typeof(StaticArrays.SOneTo(7))))
@test ArrayInterface.known_last(typeof(StaticArrays.SOneTo(7))) == 7
@test ArrayInterface.known_length(typeof(StaticArrays.SOneTo(7))) == 7

using LinearAlgebra, SparseArrays

D=Diagonal([1,2,3,4])
@test has_sparsestruct(D)
rowind,colind=findstructralnz(D)
@test [D[rowind[i],colind[i]] for i in 1:length(rowind)]==[1,2,3,4]
@test length(rowind)==4
@test length(rowind)==length(colind)

Bu = Bidiagonal([1,2,3,4], [7,8,9], :U)
@test has_sparsestruct(Bu)
rowind,colind=findstructralnz(Bu)
@test [Bu[rowind[i],colind[i]] for i in 1:length(rowind)]==[1,7,2,8,3,9,4]
Bl = Bidiagonal([1,2,3,4], [7,8,9], :L)
@test has_sparsestruct(Bl)
rowind,colind=findstructralnz(Bl)
@test [Bl[rowind[i],colind[i]] for i in 1:length(rowind)]==[1,7,2,8,3,9,4]

Tri=Tridiagonal([1,2,3],[1,2,3,4],[4,5,6])
@test has_sparsestruct(Tri)
rowind,colind=findstructralnz(Tri)
@test [Tri[rowind[i],colind[i]] for i in 1:length(rowind)]==[1,2,3,4,4,5,6,1,2,3]
STri=SymTridiagonal([1,2,3,4],[5,6,7])
@test has_sparsestruct(STri)
rowind,colind=findstructralnz(STri)
@test [STri[rowind[i],colind[i]] for i in 1:length(rowind)]==[1,2,3,4,5,6,7,5,6,7]

Sp=sparse([1,2,3],[1,2,3],[1,2,3])
@test has_sparsestruct(Sp)
rowind,colind=findstructralnz(Sp)
@test [Tri[rowind[i],colind[i]] for i in 1:length(rowind)]==[1,2,3]

@test ArrayInterface.ismutable(spzeros(1, 1))
@test ArrayInterface.ismutable(spzeros(1))

@test !fast_scalar_indexing(qr(rand(10, 10)).Q)
@test !fast_scalar_indexing(qr(rand(10, 10), Val(true)).Q)
@test !fast_scalar_indexing(lq(rand(10, 10)).Q)

using BandedMatrices

B=BandedMatrix(Ones(5,5), (-1,2))
B[band(1)].=[1,2,3,4]
B[band(2)].=[5,6,7]
@test has_sparsestruct(B)
rowind,colind=findstructralnz(B)
@test [B[rowind[i],colind[i]] for i in 1:length(rowind)]==[5,6,7,1,2,3,4]
B=BandedMatrix(Ones(4,6), (-1,2))
B[band(1)].=[1,2,3,4]
B[band(2)].=[5,6,7,8]
rowind,colind=findstructralnz(B)
@test [B[rowind[i],colind[i]] for i in 1:length(rowind)]==[5,6,7,8,1,2,3,4]

using BlockBandedMatrices
BB=BlockBandedMatrix(Ones(10,10),[1,2,3,4],[4,3,2,1],(1,0))
BB[Block(1,1)].=[1 2 3 4]
BB[Block(2,1)].=[5 6 7 8;9 10 11 12]
rowind,colind=findstructralnz(BB)
@test [BB[rowind[i],colind[i]] for i in 1:length(rowind)]==
    [1,5,9,2,6,10,3,7,11,4,8,12,
    1,1,1,1,1,1,1,1,1,1,1,1,1,1,
    1,1,1,1,1,1,1,1,1,1,1,1,1,1,
    1,1,1,1,1]

dense=collect(Ones(8,8))
for i in 1:8
    dense[:,i].=[1,2,3,4,5,6,7,8]
end
BBB=BandedBlockBandedMatrix(dense, [4, 4] ,[4, 4], (1, 1), (1, 1))
rowind,colind=findstructralnz(BBB)
@test [BBB[rowind[i],colind[i]] for i in 1:length(rowind)]==
    [1,2,3,1,2,3,4,2,3,4,5,6,7,5,6,7,8,6,7,8,
     1,2,3,1,2,3,4,2,3,4,5,6,7,5,6,7,8,6,7,8]

@testset "setindex" begin
    @testset "$(typeof(x))" for x in [
        zeros(3),
        falses(3),
        spzeros(3),
    ]
        y = setindex(x, true, 1)
        @test iszero(x)  # x is not mutated
        @test y[1] == true
        @test iszero(x[CartesianIndices(size(x)) .== [CartesianIndex(1)]])

        y2 = setindex(x, one.(x), :)
        @test iszero(x)
        @test all(isone, y2)
    end

    @testset "$(typeof(x))" for x in [
        zeros(3, 3),
        falses(3, 3),
        spzeros(3, 3),
    ]
        y = setindex(x, true, 1, 1)
        @test iszero(x)  # x is not mutated
        @test y[1, 1] == true
        @test iszero(x[CartesianIndices(size(x)) .== [CartesianIndex(1, 1)]])

        y2 = setindex(x, one.(x), :, :)
        @test iszero(x)
        @test all(isone, y2)
    end

    @testset "$(typeof(x))" for x in [
        zeros(3, 3, 3),
        falses(3, 3, 3),
    ]
        y = setindex(x, true, 1, 1, 1)
        @test iszero(x)  # x is not mutated
        @test y[1, 1, 1] == true
        @test iszero(x[CartesianIndices(size(x)) .== [CartesianIndex(1, 1, 1)]])

        y2 = setindex(x, one.(x), :, :, :)
        @test iszero(x)
        @test all(isone, y2)
    end
end

using SuiteSparse
@testset "lu_instance" begin
  for A in [
    randn(5, 5),
    @SMatrix(randn(5, 5)),
    @MMatrix(randn(5, 5)),
    sprand(50, 50, 0.5)
  ]
    @test lu_instance(A) isa typeof(lu(A))
  end
  @test lu_instance(1) === 1
end

using Random
using ArrayInterface: issingular
@testset "issingular" begin
    for T in [Float64, ComplexF64]
        R = randn(MersenneTwister(2), T, 5, 5)
        S = Symmetric(R)
        L = UpperTriangular(R)
        U = LowerTriangular(R)
        @test all(!issingular, [R, S, L, U, U'])
        R[:, 2] .= 0
        @test all(issingular, [R, L, U, U'])
        @test !issingular(S)
        R[2, :] .= 0
        @test issingular(S)
        @test all(!issingular, [UnitLowerTriangular(R), UnitUpperTriangular(R), UnitUpperTriangular(R)'])
    end
end

using ArrayInterface: zeromatrix
@test zeromatrix(rand(4,4,4)) == zeros(4*4*4,4*4*4)

using ArrayInterface: parent_type
@testset "Parent Type" begin
    x = ones(4, 4)
    @test parent_type(view(x, 1:2, 1:2)) <: typeof(x)
    @test parent_type(reshape(x, 2, :)) <: typeof(x)
    @test parent_type(transpose(x)) <: typeof(x)
    @test parent_type(Symmetric(x)) <: typeof(x)
    @test parent_type(UpperTriangular(x)) <: typeof(x)
    @test parent_type(PermutedDimsArray(x, (2,1))) <: typeof(x)
    @test parent_type(Base.Slice(1:10)) <: UnitRange{Int}
end

@testset "Range Interface" begin
    @testset "Range Constructors" begin
        @test @inferred(StaticInt(1):StaticInt(10)) == 1:10
        @test @inferred(StaticInt(1):StaticInt(2):StaticInt(10)) == 1:2:10
        @test @inferred(1:StaticInt(2):StaticInt(10)) == 1:2:10
        @test @inferred(StaticInt(1):StaticInt(2):10) == 1:2:10
        @test @inferred(StaticInt(1):2:StaticInt(10)) == 1:2:10
        @test @inferred(1:2:StaticInt(10)) == 1:2:10
        @test @inferred(1:StaticInt(2):10) == 1:2:10
        @test @inferred(StaticInt(1):2:10) == 1:2:10

        @test @inferred(StaticInt(1):StaticInt(1):StaticInt(10)) === ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), StaticInt(10))
        @test @inferred(StaticInt(1):StaticInt(1):10) === ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), 10)
        @test @inferred(1:StaticInt(1):10) === ArrayInterface.OptionallyStaticUnitRange(1, 10)
        @test length(StaticInt{-1}():StaticInt{-1}():StaticInt{-10}()) == 10

        @test UnitRange(ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), StaticInt(10))) === UnitRange(1, 10)
        @test UnitRange{Int}(ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), StaticInt(10))) === UnitRange(1, 10)

        @test AbstractUnitRange{Int}(ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), StaticInt(10))) isa ArrayInterface.OptionallyStaticUnitRange
        @test AbstractUnitRange{UInt}(ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), StaticInt(10))) isa Base.OneTo
        @test AbstractUnitRange{UInt}(ArrayInterface.OptionallyStaticUnitRange(StaticInt(2), StaticInt(10))) isa UnitRange

        @test @inferred((StaticInt(1):StaticInt(10))[StaticInt(2):StaticInt(3)]) === StaticInt(2):StaticInt(3)
        @test @inferred((StaticInt(1):StaticInt(10))[StaticInt(2):3]) === StaticInt(2):3
        @test @inferred((StaticInt(1):StaticInt(10))[2:3]) === 2:3
        @test @inferred((1:StaticInt(10))[StaticInt(2):StaticInt(3)]) === 2:3

        @test -(StaticInt{1}():StaticInt{10}()) === StaticInt{-1}():StaticInt{-1}():StaticInt{-10}()

        @test reverse(StaticInt{1}():StaticInt{10}()) === StaticInt{10}():StaticInt{-1}():StaticInt{1}()
        @test reverse(StaticInt{1}():StaticInt{2}():StaticInt{9}()) === StaticInt{9}():StaticInt{-2}():StaticInt{1}()
    end

    @test isnothing(@inferred(ArrayInterface.known_first(typeof(1:4))))
    @test isone(@inferred(ArrayInterface.known_first(Base.OneTo(4))))
    @test isone(@inferred(ArrayInterface.known_first(typeof(Base.OneTo(4)))))

    @test isnothing(@inferred(ArrayInterface.known_last(1:4)))
    @test isnothing(@inferred(ArrayInterface.known_last(typeof(1:4))))

    @test isnothing(@inferred(ArrayInterface.known_step(typeof(1:0.2:4))))
    @test isone(@inferred(ArrayInterface.known_step(1:4)))
    @test isone(@inferred(ArrayInterface.known_step(typeof(1:4))))

    @testset "length" begin
        @test @inferred(length(ArrayInterface.OptionallyStaticUnitRange(1, 0))) == 0
        @test @inferred(length(ArrayInterface.OptionallyStaticUnitRange(1, 10))) == 10
        @test @inferred(length(ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), 10))) == 10
        @test @inferred(length(ArrayInterface.OptionallyStaticUnitRange(StaticInt(0), 10))) == 11
        @test @inferred(length(ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), StaticInt(10)))) == 10
        @test @inferred(length(ArrayInterface.OptionallyStaticUnitRange(StaticInt(0), StaticInt(10)))) == 11

        @test @inferred(length(StaticInt(1):StaticInt(2):StaticInt(0))) == 0
        @test @inferred(length(StaticInt(0):StaticInt(-2):StaticInt(1))) == 0
    end

    @test @inferred(getindex(ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), 10), 1)) == 1
    @test @inferred(getindex(ArrayInterface.OptionallyStaticUnitRange(StaticInt(0), 10), 1)) == 0
    @test_throws BoundsError getindex(ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), 10), 0)
    @test_throws BoundsError getindex(ArrayInterface.OptionallyStaticStepRange(StaticInt(1), 2, 10), 0)
    @test_throws BoundsError getindex(ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), 10), 11)
    @test_throws BoundsError getindex(ArrayInterface.OptionallyStaticStepRange(StaticInt(1), 2, 10), 11)
end

@testset "Memory Layout" begin
    A = zeros(3,4,5);
    @test device(A) === ArrayInterface.CPUPointer()
    @test device((1,2,3)) === ArrayInterface.CPUIndex()
    @test device(PermutedDimsArray(A,(3,1,2))) === ArrayInterface.CPUPointer()
    @test device(view(A, 1, :, 2:4)) === ArrayInterface.CPUPointer()
    @test device(view(A, 1, :, 2:4)') === ArrayInterface.CPUPointer()
    @test device(@SArray(zeros(2,2,2))) === ArrayInterface.CPUIndex()
    @test device(@view(@SArray(zeros(2,2,2))[1,1:2,:])) === ArrayInterface.CPUIndex()
    @test device(@MArray(zeros(2,2,2))) === ArrayInterface.CPUPointer()
    @test isnothing(device("Hello, world!"))

    @test @inferred(contiguous_axis(@SArray(zeros(2,2,2)))) === ArrayInterface.Contiguous(1)
    @test @inferred(contiguous_axis(A)) === ArrayInterface.Contiguous(1)
    @test @inferred(contiguous_axis(PermutedDimsArray(A,(3,1,2)))) === ArrayInterface.Contiguous(2)
    @test @inferred(contiguous_axis(@view(PermutedDimsArray(A,(3,1,2))[2,1:2,:]))) === ArrayInterface.Contiguous(1)
    @test @inferred(contiguous_axis(transpose(@view(PermutedDimsArray(A,(3,1,2))[2,1:2,:])))) === ArrayInterface.Contiguous(2)
    @test @inferred(contiguous_axis(@view(PermutedDimsArray(A,(3,1,2))[2:3,1:2,:]))) === ArrayInterface.Contiguous(2)
    @test @inferred(contiguous_axis(@view(PermutedDimsArray(A,(3,1,2))[2:3,2,:]))) === ArrayInterface.Contiguous(-1)
    @test @inferred(contiguous_axis(@view(PermutedDimsArray(A,(3,1,2))[2:3,2,:])')) === ArrayInterface.Contiguous(-1)
    @test @inferred(contiguous_axis(@view(PermutedDimsArray(A,(3,1,2))[:,1:2,1])')) === ArrayInterface.Contiguous(1)

    @test @inferred(ArrayInterface.contiguous_axis_indicator(@SArray(zeros(2,2,2)))) === (Val(true),Val(false),Val(false))
    @test @inferred(ArrayInterface.contiguous_axis_indicator(A)) === (Val(true),Val(false),Val(false))
    @test @inferred(ArrayInterface.contiguous_axis_indicator(PermutedDimsArray(A,(3,1,2)))) === (Val(false),Val(true),Val(false))
    @test @inferred(ArrayInterface.contiguous_axis_indicator(@view(PermutedDimsArray(A,(3,1,2))[2,1:2,:]))) === (Val(true),Val(false))
    @test @inferred(ArrayInterface.contiguous_axis_indicator(@view(PermutedDimsArray(A,(3,1,2))[2,1:2,:])')) === (Val(false),Val(true))
    @test @inferred(ArrayInterface.contiguous_axis_indicator(@view(PermutedDimsArray(A,(3,1,2))[2:3,1:2,:]))) === (Val(false),Val(true),Val(false))
    @test @inferred(ArrayInterface.contiguous_axis_indicator(@view(PermutedDimsArray(A,(3,1,2))[2:3,2,:]))) === (Val(false),Val(false))
    @test @inferred(ArrayInterface.contiguous_axis_indicator(@view(PermutedDimsArray(A,(3,1,2))[2:3,2,:])')) === (Val(false),Val(false))
    @test @inferred(ArrayInterface.contiguous_axis_indicator(@view(PermutedDimsArray(A,(3,1,2))[:,1:2,1])')) === (Val(true),Val(false))

    @test @inferred(contiguous_batch_size(@SArray(zeros(2,2,2)))) === ArrayInterface.ContiguousBatch(0)
    @test @inferred(contiguous_batch_size(A)) === ArrayInterface.ContiguousBatch(0)
    @test @inferred(contiguous_batch_size(PermutedDimsArray(A,(3,1,2)))) === ArrayInterface.ContiguousBatch(0)
    @test @inferred(contiguous_batch_size(@view(PermutedDimsArray(A,(3,1,2))[2,1:2,:]))) === ArrayInterface.ContiguousBatch(0)
    @test @inferred(contiguous_batch_size(@view(PermutedDimsArray(A,(3,1,2))[2,1:2,:])')) === ArrayInterface.ContiguousBatch(0)
    @test @inferred(contiguous_batch_size(@view(PermutedDimsArray(A,(3,1,2))[2:3,1:2,:]))) === ArrayInterface.ContiguousBatch(0)
    @test @inferred(contiguous_batch_size(@view(PermutedDimsArray(A,(3,1,2))[2:3,2,:]))) === ArrayInterface.ContiguousBatch(-1)
    @test @inferred(contiguous_batch_size(@view(PermutedDimsArray(A,(3,1,2))[2:3,2,:])')) === ArrayInterface.ContiguousBatch(-1)
    @test @inferred(contiguous_batch_size(@view(PermutedDimsArray(A,(3,1,2))[:,1:2,1])')) === ArrayInterface.ContiguousBatch(0)

    @test @inferred(stride_rank(@SArray(zeros(2,2,2)))) === ArrayInterface.StrideRank((1, 2, 3))
    @test @inferred(stride_rank(A)) === ArrayInterface.StrideRank((1,2,3))
    @test @inferred(stride_rank(PermutedDimsArray(A,(3,1,2)))) === ArrayInterface.StrideRank((3, 1, 2))
    @test @inferred(stride_rank(@view(PermutedDimsArray(A,(3,1,2))[2,1:2,:]))) === ArrayInterface.StrideRank((1, 2))
    @test @inferred(stride_rank(@view(PermutedDimsArray(A,(3,1,2))[2,1:2,:])')) === ArrayInterface.StrideRank((2, 1))
    @test @inferred(stride_rank(@view(PermutedDimsArray(A,(3,1,2))[2:3,1:2,:]))) === ArrayInterface.StrideRank((3, 1, 2))
    @test @inferred(stride_rank(@view(PermutedDimsArray(A,(3,1,2))[2:3,2,:]))) === ArrayInterface.StrideRank((3, 2))
    @test @inferred(stride_rank(@view(PermutedDimsArray(A,(3,1,2))[2:3,2,:])')) === ArrayInterface.StrideRank((2, 3))
    @test @inferred(stride_rank(@view(PermutedDimsArray(A,(3,1,2))[:,1:2,1])')) === ArrayInterface.StrideRank((1, 3))
    @test @inferred(stride_rank(@view(PermutedDimsArray(A,(3,1,2))[:,2,1])')) === ArrayInterface.StrideRank((2, 1))

    @test @inferred(ArrayInterface.is_column_major(@SArray(zeros(2,2,2)))) === Val{true}()
    @test @inferred(ArrayInterface.is_column_major(A)) === Val{true}()
    @test @inferred(ArrayInterface.is_column_major(PermutedDimsArray(A,(3,1,2)))) === Val{false}()
    @test @inferred(ArrayInterface.is_column_major(@view(PermutedDimsArray(A,(3,1,2))[2,1:2,:]))) === Val{true}()
    @test @inferred(ArrayInterface.is_column_major(@view(PermutedDimsArray(A,(3,1,2))[2,1:2,:])')) === Val{false}()
    @test @inferred(ArrayInterface.is_column_major(@view(PermutedDimsArray(A,(3,1,2))[2:3,1:2,:]))) === Val{false}()
    @test @inferred(ArrayInterface.is_column_major(@view(PermutedDimsArray(A,(3,1,2))[2:3,2,:]))) === Val{false}()
    @test @inferred(ArrayInterface.is_column_major(@view(PermutedDimsArray(A,(3,1,2))[2:3,2,:])')) === Val{true}()
    @test @inferred(ArrayInterface.is_column_major(@view(PermutedDimsArray(A,(3,1,2))[:,1:2,1])')) === Val{true}()
    @test @inferred(ArrayInterface.is_column_major(@view(PermutedDimsArray(A,(3,1,2))[:,2,1])')) === Val{false}()
    @test @inferred(ArrayInterface.is_column_major(1:10)) === Val{false}()
    @test @inferred(ArrayInterface.is_column_major(2.3)) === Val{false}()

    @test @inferred(dense_dims(@SArray(zeros(2,2,2)))) === ArrayInterface.DenseDims((true,true,true))
    @test @inferred(dense_dims(A)) === ArrayInterface.DenseDims((true,true,true))
    @test @inferred(dense_dims(PermutedDimsArray(A,(3,1,2)))) === ArrayInterface.DenseDims((true,true,true))
    @test @inferred(dense_dims(@view(PermutedDimsArray(A,(3,1,2))[2,1:2,:]))) === ArrayInterface.DenseDims((true,false))
    @test @inferred(dense_dims(@view(PermutedDimsArray(A,(3,1,2))[2,1:2,:])')) === ArrayInterface.DenseDims((false,true))
    @test @inferred(dense_dims(@view(PermutedDimsArray(A,(3,1,2))[2:3,1:2,:]))) === ArrayInterface.DenseDims((false,true,false))
    @test @inferred(dense_dims(@view(PermutedDimsArray(A,(3,1,2))[2:3,:,1:2]))) === ArrayInterface.DenseDims((false,true,true))
    @test @inferred(dense_dims(@view(PermutedDimsArray(A,(3,1,2))[2:3,2,:]))) === ArrayInterface.DenseDims((false,false))
    @test @inferred(dense_dims(@view(PermutedDimsArray(A,(3,1,2))[2:3,2,:])')) === ArrayInterface.DenseDims((false,false))
    @test @inferred(dense_dims(@view(PermutedDimsArray(A,(3,1,2))[:,1:2,1])')) === ArrayInterface.DenseDims((true,false))

    B = Array{Int8}(undef, 2,2,2,2);
    doubleperm = PermutedDimsArray(PermutedDimsArray(B,(4,2,3,1)), (4,2,1,3));
    @test collect(strides(B))[collect(stride_rank(doubleperm))] == collect(strides(doubleperm))
end

using OffsetArrays
@testset "Static-Dynamic Size, Strides, and Offsets" begin
    A = zeros(3,4,5); Ap = @view(PermutedDimsArray(A,(3,1,2))[:,1:2,1])';
    S = @SArray zeros(2,3,4); Sp = @view(PermutedDimsArray(S,(3,1,2))[2:3,1:2,:]);
    M = @MArray zeros(2,3,4); Mp = @view(PermutedDimsArray(M,(3,1,2))[:,2,:])';
    Sp2 = @view(PermutedDimsArray(S,(3,2,1))[2:3,:,:]);
    Mp2 = @view(PermutedDimsArray(M,(3,1,2))[2:3,:,2])';

    @test @inferred(ArrayInterface.size(A)) === (3,4,5)
    @test @inferred(ArrayInterface.size(Ap)) === (2,5)
    @test @inferred(ArrayInterface.size(A)) === size(A)
    @test @inferred(ArrayInterface.size(Ap)) === size(Ap)

    @test @inferred(ArrayInterface.size(S)) === (StaticInt(2), StaticInt(3), StaticInt(4))
    @test @inferred(ArrayInterface.size(Sp)) === (2, 2, StaticInt(3))
    @test @inferred(ArrayInterface.size(Sp2)) === (2, StaticInt(3), StaticInt(2))
    @test @inferred(ArrayInterface.size(S)) == size(S)
    @test @inferred(ArrayInterface.size(Sp)) == size(Sp)
    @test @inferred(ArrayInterface.size(Sp2)) == size(Sp2)
    @test @inferred(ArrayInterface.size(Sp2, StaticInt(1))) === 2
    @test @inferred(ArrayInterface.size(Sp2, StaticInt(2))) === StaticInt(3)
    @test @inferred(ArrayInterface.size(Sp2, StaticInt(3))) === StaticInt(2)

    @test @inferred(ArrayInterface.size(M)) === (StaticInt(2), StaticInt(3), StaticInt(4))
    @test @inferred(ArrayInterface.size(Mp)) === (StaticInt(3), StaticInt(4))
    @test @inferred(ArrayInterface.size(Mp2)) === (StaticInt(2), 2)
    @test @inferred(ArrayInterface.size(M)) == size(M)
    @test @inferred(ArrayInterface.size(Mp)) == size(Mp)
    @test @inferred(ArrayInterface.size(Mp2)) == size(Mp2)

    @test @inferred(ArrayInterface.strides(A)) === (StaticInt(1), 3, 12)
    @test @inferred(ArrayInterface.strides(Ap)) === (StaticInt(1), 12)
    @test @inferred(ArrayInterface.strides(A)) == strides(A)
    @test @inferred(ArrayInterface.strides(Ap)) == strides(Ap)

    @test @inferred(ArrayInterface.strides(S)) === (StaticInt(1), StaticInt(2), StaticInt(6))
    @test @inferred(ArrayInterface.strides(Sp)) === (StaticInt(6), StaticInt(1), StaticInt(2))
    @test @inferred(ArrayInterface.strides(Sp2)) === (StaticInt(6), StaticInt(2), StaticInt(1))
    @test @inferred(ArrayInterface.stride(Sp2, StaticInt(1))) === StaticInt(6)
    @test @inferred(ArrayInterface.stride(Sp2, StaticInt(2))) === StaticInt(2)
    @test @inferred(ArrayInterface.stride(Sp2, StaticInt(3))) === StaticInt(1)

    @test @inferred(ArrayInterface.strides(M)) === (StaticInt(1), StaticInt(2), StaticInt(6))
    @test @inferred(ArrayInterface.strides(Mp)) === (StaticInt(2), StaticInt(6))
    @test @inferred(ArrayInterface.strides(Mp2)) === (StaticInt(1), StaticInt(6))
    @test @inferred(ArrayInterface.strides(M)) == strides(M)
    @test @inferred(ArrayInterface.strides(Mp)) == strides(Mp)
    @test @inferred(ArrayInterface.strides(Mp2)) == strides(Mp2)

    @test @inferred(ArrayInterface.offsets(A)) === (StaticInt(1), StaticInt(1), StaticInt(1))
    @test @inferred(ArrayInterface.offsets(Ap)) === (StaticInt(1), StaticInt(1))

    @test @inferred(ArrayInterface.offsets(S)) === (StaticInt(1), StaticInt(1), StaticInt(1))
    @test @inferred(ArrayInterface.offsets(Sp)) === (StaticInt(1), StaticInt(1), StaticInt(1))
    @test @inferred(ArrayInterface.offsets(Sp2)) === (StaticInt(1), StaticInt(1), StaticInt(1))

    @test @inferred(ArrayInterface.offsets(M)) === (StaticInt(1), StaticInt(1), StaticInt(1))
    @test @inferred(ArrayInterface.offsets(Mp)) === (StaticInt(1), StaticInt(1))
    @test @inferred(ArrayInterface.offsets(Mp2)) === (StaticInt(1), StaticInt(1))

    O = OffsetArray(A, 3, 7, 10);
    Op = PermutedDimsArray(O,(3,1,2));
    @test @inferred(ArrayInterface.offsets(O)) === (4, 8, 11)
    @test @inferred(ArrayInterface.offsets(Op)) === (11, 4, 8)

    @test @inferred(ArrayInterface.offsets((1,2,3))) === (StaticInt(1),)
end

@test ArrayInterface.can_avx(ArrayInterface.can_avx) == false

@testset "can_change_size" begin
    @test ArrayInterface.can_change_size([1])
    @test ArrayInterface.can_change_size(Vector{Int})
    @test ArrayInterface.can_change_size(Dict{Symbol,Any})
    @test !ArrayInterface.can_change_size(Base.ImmutableDict{Symbol,Int64})
    @test !ArrayInterface.can_change_size(Tuple{})
end

@testset "known_length" begin
    @test ArrayInterface.known_length(@inferred(ArrayInterface.indices(SOneTo(7)))) == 7
    @test ArrayInterface.known_length(1:2) == nothing
    @test ArrayInterface.known_length((1,)) == 1
    @test ArrayInterface.known_length((a=1,b=2)) == 2
    @test ArrayInterface.known_length([]) == nothing

    x = view(SArray{Tuple{3,3,3}}(ones(3,3,3)), :, SOneTo(2), 2)
    @test @inferred(ArrayInterface.known_length(x)) == 6
    @test @inferred(ArrayInterface.known_length(x')) == 6

    itr = StaticInt(1):StaticInt(10)
    @inferred(ArrayInterface.known_length((i for i in itr))) == 10
end

@testset "indices" begin
    A23 = ones(2,3); SA23 = @SMatrix ones(2,3);
    A32 = ones(3,2); SA32 = @SMatrix ones(3,2);
    @test @inferred(ArrayInterface.indices((A23, A32))) == 1:6
    @test @inferred(ArrayInterface.indices((SA23, A32))) == 1:6
    @test @inferred(ArrayInterface.indices((A23, SA32))) == 1:6
    @test @inferred(ArrayInterface.indices((SA23, SA32))) == 1:6
    @test @inferred(ArrayInterface.indices(A23)) == 1:6
    @test @inferred(ArrayInterface.indices(SA23)) == 1:6
    @test @inferred(ArrayInterface.indices(A23, 1)) == 1:2
    @test @inferred(ArrayInterface.indices(SA23, StaticInt(1))) === Base.Slice(StaticInt(1):StaticInt(2))
    @test @inferred(ArrayInterface.indices((A23, A32), (1, 2))) == 1:2
    @test @inferred(ArrayInterface.indices((SA23, A32), (StaticInt(1), 2))) === Base.Slice(StaticInt(1):StaticInt(2))
    @test @inferred(ArrayInterface.indices((A23, SA32), (1, StaticInt(2)))) === Base.Slice(StaticInt(1):StaticInt(2))
    @test @inferred(ArrayInterface.indices((SA23, SA32), (StaticInt(1), StaticInt(2)))) === Base.Slice(StaticInt(1):StaticInt(2))
    @test @inferred(ArrayInterface.indices((A23, A23), 1)) == 1:2
    @test @inferred(ArrayInterface.indices((SA23, SA23), StaticInt(1))) === Base.Slice(StaticInt(1):StaticInt(2))
    @test @inferred(ArrayInterface.indices((SA23, A23), StaticInt(1))) === Base.Slice(StaticInt(1):StaticInt(2))
    @test @inferred(ArrayInterface.indices((A23, SA23), StaticInt(1))) === Base.Slice(StaticInt(1):StaticInt(2))
    @test @inferred(ArrayInterface.indices((SA23, SA23), StaticInt(1))) === Base.Slice(StaticInt(1):StaticInt(2))
    @test_throws AssertionError ArrayInterface.indices((A23, ones(3, 3)), 1)
    @test_throws AssertionError ArrayInterface.indices((A23, ones(3, 3)), (1, 2))
    @test_throws AssertionError ArrayInterface.indices((SA23, ones(3, 3)), StaticInt(1))
    @test_throws AssertionError ArrayInterface.indices((SA23, ones(3, 3)), (StaticInt(1), 2))
    @test_throws AssertionError ArrayInterface.indices((SA23, SA23), (StaticInt(1), StaticInt(2)))

    @test size(similar(ones(2, 4), ArrayInterface.indices(ones(2, 4), 1), ArrayInterface.indices(ones(2, 4), 2))) == (2, 4)
    @test axes(ArrayInterface.indices(ones(2,2))) === (StaticInt(1):4,)
    @test axes(Base.Slice(StaticInt(2):4)) === (Base.IdentityUnitRange(StaticInt(2):4),)
    @test Base.axes1(ArrayInterface.indices(ones(2,2))) === StaticInt(1):4
    @test Base.axes1(Base.Slice(StaticInt(2):4)) === Base.IdentityUnitRange(StaticInt(2):4)
end

@testset "StaticInt" begin
    @test iszero(StaticInt(0))
    @test !iszero(StaticInt(1))
    @test @inferred(one(StaticInt(1))) === StaticInt(1)
    @test @inferred(zero(StaticInt(1))) === StaticInt(0)
    @test @inferred(one(StaticInt)) === StaticInt(1)
    @test @inferred(zero(StaticInt)) === StaticInt(0)
    @test eltype(one(StaticInt)) <: Int

    x = StaticInt(1)
    @test @inferred(Bool(x)) isa Bool
    @test @inferred(BigInt(x)) isa BigInt
    @test @inferred(Integer(x)) === x
    # test for ambiguities and correctness
    for i ∈ [StaticInt(0), StaticInt(1), StaticInt(2), 3]
        for j ∈ [StaticInt(0), StaticInt(1), StaticInt(2), 3]
            i === j === 3 && continue
            for f ∈ [+, -, *, ÷, %, <<, >>, >>>, &, |, ⊻, ==, ≤, ≥]
                (iszero(j) && ((f === ÷) || (f === %))) && continue # integer division error
                @test convert(Int, @inferred(f(i,j))) == f(convert(Int, i), convert(Int, j))
            end
        end
        i == 3 && break
        for f ∈ [+, -, *, /, ÷, %, ==, ≤, ≥]
            x = f(convert(Int, i), 1.4)
            y = f(1.4, convert(Int, i))
            @test convert(typeof(x), @inferred(f(i, 1.4))) === x
            @test convert(typeof(y), @inferred(f(1.4, i))) === y # if f is division and i === StaticInt(0), returns `NaN`; hence use of ==== in check.
        end
    end
end

@testset "insert/deleteat" begin
    @test @inferred(ArrayInterface.insert([1,2,3], 2, -2)) == [1, -2, 2, 3]
    @test @inferred(ArrayInterface.deleteat([1, 2, 3], 2)) == [1, 3]

    @test @inferred(ArrayInterface.deleteat([1, 2, 3], [1, 2])) == [3]
    @test @inferred(ArrayInterface.deleteat([1, 2, 3], [1, 3])) == [2]
    @test @inferred(ArrayInterface.deleteat([1, 2, 3], [2, 3])) == [1]

    @test @inferred(ArrayInterface.insert((2,3,4), 1, -2)) == (-2, 2, 3, 4)
    @test @inferred(ArrayInterface.insert((2,3,4), 2, -2)) == (2, -2, 3, 4)
    @test @inferred(ArrayInterface.insert((2,3,4), 3, -2)) == (2, 3, -2, 4)

    @test @inferred(ArrayInterface.deleteat((2, 3, 4), 1)) == (3, 4)
    @test @inferred(ArrayInterface.deleteat((2, 3, 4), 2)) == (2, 4)
    @test @inferred(ArrayInterface.deleteat((2, 3, 4), 3)) == (2, 3)
    @test ArrayInterface.deleteat((1, 2, 3), [1, 2]) == (3,)
end

include("indexing.jl")

ext: jl


text:
using ProgressMeter
using PmapProgressMeter

function initial_states(problems::Vector{OPCSP}; rng_offset::Int=1000)
    states = Array(OPCSPState, length(problems))
    for j in 1:length(states)
        p = problems[j]
        states[j] = rand!(MersenneTwister(j+rng_offset), create_state(p), initial_belief(p))
    end
    return states
end

function test_run(p::OPCSP, is::OPCSPState, solver::Union{MCTS.DPWSolver,MCTS.AgUCTSolver}, rng::AbstractRNG=MersenneTwister())
    policy = MCTSAdapter(POMDPs.solve(solver, OPCSPBeliefMDP(p)))
    sim = POMDPToolbox.HistoryRecorder(rng=rng, initial_state=is)
    simr = simulate(sim, p, policy, OPCSPUpdater(p), initial_belief(p))
    path = Int[s.i for s in sim.state_hist]
    r = reward(p, is.d, path)
    @assert abs(simr + p.r[p.stop] + is.d[p.stop] - r) < 1e-6
    return r
end

function test_run(p::OPCSP, is::OPCSPState, policy::Policy, rng::AbstractRNG=MersenneTwister())
    sim = POMDPToolbox.HistoryRecorder(rng=rng, initial_state=is)
    simr = simulate(sim, p, policy, updater(policy), initial_belief(p))
    path = Int[s.i for s in sim.state_hist]
    r = reward(p, is.d, path)
    @assert abs(simr + p.r[p.stop] + is.d[p.stop] - r) < 1e-6
    return r
end

function test_run(p::OPCSP, is::OPCSPState, solver::Solver, rng::AbstractRNG=MersenneTwister())
    policy = POMDPs.solve(solver, p)
    return test_run(p, is, policy, rng)
end

function test_run(p::OPCSP, is::OPCSPState, solver::OPSolver, rng::AbstractRNG=MersenneTwister())
    path = solve_op(solver, p)
    return reward(p, is.d, path)
end

function test_run(p::OPCSP, is::OPCSPState, solver::Cheater, rng::AbstractRNG=MersenneTwister())
    path = cheat(p, is.d)
    return reward(p, is.d, path)
end

function test_run(p::OPCSP, is::OPCSPState, solver::OnlyUnobservableUncertainty, rng::AbstractRNG=MersenneTwister())
    path = cheat_observe_others(p, is.d)
    return reward(p, is.d, path)
end

function evaluate_performance(problems::Vector{OPCSP}, iss::Vector{OPCSPState}, solver; rng_offset::Int=100, parallel=true)
    # rewards = SharedArray(Float64, length(problems))
    N = length(problems)
    if parallel
        prog = Progress( N, dt=0.1, barlen=50, output=STDERR)
        rewards = pmap(test_run,
                       prog,
                       problems,
                       iss,
                       [solver for i in 1:N],
                       [MersenneTwister(j+rng_offset) for j in 1:N])
    else
        rewards = Array(Float64, length(problems))
        @showprogress for j in 1:length(problems)
            rewards[j] = test_run(problems[j], iss[j], solver, MersenneTwister(j+rng_offset))
        end
    end
    return sdata(rewards)
end


text:
module OrbitalIndices
export OrbitalIndex, @nl_str

using ArgCheck

""" Indexing with respect to first two quantum numbers """
immutable OrbitalIndex{I <: Integer}
    """ Primary quantum number """
    n::I
    """ Orbital momentum number """
    l::I

    OrbitalIndex{I}(n::I, l::I) where I <: Integer = begin
        @argcheck n ≥ 0
        @argcheck 0 ≤ l
        @argcheck l ≤ n
        new(n, l)
    end
end

OrbitalIndex(n::Integer, l::Integer) = begin
    @argcheck n ≥ 0
    @argcheck 0 ≤ l
    @argcheck l ≤ n
    OrbitalIndex{promote_type(typeof(n), typeof(l))}(n, l)
end

Base.isequal(a::OrbitalIndex, b::OrbitalIndex) = a.n ==  b.n && a.l == b.l
Base.isless(a::OrbitalIndex, b::OrbitalIndex) = a.n < b.n || (a.n == b.n && a.l < b.l)

Base.convert(::Type{Integer}, i::OrbitalIndex) =
    convert(typeof(i).parameters[1], (i.n * i.n + i.n) / 2) + i.l + 1
Base.convert(::Type{OrbitalIndex}, i::Integer) = begin
    n = trunc(typeof(i), (-1 + sqrt(8i - 7))/2)
    @assert n * (n + 1) ≤ 2i ≤ (n + 1) * (n + 2)
    OrbitalIndex{typeof(i)}(n, i - 1 - convert(typeof(i), (n * n + n) / 2))
end

Base.print(io::IO, orb::OrbitalIndex) = print(io, "n=", orb.n, ", l=", orb.l)

macro nl_str(s)
    args = [parse(Int64, u) for u in split(s, ",")]
    @argcheck length(args) == 2
    quote
        OrbitalIndex($(args[1]), $(args[2]))
    end
end

""" Unit range of orbitals with quantum numbers n, l """
immutable OrbitalIndexUnitRange{I <: Integer} <: AbstractUnitRange{OrbitalIndex{I}}
    start::OrbitalIndex{I}
    stop::OrbitalIndex{I}
end

""" Step range of through orbitals with quantum numbers n, l

    At each iteration, the step (n', l') is added to the current index. If l' + l₀ > n' +
    n₀, then, the iteration wraps to n₁ = n' + n₀ and l₁ = l' + l₀ - n' - n₀.
"""
immutable OrbitalIndexRange{I <: Integer} <: OrdinalRange{OrbitalIndex{I}, OrbitalIndex{I}}
  start::OrbitalIndex{I}
  step::OrbitalIndex{I}
  stop::OrbitalIndex{I}
end

function Base.colon(start::OrbitalIndex, stop::OrbitalIndex)
	OrbitalIndexUnitRange(promote(start, stop)...)
end

function Base.colon(start::OrbitalIndex, step::OrbitalIndex, stop::OrbitalIndex)
    OrbitalIndexRange(promote(start, step, stop)...)
end

Base.length(r::OrbitalIndexUnitRange) = convert(Integer, r.stop) - convert(Integer, r.start)
function Base.length(r::OrbitalIndexRange)
    index = 0
    for u in r
        index += 1
    end
    index
end

Base.start(r::Union{OrbitalIndexUnitRange, OrbitalIndexRange}) = r.start
function Base.next(::OrbitalIndexUnitRange, state::OrbitalIndex)
    if state.l == state.n
        nstate = OrbitalIndex(state.n + 1, 0)
    else
        nstate = OrbitalIndex(state.n, state.l + 1)
    end
    state, nstate
end
function Base.next(iter::OrbitalIndexRange, state::OrbitalIndex)
    state, OrbitalIndex(state.n + iter.step.n, state.l + iter.step.l)
end
Base.done(iter::OrbitalIndexUnitRange, state::OrbitalIndex) = state > iter.stop
Base.done(iter::OrbitalIndexRange, state::OrbitalIndex) = state > iter.stop
end


text:
using FastRationals
using BenchmarkTools

function hilbert_matrix(::Type{T}, n::Int=7) where T<:Signed
    num = one(FastRational{T})
    return [num // (T(idx) + T(inx) - one1) for idx in 1:n, inx in 1:n]
end


text:
value(x::SIUnits.SIQuantity) = x.val # While not in master: https://github.com/Keno/SIUnits.jl/issues/14

Height = quantity(Float64, SIUnits.Meter)
Volume = quantity(Float64, SIUnits.Meter^3)
Discharge = quantity(Float64, SIUnits.Meter^3/SIUnits.Second)
Power = quantity(Float64, SIUnits.Watt)

_from_unitful(x::Height) = value(x)
_from_unitful(x::Volume) = value(x) / 10^6
_from_unitful(x::Discharge) = value(x) / 10^6
_from_unitful(x::Power) = value(x) / 10^6
_from_unitful(x::Tuple{Volume, Volume}) = (_from_unitful(x[1]), _from_unitful(x[2]))

_from_unitful(x::Float64) = x
_from_unitful(x::Tuple{Float64, Float64}) = x

"""
This set of functions takes as input a yearly time series and output the beginning or the end of the hydrological season.
"""
immutable SeasonsSeparator
  date_begin::Dict{Symbol, Function}
  date_end::Dict{Symbol, Function}

  function SeasonsSeparator(date_begin::Dict{Symbol, Function}, date_end::Dict{Symbol, Function})
    # Check if IDs match. The order does not necessarily match, though! Hence the complex code with looping.
    beg_ids = collect(keys(date_begin))
    end_ids = collect(keys(date_end))

    for id in beg_ids
      if ! in(id, end_ids)
        error("The seasons do not match: found a beginning for season " * string(id) * ", but no end.")
      end
    end
    for id in end_ids
      if ! in(id, beg_ids)
        error("The seasons do not match: found an end for season " * string(id) * ", but no beginning.")
      end
    end

    return new(date_begin, date_end)
  end
end

"""
A seasons separator for Western Europe, based on hydrological years (starts in October, ends in September).
           October   April	|	May			September
             1st     30th   | 1st       30th
                  Wet       |       Dry
"""
HydrologicalWesternEurope = SeasonsSeparator(
  Dict(:Wet => y -> DateTime(y, October, 1),
       :Dry => y -> DateTime(y + 1, May, 1)),
  Dict(:Wet => y -> DateTime(y + 1, April, Dates.daysinmonth(y + 1, April)),
       :Dry => y -> DateTime(y + 1, September, Dates.daysinmonth(y + 1, September)))
)

seasons(seasonsSep::SeasonsSeparator) = collect(keys(seasonsSep.date_begin))
isSeasonKnown(seasonsSep::SeasonsSeparator, season::Symbol) = in(season, seasons(seasonsSep))
function ensureSeasonKnown(seasonsSep::SeasonsSeparator, season::Symbol)
  if ! isSeasonKnown(seasonsSep, season)
    error("Unknown season " * string(seasons) * ". " *
          "The defined seasons for this separator are: " * string(seasons(seasonsSep)))
  end
end
function begins(seasonsSep::SeasonsSeparator, season::Symbol, chosenYear::Int)
  ensureSeasonKnown(seasonsSep, season)
  return seasonsSep.date_begin[season](chosenYear)
end
function ends(seasonsSep::SeasonsSeparator, season::Symbol, chosenYear::Int)
  ensureSeasonKnown(seasonsSep, season)
  return seasonsSep.date_end[season](chosenYear)
end
length(seasonsSep::SeasonsSeparator, season::Symbol, chosenYear::Int) =
  Day(length(begins(seasonsSep, season, chosenYear) : ends(seasonsSep, season, chosenYear)))
function length(seasonsSep::SeasonsSeparator, season::Symbol, chosenYear::Int, period::Period)
  daily = Dates.days(length(seasonsSep, season, chosenYear))
  days_per_period = Dates.days(period)
  return round(Int, daily / days_per_period)
end

"""
Computes the length of each season while guaranteeing that, for the given `period`, all seasons have the same number of
time steps, whatever the value of `chosenYear` is.
"""
function seasonBoundariesGuarantee(seasonsSep::SeasonsSeparator, period::Period, chosenYear::Int)
  ## Idea: start from a basic decomposition based on the actual year; then reassign time steps at the borders so
  ## that each season is closer to its preferred length.

  # Compute the preferred lengths for an arbitrary year (here, 2000, which is not leap).
  preferred = [season => length(seasonsSep, season, 2000, period) for season in seasons(seasonsSep)]

  # Base beginnings and ends, with season lengths.
  bs = [season => begins(seasonsSep, season, chosenYear) for season in seasons(seasonsSep)] # Begin
  es = [season =>   ends(seasonsSep, season, chosenYear) for season in seasons(seasonsSep)] # End
  ls = [season => length(seasonsSep, season, chosenYear, period) for season in seasons(seasonsSep)] # Length

  # Quantify the delta between those values and the preferred ones: if some periods are missing, decide to take them
  # before the beginning of the first season. Otherwise, the delta is not important for now: it will be compensated
  # at the end of the seasons with the next loop.
  total_delta = 0
  for season in seasons(seasonsSep)
    delta = ls[season] - preferred[season]
    if delta < 0 # Length too small for this season: should have started before.
      total_delta += abs(delta)
    elseif delta > 0 # Too long: should have started after, i.e. compensate some shorter seasons.
      total_delta -= abs(delta)
    end
  end
  if total_delta > 0
    season = minimum(bs)[1]
    bs[season] -= total_delta * period
  end

  # Now, adapt the other boundaries in between: start from the first season, then add up the periods up
  # with the reference length.
  ordered_seasons = [pair[2] for pair in sort(collect(zip(values(bs), keys(bs))))]
  es[ordered_seasons[1]] = bs[ordered_seasons[1]] + (preferred[ordered_seasons[1]] - 1) * period
  for i in 2:length(ordered_seasons)
    prev_season = ordered_seasons[i - 1]
    season = ordered_seasons[i]
    bs[season] = es[prev_season] + period
    es[season] = bs[season] + (preferred[season] - 1) * period
  end

  # Done!
  return bs, es
end

"""
Separates seasons for the given time series `ts`. The seasons begin and end as imposed by `seasonsSep`.
The separation is performed for year `chosenYear`. The output is a dictionary, where the time series are indexed
by the season identifier.

When `guaranteeTimeStep` is `true`, then the seasons are such that all years have the same number of
time steps (whose duration is the same as that of `ts`). This option might introduce slight shifts in the season
(i.e. it may begin or end at dates slightly different from `seasonsSep`).
"""
function separate(ts::TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}, seasonsSep::SeasonsSeparator,
                  chosenYear::Int=year(minimum(timestamp(ts))); guaranteeTimeStep::Bool=false)
  # Simple case: just do it.
  if ! guaranteeTimeStep
    rets = Dict{Symbol, TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}}()
    for season in seasons(seasonsSep)
      rets[season] = to(from(ts, begins(seasonsSep, season, chosenYear)), ends(seasonsSep, season, chosenYear))
    end
    return rets
  end

  # Make the seasons dictionary with the time step guarantee.
  (bs, es) = seasonBoundariesGuarantee(seasonsSep, meta(ts), chosenYear)
  rets = Dict{Symbol, TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}}()
  for season in seasons(seasonsSep)
    rets[season] = to(from(ts, bs[season]), es[season])
  end
  return rets
end

function separate(ts::TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}, seasonsSep::SeasonsSeparator,
                  season::Symbol, chosenYear::Int=year(minimum(timestamp(ts))); guaranteeTimeStep::Bool=false)
  # Simple case: just do it.
  if ! guaranteeTimeStep
    rets = Dict{Symbol, TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}}()
    return to(from(ts, begins(seasonsSep, season, chosenYear)), ends(seasonsSep, season, chosenYear))
  end

  # Make the seasons dictionary with the time step guarantee.
  (bs, es) = seasonBoundariesGuarantee(seasonsSep, meta(ts), chosenYear)
  return to(from(ts, bs[season]), es[season])
end

"""
Separates seasons for the given array of time series.
"""
separate(TS::Array{TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}, 1}, seasons::SeasonsSeparator, season::Symbol, year::Int) =
  [separate(ts, seasons, season, year) for ts in TS]
separate(TS::Array{TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}, 1}, seasons::SeasonsSeparator, season::Symbol) =
  [separate(ts, seasons, season) for ts in TS]
separate(TS::Array{TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}, 1}, seasons::SeasonsSeparator, year::Int) =
  [separate(ts, seasons, year) for ts in TS]
separate(TS::Array{TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}, 1}, seasons::SeasonsSeparator) =
  [separate(ts, seasons) for ts in TS]

"""
Collects seasons for the given time series: throughout the time series, it collects the asked season (e.g., all time steps within the wet season);
for one-year series, it is equivalent to `separate`.
"""
function Base.collect(ts::TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}, seasons::SeasonsSeparator, season::Symbol)
  minY = year(minimum(timestamp(ts)))
  maxY = year(maximum(timestamp(ts)))
  return vcat([separate(ts, seasons, season, year) for year in minY:maxY]...)
end

"""
A non-leap year iterator. State: current year in iteration, number of already generated years.

It is more efficient than filtering over years with the built-in function `isleapyear` (Julia 0.4.6):
  @time collect(NonLeapYearIterator(1, 1_000))
       0.000059 seconds (648 allocations: 26.734 KB)
  @time take(filter(x -> ! isleapyear(x), 1:1_500), 1_000)
       0.001695 seconds (1.03 k allocations: 28.041 KB)
"""
immutable NonLeapYearIterator
  start::Int
  n::Int
end

typealias NonLeapYearIteratorState Tuple{Int64, Int64}
Base.start(nly::NonLeapYearIterator) = (nly.start, 0)
function Base.next(nly::NonLeapYearIterator, state::NonLeapYearIteratorState)
  y = state[1] + 1
  if isleapyear(y) # Never two consecutive leap years.
    y = y + 1
  end
  return y, (y, state[2] + 1)
end
Base.done(nly::NonLeapYearIterator, state::NonLeapYearIteratorState) = state[2] >= nly.n

"""
These types allows considering an array of time series as scenarios within the `feasibilityProbability` solver.
"""
immutable ArrayAsIterableReservoirScenarios
  scenarios::Array{TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}, 1}
end
Base.start(ars::ArrayAsIterableReservoirScenarios) = 1
Base.next(ars::ArrayAsIterableReservoirScenarios, state::Int) =
  ArrayAsIterableReservoirScenario(values(ars.scenarios[state])), state + 1
Base.done(ars::ArrayAsIterableReservoirScenarios, state::Int) = state >= length(ars.scenarios)
eachscenario(lts::Array{TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}, 1}) =
  ArrayAsIterableReservoirScenarios(lts)

immutable ArrayAsIterableReservoirScenario # Define a new type to avoid overloading totalInflow() with too generic types.
  scenarios::Array{Float64, 1}
end
totalInflow(ars::ArrayAsIterableReservoirScenario) = ars.scenarios

period2seconds(p::Period) = Dates.days(p) * 24 * 60 * 60

"""
Ensures the given `array` has exactly the length `givenLength`, either by removing parts of it, or by copying
the whole array an integer number of times (otherwise, an error is thrown).
"""
function toLength(array::Array{Float64, 1}, givenLength::Int)
  if length(array) == givenLength
    return array
  elseif length(array) > givenLength
    return array[1:givenLength]
  elseif length(array) < givenLength
    n_copies = try
      Int(givenLength // length(array))
    catch
      exact_str = string(givenLength) * "//" * string(length(array)) # Direct string(//) will reduce the fraction.
      error("To match the expected length, would need to copy a noninteger number of times " *
            "(" * exact_str            * " ~= " * string(givenLength / length(array)) * ").")
      #           Exact version (like 1//2)       Approximation (like 0.5)
    end
    return vec(repmat(array, 1, n_copies))
  end
end


text:
sorted_unique_count(m) = isempty(m) ? 0 : 1+sum(m[2:end].≠m[1:end-1])

inverse_power_intensity_function(α::Real) = m -> 1/sorted_unique_count(m)^α
all_or_nothing_intensity_function(α::Real, β::Real) = m -> all(i -> i == first(m), m) ? α : β

max_stars_repo_name: LilithHafner/FBD

max_issues_repo_name: LilithHafner/FBD

max_forks_repo_name: LilithHafner/FBD


text:
using ArgParse, Match
s = ArgParseSettings()
@add_arg_table! s begin
    "--configuration";
    "-T";                      default=Float32;          eval_arg=true
    "--Nside";                 default=128;              arg_type=Int
    "--θpix";                  default=3;                arg_type=Int
    "--pol";                   default=:P;               arg_type=Symbol
    "--beamFWHM";              default=1.;               arg_type=Float64
    "--μKarcminT";             default=1.;               arg_type=Float64
    "--bandpass_mask";         default="LowPass(3000)"
    "--pixel_mask_kwargs";     eval_arg=true
    "--storage";               default="Array";          range_tester=in(("Array","CuArray"))
    "--ϕstart";                default="quasi_sample";   range_tester=in(("quasi_sample","bestfit","0"))
    "--Nϕ_fac";                default=2f0;              arg_type=Float32
    "--filename"
    "--resume";                action=:store_true
    "--symp_kwargs";           default=[(N=25, ϵ=0.01)]; eval_arg=true
    "--nchains";               default=1;                arg_type=Int
    "--nsamps_per_chain";      default=2000;             arg_type=Int
    "--nlenseflow_ode";        default=10;               arg_type=Int
    "--seed";                  arg_type=Int;
    "--nchunk";                default=30;               arg_type=Int
    "--nsavemaps";             default=5;                arg_type=Int
    "--nburnin_always_accept"; default=50;               arg_type=Int
    "--nburnin_fixθ";          default=10;               arg_type=Int
    "--sampled_params";        default=[:Aϕ,:r];         eval_arg=true
    "--rfid";                  default=0.05;             arg_type=Float64
end
args = parse_args(ARGS, s)

###

# configurations from the paper
merge!(args, @match args["configuration"] begin
    "2PARAM" => begin
        Dict(
            "pol"               => :P,
            "Nside"             => 256,
            "θpix"              => 2,
            "μKarcminT"         => 1/√2,
            "beamFWHM"          => 2,
            "bandpass_mask"     => "LowPass(5000)",
            "pixel_mask_kwargs" => (edge_padding_deg=0.4, apodization_deg=0.6, edge_rounding_deg=0.1, num_ptsrcs=0),
            "sampled_params"    => [:Aϕ,:r],
            "rfid"              => 0.04,
            "symp_kwargs"       => [(N=25, ϵ=0.02)]
        )
    end
    "MANY" => begin
        if !(args["rfid"] in [0,0.02,0.04])
            @warn "rfid=$(args["rfid"]) inconsistent with configuration $(args["configuration"])"
        end
        Dict(
            "pol"               => :P,
            "Nside"             => 256,
            "θpix"              => 3,
            "μKarcminT"         => 1/√2,
            "beamFWHM"          => 3,
            "bandpass_mask"     => "LowPass(3500)",
            "pixel_mask_kwargs" => (edge_padding_deg=0.6, apodization_deg=0.9, edge_rounding_deg=0.1, num_ptsrcs=0),
            "sampled_params"    => [:r],
            "symp_kwargs"       => [(N=25, ϵ=0.02)]
        )
    end
    "BIG" => begin
        if !(args["rfid"] in [0,0.01,0.02])
            @warn "rfid=$(args["rfid"]) inconsistent with configuration $(args["configuration"])"
        end
        Dict(
            "pol"               => :P,
            "Nside"             => 512,
            "θpix"              => 3,
            "μKarcminT"         => 1/√2,
            "beamFWHM"          => 3,
            "bandpass_mask"     => "LowPass(3500)",
            "pixel_mask_kwargs" => (edge_padding_deg=1.2, apodization_deg=1.8, edge_rounding_deg=0.2, num_ptsrcs=0),
            "sampled_params"    => [:r],
            "symp_kwargs"       => [(N=25, ϵ=0.02)]
        )
    end
    _::Nothing => Dict()
    c => error("Unrecognized configuration: $c")
end)

if args["storage"]=="CuArray"
    using CuArrays
end
args["storage"] = eval(Symbol(args["storage"]))
args["ϕstart"] = args["ϕstart"] == "0" ? 0 : Symbol(args["ϕstart"]);

## if running multiple chains & GPUs, set each process to use one GPU

using Distributed

if args["nchains"]>1 && args["storage"]==CuArray

    using CUDAdrv
    length(devices()) == args["nchains"] || error("nchains=$(args["nchains"]) but $(length(devices())) GPUs found.")
    addprocs(length(devices()))
    @everywhere workers() begin
        using Pkg
        # activate the same environment on the workers as the master
        Pkg.activate($(Pkg.API.Context().env.project_file))
        using CUDAnative, CuArrays
        # until the fixes to https://github.com/JuliaGPU/CuArrays.jl/issues/589 hit a release:
        CuArrays.CURAND.seed!(rand(0:typemax(Int)))
    end

    # assign devices
    asyncmap((zip(workers(), devices()))) do (p, d)
        remotecall_wait(p) do
            @info "Worker $p uses $d"
            device!(d)
        end
    end

end

### load dataset

@info "Loading code..."
@time @everywhere begin
    using CMBLensing
    using CMBLensing: grid_and_sample
end

@info "Loading dataset..."
@unpack f, f̃, ϕ, ds = @time load_sim_dataset(
    Nside             = args["Nside"],
    T                 = args["T"],
    θpix              = args["θpix"],
    pol               = args["pol"],
    L                 = LenseFlow{RK4Solver{args["nlenseflow_ode"]}},
    seed              = args["seed"],
    rfid              = args["rfid"],
    beamFWHM          = args["beamFWHM"],
    μKarcminT         = args["μKarcminT"],
    pixel_mask_kwargs = args["pixel_mask_kwargs"],
    bandpass_mask     = eval(Meta.parse(args["bandpass_mask"])),
    Nϕ_fac            = args["Nϕ_fac"],
);

### baylens2-specific Gibbs θ-pass

@everywhere function baylens2_gibbs_pass_θ(;kwargs...)

    @unpack f°,ϕ°,θ,ds,θrange,progress = kwargs

    lnP_Aϕ, θ_Aϕ, lnP_r, θ_r = (), (), (), ()

    # sample Aϕ
    if :Aϕ in keys(θrange)
        lnP_Aϕ, θ_Aϕ = grid_and_sample(
            θ′->lnP(:mix,f°,ϕ°,merge(θ,θ′),ds),
            (Aϕ=θrange.Aϕ,),
            progress=(progress==:verbose)
        )
        lnP_Aϕ = (Aϕ=lnP_Aϕ,)
        θ = merge(θ,θ_Aϕ)
    end

    # sample r
    if :r in keys(θrange)
        lnP_r, θ_r = grid_and_sample(
            θ′->lnP(:mix,f°,ϕ°,merge(θ,θ′),ds),
            (r=θrange.r,),
            progress=(progress==:verbose)
        )
        lnP_r = (r=lnP_r,)
        θ = merge(θ,θ_r)
    end

    (;lnP_Aϕ..., lnP_r...), θ

end

## run sampler

θrange₀ = (
    Aϕ = range(0.75, 1.25, length=50),
    r  = range(sqrt(1e-6), sqrt(1e-1), length=50).^2
)

if !isinteractive()

    @info "Starting chain..."
    @unpack (rundat, chains) = sample_joint(

        ds,

        nchains = args["nchains"],

        nsamps_per_chain      = args["nsamps_per_chain"],
        storage               = args["storage"],
        chains                = args["resume"] ? :resume : nothing,
        nchunk                = args["nchunk"],
        nsavemaps             = args["nsavemaps"],
        nburnin_always_accept = args["nburnin_always_accept"],
        nburnin_fixθ          = args["nburnin_fixθ"],

        θrange = (;(k=>v for (k,v) in pairs(θrange₀) if k in args["sampled_params"])...),
        gibbs_pass_θ = baylens2_gibbs_pass_θ,

        ϕstart = args["ϕstart"],
        Nϕ_fac = args["Nϕ_fac"],

        wf_kwargs   = (tol=1e-1, nsteps=200),
        symp_kwargs = args["symp_kwargs"],
        MAP_kwargs  = (αmax=0.3, nsteps=5),

        progress = :summary,
        interruptable = true,

        filename = args["filename"]

    );

end


text:
const graphics = (W"Graphics", W"Graphics3D", W"Legended")

Base.Multimedia.showable(::MIME"image/svg+xml", w::MathLink.WExpr) = in(w.head, graphics)
Base.show(io::IO, ::MIME"image/svg+xml", w::MathLink.WExpr) = print(io, weval(W"ExportString"(w, "SVG")))

Base.Multimedia.showable(::MIME"image/png", w::MathLink.WExpr) = in(w.head, graphics)
function Base.show(io::IO, ::MIME"image/png", w::MathLink.WExpr)
    x = weval(W"ExportByteArray"(w, "PNG"))
    write(io, x.args[1].args)
end


text:
using Documenter

makedocs(sitename="RANSACVisualizer.jl",
        pages=["Home" => "index.md"])

deploydocs(repo="github.com/cserteGT3/RANSACVisualizer.jl.git")


text:
using SchumakerSpline
using Test

# Run tests

println("Test with Float64s")
@time @test include("Test_with_Floats.jl")
println("Test with Ints")
@time @test include("Test_with_Ints.jl")
println("Test with Dates")
@time @test include("Test_with_Dates.jl")
println("Test Intercepts")
@time @test include("test_intercepts.jl")
println("Test Splicing of Splines")
@time @test include("test_splice_splines.jl")
println("Test Plots")
@time @test include("test_plots.jl")
println("Test Extrapolation")
@time @test include("test_extrapolation.jl")


text:
import FileIO#, LinearAlgebra, JSON
dir = @__DIR__

erlangDParams = FileIO.load(dir*"/erlangParamsData/erlangDParams.jld2", "erlangDParams")

"""
    cme_params

A Dict storing parameters of concentrated matrix exponential distributions. See https://github.com/ghorvath78/iltcme.
"""
cme_params = FileIO.load(dir*"/CMEParamsData/CMEParams.jld2", "CMEParams")

"""
    AbstractMatrixExponential
"""
abstract type AbstractMatrixExponential end

"""
MatrixExponential constructor method

    MatrixExponential(
        a::Array{Float64,2},
        S::Array{Float64,2},
        s::Array{Float64,1},
        D::Array{Float64}=[0],
    )

Inputs:
 - `a` a 1 by p Array of reals
 - `S` a p by p Array of reals
 - `s` a p by 1 Array of reals
 - `D` an optional argument, if empty then me.D is the identity,
    else is a p by p matrix
 Throws an error if the dimensions are inconsistent.
"""
struct MatrixExponential <:AbstractMatrixExponential
    a::Array{Float64,2}
    S::Array{Float64,2}
    s::Array{Float64,1}
    D::Array{Float64,2}
    function MatrixExponential(
        a::Array{Float64,2},
        S::Array{Float64,2},
        s::Array{Float64,1},
        D::Array{Float64,2},
    )

        s1 = size(a,1)
        s2 = size(a,2)
        s3 = size(S,1)
        s4 = size(S,2)
        s5 = size(s,1)
        s6 = size(s,2)
        s7 = size(D,1)
        s8 = size(D,2)
        checksquare(D)
        checksquare(S)
        test = (s1!=1) || (s6!=1) || any(([s2;s3;s4;s7;s8].-s5).!=0)
        test && throw(DimensionMismatch("Dimensions of ME representation not consistent"))
        return new(a,S,s,D)
    end
end
MatrixExponential(a::Array{Float64,2},S::Array{Float64,2},s::Array{Float64,1}) =
    MatrixExponential(a,S,s,Matrix{Float64}(LinearAlgebra.I(size(S,1))))

struct CMEOperator <: AbstractArray{Float64,2}
    S::Array{Float64,2}
end
Base.size(S::CMEOperator) = size(S.S)
Base.size(S::CMEOperator,i) = size(S.S,i)
Base.getindex(S::CMEOperator,i) = S.S[i]
Base.getindex(S::CMEOperator,i,j) = S.S[i,j]
Base.setindex!(S::CMEOperator,x,i) = (S.S[i]=x)
Base.setindex!(S::CMEOperator,x,i,j) = (S.S[i,j]=x)
Base.:*(S::CMEOperator,x::Real) = CMEOperator(S.S*x)
Base.:*(x::Real,S::CMEOperator) = CMEOperator(x*S.S)
function Base.:-(D::LinearAlgebra.Diagonal{Bool, Vector{Bool}},S::CMEOperator)
    !(size(D)==size(S))&&throw(DimensionMismatch("matrices must be same size"))
    S[1,1] = 1.0-S[1,1]
    two_n_blocks = (size(S,1)-1)
    for n in 2:2:two_n_blocks
        S[n,n]      = 1.0-S[n,n]
        S[n+1,n]    = -S[n+1,n]
        S[n+1,n+1]  = 1.0-S[n+1,n+1]
        S[n,n+1]    = -S[n,n+1]
    end
    return S
end

function LinearAlgebra.exp(Q::CMEOperator)
    S = Q.S
    expS = zeros(size(S))
    tmp = exp(S[1,1])
    expS[1,1] = tmp
    two_n_blocks = (size(S,1)-1)
    for n in 2:2:two_n_blocks
        kωt = S[n+1,n]
        expS[n,n]       =  tmp*cos(kωt)    # 1,1 top-left
        expS[n,n+1]     = -tmp*sin(kωt)    # 1,2 top-right
        expS[n+1,n]     =  tmp*sin(kωt)    # 2,1 bottom-left
        expS[n+1,n+1]   =  tmp*cos(kωt)    # 2,2 bottom-right
    end
    return CMEOperator(expS)
end

function LinearAlgebra.inv(Q::CMEOperator)
    S = Q.S
    invS = zeros(size(S))
    invS[1,1] = 1.0./S[1,1]

    two_n_blocks = (size(S,1)-1)
    for n in 2:2:two_n_blocks
        det_block = S[n,n]*S[n+1,n+1] - S[n+1,n]*S[n,n+1]
        invS[n,n]       =  S[n+1,n+1]/det_block # 1,1 top-left
        invS[n,n+1]     = -S[n,n+1]/det_block   # 1,2 top-right
        invS[n+1,n]     = -S[n+1,n]/det_block   # 2,1 bottom-left
        invS[n+1,n+1]   =  S[n,n]/det_block     # 2,2 bottom-right
    end
    return invS
end

_order(ME::AbstractMatrixExponential) = length(ME.a)

"""
Similar to MatrixExponential but has a neater implementation of `orbit()` due to nice structure.
"""
struct ConcentratedMatrixExponential <: AbstractMatrixExponential
    a::Array{Float64,2}
    S::CMEOperator
    s::Array{Float64,1}
    D::Array{Float64,2}
    function ConcentratedMatrixExponential(
        a::Array{Float64,2},
        S::CMEOperator,
        s::Array{Float64,1},
        D::Array{Float64,2},
    )
        s1 = size(a,1)
        s2 = size(a,2)
        s3 = size(S,1)
        s4 = size(S,2)
        s5 = size(s,1)
        s6 = size(s,2)
        s7 = size(D,1)
        s8 = size(D,2)
        checksquare(D)
        checksquare(S)
        test = (s1!=1) || (s6!=1) || any(([s2;s3;s4;s7;s8].-s5).!=0)
        if test
            error("Dimensions of ME representation not consistent")
        else
            return new(a,S,s,D)
        end
    end
end

function ConcentratedMatrixExponential(a::Array{Float64,2},S::CMEOperator,
    s::Array{Float64,1})
    return ConcentratedMatrixExponential(a,S,s,Matrix{Float64}(LinearAlgebra.I(size(S,1))))
end

pdf(me::AbstractMatrixExponential) = x->only(me.a*exp(me.S*x)*me.s)
"""
    pdf([a::Array{Float64, 2},] me::AbstractMatrixExponential)

Return a the probability density function of a matrix exponential as a function.

If `a` is not specified then the initial vector of the `me` is `me.a`

# Arguments:
- `a`: a row vector with length `size(me.S,1)`
"""
pdf(a::Array{Float64,2}, me::AbstractMatrixExponential) =
    (length(a)==size(me.S,1)) ? (x->only(a*exp(me.S*x)*me.s)) : throw(
        DomainError("a and me.S must have compatible size"))

pdf(me::AbstractMatrixExponential, x::Real) = pdf(me)(x)
"""
    pdf(a::Array{Float64, 2}, me::AbstractMatrixExponential, x::Real)

Return a the probability density function of a matrix exponential evaluated at `x`.
"""
pdf(a::Array{Float64,2}, me::AbstractMatrixExponential, x::Real) = pdf(a,me)(x)

pdf(me::AbstractMatrixExponential, x::Array{Float64}) = pdf(me).(x)
pdf(a::Array{Float64,2}, me::AbstractMatrixExponential, x::Array{Float64}) = pdf(a,me).(x)

ccdf(me::AbstractMatrixExponential) = x->sum(me.a*exp(me.S*x))
"""
    ccdf(a::Array{Float64, 2}, me::AbstractMatrixExponential)

Return a the complimentary cumulative distribution function (1-cdf) of a matrix exponential as a function.

If `a` is not specified then the initial vector of the `me` is `me.a`

# Arguments:
- `a`: a row vector with length `size(me.S,1)`
"""
ccdf(a::Array{Float64,2}, me::AbstractMatrixExponential) = (length(a)==size(me.S,1)) ? (x->sum(a*exp(me.S*x))) : throw(
    DomainError("a and me.S must have compatible size"))

ccdf(me::AbstractMatrixExponential, x::Real) = ccdf(me)(x)
"""
    ccdf(a::Array{Float64, 2}, me::AbstractMatrixExponential, x::Real)

Return a the complimentary cumulative distribution function (1-cdf) of a matrix exponential evaluated at `x`.
"""
ccdf(a::Array{Float64,2}, me::AbstractMatrixExponential, x::Real) = ccdf(a,me)(x)

ccdf(me::AbstractMatrixExponential, x::Array{Float64}) = ccdf(me).(x)
ccdf(a::Array{Float64,2}, me::AbstractMatrixExponential, x::Array{Float64}) = ccdf(a,me).(x)

cdf(me::AbstractMatrixExponential) = x->1-ccdf(me,x)
"""
    cdf(a::Array{Float64, 2}, me::AbstractMatrixExponential) = begin

Return a the cumulative distribution function of a matrix exponential as a function.

If `a` is not specified then the initial vector of the `me` is `me.a`

# Arguments:
- `a`: a row vector with length `size(me.S,1)`
"""
cdf(a::Array{Float64,2}, me::AbstractMatrixExponential) = (length(a)==size(me.S,1)) ? (x->1-ccdf(a,me)(x)) : throw(
    DomainError("a and me.S must have compatible size"))

cdf(me::AbstractMatrixExponential, x::Real) = cdf(me)(x)
"""
    cdf(a::Array{Float64, 2}, me::AbstractMatrixExponential, x::Real)

Return a the cumulative distribution function of a matrix exponential evaluated at `x`.
"""
cdf(a::Array{Float64,2}, me::AbstractMatrixExponential, x::Real) = cdf(a,me)(x)

cdf(me::AbstractMatrixExponential, x::Array{Float64}) = cdf(me).(x)
cdf(a::Array{Float64,2}, me::AbstractMatrixExponential, x::Array{Float64}) = cdf(a,me).(x)

"""
    build_me(params; mean::Real = 1)

Return a ConcentratedMatrixExponential as defined by the dictionary params. i.e. see CMEParams[3]
for an order 3 CME.
"""
function build_me(params; mean::Real = 1)
    N = 2*params["n"]+1
    α = zeros(1,N)
    α[1] = params["c"]
    a = params["a"]
    b = params["b"]
    ω =  params["omega"]
    for k in 1:params["n"]
        kω = k*ω
        α[2*k] = (1/2)*( a[k]*(1+kω) - b[k]*(1-kω) )/(1+kω^2)
        α[2*k+1] = (1/2)*( a[k]*(1-kω) + b[k]*(1+kω) )/(1+kω^2)
    end
    α = α./sum(α)
    Q = zeros(N,N)
    Q[1,1] = -1
    for k in 1:params["n"]
        kω = k*ω
        idx = 2*k:(2*k+1)
        Q[idx,idx] = [-1 -kω; kω -1]
    end
    Q = CMEOperator(Q.*-sum(α/Q)./mean)
    q = -sum(Q,dims=2)[:]
    return ConcentratedMatrixExponential(α,Q,q,params["D"])
end

ConcentratedMatrixExponential(order::Int; mean::Float64 = 1.0) =
    build_me(cme_params[order], mean=mean)

function fast_exp(S::Matrix)

end

function build_erlang(order; mean::Float64 = 1.0)
    α = zeros(1,order) # inital distribution
    α[1] = 1
    λ = order/mean
    Q = zeros(order,order)
    Q = Q + LinearAlgebra.diagm(0=>repeat(-[λ],order), 1=>repeat([λ],order-1))
    q = -sum(Q,dims=2)[:]
    D = Array{Float64}(LinearAlgebra.I(order))[end:-1:1,:]
    return MatrixExponential(α,Q,q,D)
end

function orbit(me::AbstractMatrixExponential)
    function _orbit(t)
        num = me.a*exp(me.S*t)
        denom = sum(num)
        return num./denom
    end
    return t->_orbit(t)
end
orbit(me::AbstractMatrixExponential,t::Float64) =
    orbit(me)(t)

function orbit(me::ConcentratedMatrixExponential)
    params = cme_params[_order(me)]
    mean = -sum(me.a/me.S)
    function _orbit(t)
        position = zeros(size(me.a))
        position[1] = me.a[1]
        for k in 1:params["n"]
            kωt = k*params["omega"]*t
            idx = 2*k
            idx2 = idx+1
            position[idx] = me.a[idx]*cos(kωt) + me.a[idx2]*sin(kωt)
            position[idx2] = -me.a[idx]*sin(kωt) + me.a[idx2]*cos(kωt)
        end
        position = position./sum(position)
        return position
    end
    mu1 = params["mu1"]
    return t->_orbit(mu1*t/mean)
end

function expected_orbit_from_pdf(pdf::Function,me::AbstractMatrixExponential,a::Float64,b::Float64,evals::Int=10)
    # evals is an integer specifying how many points to eval the function at
    # params is a cme_params dictionary entry, i.e. cme_params[3]
    (b<=a)&&throw(DomainError("must have b>a"))
    delta = b-a # the orbit repeats after this time
    edges = range(0,delta,length=evals+1) # points at which to evaluate the fn
    h = delta/evals

    orbit_LHS = me.a
    pdf_LHS = pdf(a)

    E_orbit = zeros(size(me.a))
    for t in edges[2:end]
        orbit_RHS = orbit(me,t)
        orbit_estimate = (orbit_LHS+orbit_RHS)./2.0

        pdf_RHS = pdf(a+t)
        prob_estimate = (pdf_RHS+pdf_LHS)/2.0*h

        orbit_LHS = copy(orbit_RHS)
        pdf_LHS = copy(pdf_RHS)

        E_orbit += orbit_estimate*prob_estimate
    end
    return E_orbit
end

function expected_orbit_from_cdf(cdf::Function,me::AbstractMatrixExponential,a::Float64,b::Float64,evals::Int=10)
    # evals is an integer specifying how many points to eval the function at
    # params is a cme_params dictionary entry, i.e. cme_params[3]
    (b<=a)&&throw(DomainError("must have b>a"))
    delta = b-a # the orbit repeats after this time
    edges = range(0,delta,length=evals+1) # points at which to evaluate the fn

    orbit_LHS = me.a
    cdf_LHS = cdf(a)

    E_orbit = zeros(size(me.a))
    for t in edges[2:end]
        orbit_RHS = orbit(me,t)
        orbit_estimate = (orbit_LHS+orbit_RHS)./2

        cdf_RHS = cdf(a+t)
        prob_estimate = (cdf_RHS-cdf_LHS)

        orbit_LHS = copy(orbit_RHS)
        cdf_LHS = copy(cdf_RHS)

        E_orbit += orbit_estimate*prob_estimate
    end
    return E_orbit
end

## not used, but kept them because they might be / are cool
# function renewalProperties(me::AbstractMatrixExponential)
#     density(t) = begin
#         Q = me.S
#         q = me.s
#         α = me.a
#         e = ones(size(q))
#         (α*exp((Q+q*α)*t)*q)[1]
#     end
#     mean(t) = begin
#         Q = me.S
#         q = me.s
#         α = me.a
#         e = ones(size(q))
#         temp1 = α*-Q^-1
#         temp2 = temp1*e
#         temp3 = temp1./temp2
#         temp4 = Q + (q*α)
#         ((t./temp2) - α*(I - exp(temp4*t))*(temp4 + e*temp3)^-1*q)[1]
#     end
#     ExpectedOrbit(t) = begin
#         Q = me.S
#         q = me.s
#         α = me.a
#         e = ones(size(q))
#         (α*exp((Q + q*α)*t))[1]
#     end
#     return (density=density,mean=mean,ExpectedOrbit=ExpectedOrbit)
# end


text:
a = collect(1:20)
b = shuffle(a)

print("Unshuffled Array:\n    ")
println(a)
print("Shuffled Array:\n    ")
println(b)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using DataDrivenDiffEq
using OrdinaryDiffEq
using Plots
gr()

# Define measurements from unstable system with known control input
X = [4 2 1 0.5 0.25; 7 0.7 0.07 0.007 0.0007]
U = [-4 -2 -1 -0.5]
B = Float32[1; 0]

# See fail with unknown input
sys = DMDc(X, U)
# But with a little more knowledge
sys = DMDc(X, U, B = B)

# Extract the DMD from inside DMDc
get_dynamics(sys)
# Acess all the other stuff
eigen(sys)
eigvals(sys)
eigvecs(sys)
isstable(sys)
# Get unforced dynamics
dudt_ = dynamics(sys)
prob = DiscreteProblem(dudt_, X[:, 1], (0., 10.))
sol_unforced = solve(prob,  FunctionMap())
plot(sol_unforced)
sol_unforced[:,:]
# Create a system with cos control input to stabilize
dudt_ = dynamics(sys, control = (u, p, t) -> -0.5u[1])
prob = DiscreteProblem(dudt_, X[:, 1], (0., 10.))
sol = solve(prob, FunctionMap())

plot!(sol)


text:
# Wodarz, D., Nowak, M.
# Specific therapy regimes could lead to long-term immunological control of HIV
# https://doi.org/10.1073/pnas.96.25.14464
# Page 1
using Logging

using StructuralIdentifiability

logger = Logging.SimpleLogger(stdout, Logging.Info)
global_logger(logger)

ode = @ODEmodel(
    x'(t) = lm - d * x(t) - beta * x(t) * v(t),
    y'(t) = beta * x(t) * v(t) - a * y(t),
    v'(t) = k * y(t) - u * v(t),
    w'(t) = c * x(t) * y(t) * w(t) - c * q * y(t) * w(t) - b * w(t),
    z'(t) = c * q * y(t) * w(t) - h * z(t),
    y1(t) = w(t),
    y2(t) = z(t)
)

@time println(assess_identifiability(ode))


text:
################################################################################
## Copyright 2016. Juri Hinz and Jeremy Yee.
## Conditional expectation matrices for fast Bellman recursion
################################################################################

## Type containing the skeleton matrix and permutation matrices

immutable PermMat

    ## Matrices
    constMat::Array{Float64} # (gnum, gnum) non-random component of disturbances
    permMat::Array{Float64, 3} # (gnum, gnum, rnum + 1) permutation matrices
    rIndex::Array{Int64, 2} # location of random entries (row, column)
    rnum::Int64 # number of random entries in disturbances

    function PermMat(css::CSS, rIndex::Array{Int64, 2}, search::Function)

        ## Create search structure
        transgrid = transpose(css.grid)
        tree = search(transgrid)

        ## Create the matrices
        rnum = size(rIndex)[1];
        permMat = zero(Array{Float64}(css.gnum, css.gnum, rnum + 1))

        ## Assigning entries to matrices
        for i = 1:css.dnum
            nearest = knn(tree, css.disturb[:, :, i] * transgrid, 1, true)[1]
            for j = 1:css.gnum
                permMat[j, nearest[j], rnum + 1] +=  css.weight[i]
                for k = 1:rnum
                    permMat[j, nearest[j], k] +=
                        css.disturb[rIndex[k, 1], rIndex[k, 2], i] * css.weight[i]
                end
            end
        end

        ## Extract the disturbance skeleton
        constMat = css.disturb[:, :, 1]
        for i = 1:rnum
            constMat[rIndex[i, 1], rIndex[i, 2]] = 0
        end

        return new(constMat, permMat, rIndex, rnum)

    end

end


text:
"""
   Traits for a theoretical spectrograph
Author: Eric Ford
Created: August 2020
"""

""" Delegates loading of code specifying types essential to the package.  """

import ..RvSpectMLBase.InstrumentsCommon: min_pixel, max_pixel, min_order, max_order, min_pixel_in_order, max_pixel_in_order, get_pixel_range
import ..RvSpectMLBase.InstrumentsCommon: orders_to_use_default, min_col_default, max_col_default

min_pixel(::TheoreticalInstrument1D) = 1
max_pixel(::TheoreticalInstrument1D) = 128*8192
#min_order(::TheoreticalInstrument1D) = 1
#max_order(::TheoreticalInstrument1D) = 1
min_pixel_in_order(::TheoreticalInstrument1D) = 1
max_pixel_in_order(::TheoreticalInstrument1D) = 8192
#orders_to_use_default(inst::TheoreticalInstrument1D) = 1
#min_col_default(::TheoreticalInstrument1D ) = 1
#max_col_default(inst::TheoreticalInstrument1D ) = max_pixel_in_order(inst)
#min_col_default(::TheoreticalInstrument1D, ord::Integer) = min_col_default(inst)
#max_col_default(inst::TheoreticalInstrument1D, ord::Integer) = max_col_default(inst)

min_order(::TheoreticalInstrument2D) = 1
max_order(inst::TheoreticalInstrument2D) = length(inst.λ_min)
min_pixel(::TheoreticalInstrument2D) = 1
max_pixel(inst::TheoreticalInstrument2D) = inst.pixels_per_order*length(inst.λ_min)
min_pixel_in_order(::TheoreticalInstrument2D) = 1
max_pixel_in_order(inst::TheoreticalInstrument2D) = inst.pixels_per_order

orders_to_use_default(inst::TheoreticalInstrument2D) = min_order(inst):max_order(inst)
min_col_default(::TheoreticalInstrument2D, ord::Integer) = 1
max_col_default(inst::TheoreticalInstrument2D, ord::Integer) = max_pixel_in_order(inst)

function get_pixel_range(inst::TheoreticalInstrument2D, ord::Integer)
    min_col_default(inst, ord):max_col_default(inst, ord)
end

import ..RvSpectMLBase.InstrumentsCommon: metadata_symbols_default, metadata_strings_default
metadata_symbols_default(::AnyTheoreticalInstrument) = Symbol[] #:bjd, :target, :ssbz]
metadata_strings_default(::AnyTheoreticalInstrument) = String[] #"OBSJD", "SKY-OBJ", "SSBZ000"]

import ..RvSpectMLBase.InstrumentsCommon: default_ccf_mask_v_width
default_ccf_mask_v_width(::AnyTheoreticalInstrument) = 500.0  #
export default_ccf_mask_v_width

import ..RvSpectMLBase.InstrumentsCommon: get_inst_module
get_inst_module(::AnyTheoreticalInstrument) = TheoreticalInstrument
export get_inst_module

#import ..RvSpectMLBase: filter_line_list, find_worst_telluric_in_each_chunk
# No tellurics in theoretical spectra
#export filter_line_list, find_worst_telluric_in_each_chunk

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
__precompile__()

module MixedModels

using Compat, DataArrays, GLM, DataFrames, Distributions, NLopt, Showoff, StatsBase
using StatsFuns: log2π

import StatsBase: coef, coeftable, df, deviance, fit!, fitted, loglikelihood,
    model_response, nobs, vcov
import Base: cond, std
import Distributions: Bernoulli, Binomial, Poisson, Gamma
import GLM: LogitLink, LogLink, InverseLink
import DataFrames: @~

export
       @~,
       Bernoulli,
       Binomial,
       Poisson,
       Gamma,
       LogitLink,
       LogLink,
       InverseLink,
       GeneralizedLinearMixedModel,
       LinearMixedModel,
       MixedModel,
       OptSummary,
       ReMat,
       ScalarReMat,
       VarCorr,
       VectorReMat,

       bootstrap!,
       cfactor!,
       coef,
       coeftable,
       cond,
       df,
       deviance,
       fit!,
       fitted,
       fixef,      # extract the fixed-effects parameter estimates
       getθ,
       glmm,       # define a GeneralizedLinearMixedModel
       LaplaceDeviance, # Laplace approximation to GLMM deviance
       lmm,        # create a LinearMixedModel from a formula/data specification
       loglikelihood,
       lowerbd,    # lower bounds on the covariance parameters
       model_response,
       nobs,
       objective,  # the objective function in fitting a model
       pwrss,      # penalized, weighted residual sum-of-squares
       pirls!,     # use Penalized Iteratively Reweighted Least Squares to obtain conditional modes of random effects
       ranef,      # extract the conditional modes of the random effects
       refit!,     # install a response and refit the model
       remat,      # factory for construction of ReMat objects
       sdest,      # the estimate of the standard deviation of the per-observation noise
       setθ!,
       simulate!,  # simulate a new response in place
       std,
       varest,     # estimate of the residual variance
       vcov

abstract MixedModel <: RegressionModel # model with fixed and random effects

import Base: ==, *

include("blockmats.jl")
include("linalg.jl")
include("cfactor.jl")
include("remat.jl")
include("paramlowertriangular.jl")
include("inject.jl")
include("pls.jl")
include("logdet.jl")
include("simulate.jl")
include("PIRLS.jl")
include("glm.jl")
include("mixedmodel.jl")

end # module


text:
using StatsBase

"""
Normalize in logspace while avoiding numerical underflow
L is the log normalization constant

    logPost, L = normalizeLogspace(logprior .+ loglik)
    post = exp(logPost)

"""
function normalizeLogspace{T<:Real}(x::AbstractArray{T,1})
    L = logsumexp(x)
    y = broadcast(-, x, L)
    (y, L)
end


text:
#Demo 11: Measure frequency/period
#-------------------------------------------------------------------------------

using CMDimCircuits
CMDimCircuits.@using_CData()

#Get a demo display:
include(CMDimCircuits.demoplotcfgscript)
#Normally use something like:
#CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
LBLAX_AMPV = "Amplitude (V)"
LBLAX_PERNS = "Period (ns)"
LBLAX_FREQ = "Frequency (Hz)"
LBLAX_TIME = "Time (s)"
line1 = cons(:a, line = set(style=:solid))
lshiftattr = cons(:a,
	line = set(style=:solid, color=2), glyph = set(shape=:+),
)
lnoshiftattr = cons(:a,
	line = set(style=:solid, color=3), glyph = set(shape=:x),
)

#==Input data
===============================================================================#
fnom = 1e9
osr = 50 #(fundamental) oversampling ratio
ncycles = 20
fshift=0.1 #ratio
fmod = fnom/ncycles #Frequency of modulation
xrise = CrossType(:rise)

#==Computations
===============================================================================#
T = 1/fnom
t = DataF1(0:(T/osr):(ncycles*T))
finst = fnom*(1+fshift*cos(2pi*fmod*t)) #instantaneous frequency
y = cos(2pi*finst*t)

period_shift = measperiod(y, xing=xrise, shiftx=true)/1e-9
period_noshift = measperiod(y, xing=xrise, shiftx=false)/1e-9
freq_shift = measfreq(y, xing=xrise, shiftx=true)
freq_noshift = measfreq(y, xing=xrise, shiftx=false)

#==Generate plot
===============================================================================#
plot = cons(:plot, nstrips=3,
	ystrip1 = set(axislabel=LBLAX_AMPV, striplabel="Signal"),
	ystrip2 = set(axislabel=LBLAX_PERNS, striplabel="Instantaneous Period"),
	ystrip3 = set(axislabel=LBLAX_FREQ, striplabel="Instantaneous Frequency"),
	xaxis = set(label=LBLAX_TIME),
)
push!(plot,
	cons(:wfrm, y, label="y(t)", strip=1),
	cons(:wfrm, period_shift, lshiftattr, label="shift", strip=2),
	cons(:wfrm, period_noshift, lnoshiftattr, label="no shift", strip=2),
	cons(:wfrm, freq_shift, lshiftattr, label="shift", strip=3),
	cons(:wfrm, freq_noshift, lnoshiftattr, label="no shift", strip=3),
)

pcoll = push!(cons(:plotcoll, title="Period & Frequency Measurements"), plot)
	pcoll.ncolumns = 1

#==Display results in pcoll
===============================================================================#
EasyPlot.displaygui(pcoll)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll