text:
import PredictMD
import StatsBase
import Test

Test.@test_throws(ArgumentError, PredictMD.get_number_in_each_fold(1, 2))

Test.@test_throws(DimensionMismatch, PredictMD.get_top_level_num_folds(PredictMD.CrossValidation{Int}(PredictMD.CrossValidation{Int}[PredictMD.CrossValidation{Int}(; all_indices = [1,2,3,4,5,6,7,8,9,10], num_folds_per_level = (2,2,2))], Vector{Int}[])))

Test.@test(PredictMD.get_top_level_num_folds(PredictMD.CrossValidation{Int}(PredictMD.CrossValidation{Int}[], Vector{Int}[])) == 0)

Test.@testset "CrossValidation" begin
    Test.@testset "get_number_in_each_fold" begin
        @debug("get_number_in_each_fold")
        Test.@test( sum(PredictMD.get_number_in_each_fold(100, 1)) == 100 )
        Test.@test( PredictMD.get_number_in_each_fold(100, 1) == [100] )
        Test.@test( sum(PredictMD.get_number_in_each_fold(100, 2)) == 100 )
        Test.@test( PredictMD.get_number_in_each_fold(100, 2) == [50, 50] )
        Test.@test( sum(PredictMD.get_number_in_each_fold(100, 3)) == 100 )
        Test.@test( PredictMD.get_number_in_each_fold(100, 3) == [34, 33, 33] )
        Test.@test( sum(PredictMD.get_number_in_each_fold(100, 4)) == 100 )
        Test.@test( PredictMD.get_number_in_each_fold(100, 4) == [25, 25, 25, 25] )
        Test.@test( sum(PredictMD.get_number_in_each_fold(100, 5)) == 100 )
        Test.@test( PredictMD.get_number_in_each_fold(100, 5) == [20, 20, 20, 20, 20] )
        Test.@test( sum(PredictMD.get_number_in_each_fold(100, 6)) == 100 )
        Test.@test( PredictMD.get_number_in_each_fold(100, 6) == [17, 17, 17, 17, 16, 16] )
        Test.@test( sum(PredictMD.get_number_in_each_fold(100, 7)) == 100 )
        Test.@test( PredictMD.get_number_in_each_fold(100, 7) == [15, 15, 14, 14, 14, 14, 14] )
        Test.@test( sum(PredictMD.get_number_in_each_fold(100, 8)) == 100 )
        Test.@test( PredictMD.get_number_in_each_fold(100, 8) == [13, 13, 13, 13, 12, 12, 12, 12] )
        Test.@test( sum(PredictMD.get_number_in_each_fold(100, 9)) == 100 )
        Test.@test( PredictMD.get_number_in_each_fold(100, 9) == [12, 11, 11, 11, 11, 11, 11, 11, 11] )
        Test.@test( sum(PredictMD.get_number_in_each_fold(100, 10)) == 100 )
        Test.@test( PredictMD.get_number_in_each_fold(100, 10) == [10, 10, 10, 10, 10, 10, 10, 10, 10, 10] )
    end

    Test.@testset "get_indices_in_each_fold" begin
        @debug("get_indices_in_each_fold")
        Test.@test( PredictMD.get_indices_in_each_fold([1,2,3,4,5,6,7,8,9,10], 3) == [[1,2,3,4], [5,6,7], [8,9,10]] )
        Test.@test( PredictMD.get_indices_in_each_fold([5,6,7,8,9,10], 3) == [[5,6], [7,8], [9,10]] )
        Test.@test( PredictMD.get_indices_in_each_fold([1,2,3,4,8,9,10], 3) == [[1,2,3], [4,8], [9,10]] )
        Test.@test( PredictMD.get_indices_in_each_fold([1,2,3,4,5,6,7], 3) == [[1,2,3], [4,5], [6,7]] )
    end

    Test.@testset "get_leavein_indices on vector of integers" begin
        @debug("get_leavein_indices")
        Test.@test( PredictMD.get_leavein_indices([1,2,3,4,5,6,7,8,9,10], 3, 1) == [5,6,7,8,9,10] )
        Test.@test( PredictMD.get_leavein_indices([1,2,3,4,5,6,7,8,9,10], 3, 2) == [1,2,3,4,8,9,10] )
        Test.@test( PredictMD.get_leavein_indices([1,2,3,4,5,6,7,8,9,10], 3, 3) == [1,2,3,4,5,6,7] )
    end

    Test.@testset "get_leaveout_indices on vector of integers" begin
        @debug("get_leaveout_indices")
        Test.@test( PredictMD.get_leaveout_indices([1,2,3,4,5,6,7,8,9,10], 3, 1) == [1,2,3,4] )
        Test.@test( PredictMD.get_leaveout_indices([1,2,3,4,5,6,7,8,9,10], 3, 2) == [5,6,7] )
        Test.@test( PredictMD.get_leaveout_indices([1,2,3,4,5,6,7,8,9,10], 3, 3) == [8,9,10] )
    end

    Test.@testset "nested cross validation, integer indices, small" begin
        @debug("nested cross validation, integer indices, small")
        cv = PredictMD.CrossValidation{Int}(; all_indices = [1,2,3,4,5,6,7,8,9,10], num_folds_per_level = (2,2,2))
        Test.@test( isa(cv, PredictMD.CrossValidation{Int}) )
        Test.@test( !PredictMD.isleaf(cv) )
        Test.@test( PredictMD.get_top_level_num_folds(cv) == 2 )
        Test.@test( PredictMD.get_leaveout_indices(cv, 1) == [1,2,3,4,5] )
        Test.@test( PredictMD.get_leaveout_indices(cv, 2) == [6,7,8,9,10] )
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_indices(cv) )

        cv_1 = PredictMD.get_leavein_cv(cv, 1)
        Test.@test( !PredictMD.isleaf(cv_1) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_1) == 2 )
        Test.@test( PredictMD.get_leaveout_indices(cv_1, 1) == [6,7,8] )
        Test.@test( PredictMD.get_leaveout_indices(cv_1, 2) == [9,10] )
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_indices(cv_1) )

        cv_2 = PredictMD.get_leavein_cv(cv, 2)
        Test.@test( !PredictMD.isleaf(cv_2) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_2) == 2 )
        Test.@test( PredictMD.get_leaveout_indices(cv_2, 1) ==  [1,2,3])
        Test.@test( PredictMD.get_leaveout_indices(cv_2, 2) ==  [4,5])
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_indices(cv_2) )

        cv_1_1 = PredictMD.get_leavein_cv(cv_1, 1)
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_indices(cv_1_1) )
        Test.@test( !PredictMD.isleaf(cv_1_1) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_1_1) == 2 )
        Test.@test( PredictMD.get_leaveout_indices(cv_1_1, 1) == [9] )
        Test.@test( PredictMD.get_leaveout_indices(cv_1_1, 2) == [10] )

        cv_1_2 = PredictMD.get_leavein_cv(cv_1, 2)
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_indices(cv_1_2) )
        Test.@test( !PredictMD.isleaf(cv_1_2) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_1_2) == 2 )
        Test.@test( PredictMD.get_leaveout_indices(cv_1_2, 1) == [6, 7] )
        Test.@test( PredictMD.get_leaveout_indices(cv_1_2, 2) == [8] )

        cv_2_1 = PredictMD.get_leavein_cv(cv_2, 1)
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_indices(cv_2_1) )
        Test.@test( !PredictMD.isleaf(cv_2_1) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_2_1) == 2 )
        Test.@test( PredictMD.get_leaveout_indices(cv_2_1, 1) == [4] )
        Test.@test( PredictMD.get_leaveout_indices(cv_2_1, 2) == [5] )

        cv_2_2 = PredictMD.get_leavein_cv(cv_2, 2)
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_indices(cv_2_2) )
        Test.@test( !PredictMD.isleaf(cv_2_2) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_2_2) == 2 )
        Test.@test( PredictMD.get_leaveout_indices(cv_2_2, 1) == [1,2] )
        Test.@test( PredictMD.get_leaveout_indices(cv_2_2, 2) == [3] )

        cv_1_1_1 = PredictMD.get_leavein_cv(cv_1_1, 1)
        Test.@test( PredictMD.isleaf(cv_1_1_1) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_1_1_1) == 0 )
        Test.@test_throws( ArgumentError, PredictMD.get_leaveout_indices(cv_1_1_1, 1) )
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_cv(cv_1_1_1, 1) )
        Test.@test( PredictMD.get_leavein_indices(cv_1_1_1) == [10] )

        cv_1_1_2 = PredictMD.get_leavein_cv(cv_1_1, 2)
        Test.@test( PredictMD.isleaf(cv_1_1_2) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_1_1_2) == 0 )
        Test.@test_throws( ArgumentError, PredictMD.get_leaveout_indices(cv_1_1_2, 1) )
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_cv(cv_1_1_2, 1) )
        Test.@test( PredictMD.get_leavein_indices(cv_1_1_2) == [9] )

        cv_1_2_1 = PredictMD.get_leavein_cv(cv_1_2, 1)
        Test.@test( PredictMD.isleaf(cv_1_2_1) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_1_2_1) == 0 )
        Test.@test_throws( ArgumentError, PredictMD.get_leaveout_indices(cv_1_2_1, 1) )
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_cv(cv_1_2_1, 1) )
        Test.@test( PredictMD.get_leavein_indices(cv_1_2_1) == [8] )

        cv_1_2_2 = PredictMD.get_leavein_cv(cv_1_2, 2)
        Test.@test( PredictMD.isleaf(cv_1_2_2) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_1_2_2) == 0 )
        Test.@test_throws( ArgumentError, PredictMD.get_leaveout_indices(cv_1_2_2, 1) )
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_cv(cv_1_2_2, 1) )
        Test.@test( PredictMD.get_leavein_indices(cv_1_2_2) == [6,7] )

        cv_2_1_1 = PredictMD.get_leavein_cv(cv_2_1, 1)
        Test.@test( PredictMD.isleaf(cv_2_1_1) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_2_1_1) == 0 )
        Test.@test_throws( ArgumentError, PredictMD.get_leaveout_indices(cv_2_1_1, 1) )
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_cv(cv_2_1_1, 1) )
        Test.@test( PredictMD.get_leavein_indices(cv_2_1_1) == [5] )

        cv_2_1_2 = PredictMD.get_leavein_cv(cv_2_1, 2)
        Test.@test( PredictMD.isleaf(cv_2_1_2) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_2_1_2) == 0 )
        Test.@test_throws( ArgumentError, PredictMD.get_leaveout_indices(cv_2_1_2, 1) )
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_cv(cv_2_1_2, 1) )
        Test.@test( PredictMD.get_leavein_indices(cv_2_1_2) == [4] )

        cv_2_2_1 = PredictMD.get_leavein_cv(cv_2_2, 1)
        Test.@test( PredictMD.isleaf(cv_2_2_1) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_2_2_1) == 0 )
        Test.@test_throws( ArgumentError, PredictMD.get_leaveout_indices(cv_2_2_1, 1) )
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_cv(cv_2_2_1, 1) )
        Test.@test( PredictMD.get_leavein_indices(cv_2_2_1) == [3] )

        cv_2_2_2 = PredictMD.get_leavein_cv(cv_2_2, 2)
        Test.@test( PredictMD.isleaf(cv_2_2_2) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_2_2_2) == 0 )
        Test.@test_throws( ArgumentError, PredictMD.get_leaveout_indices(cv_2_2_2, 1) )
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_cv(cv_2_2_2, 1) )
        Test.@test( PredictMD.get_leavein_indices(cv_2_2_2) == [1,2] )
    end

    Test.@testset "nested cross validation, range indices, small" begin
        @debug("nested cross validation, range indices, small")
        cv_integers = PredictMD.CrossValidation{Int}(; all_indices = [1,2,3,4,5,6,7,8,9,10], num_folds_per_level = (2,2,2))
        Test.@test( isa(cv_integers, PredictMD.CrossValidation{Int}) )
        cv_ranges = PredictMD.CrossValidation{UnitRange{Int}}(cv_integers)
        Test.@test( isa(cv_ranges, PredictMD.CrossValidation{UnitRange{Int}}) )
        Test.@test( !PredictMD.isleaf(cv_ranges) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_ranges) == 2 )
        Test.@test( PredictMD.get_leaveout_indices(cv_ranges, 1) == [1:5] )
        Test.@test( PredictMD.get_leaveout_indices(cv_ranges, 2) == [6:10] )
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_indices(cv_ranges) )

        cv_1 = PredictMD.get_leavein_cv(cv_ranges, 1)
        Test.@test( !PredictMD.isleaf(cv_1) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_1) == 2 )
        Test.@test( PredictMD.get_leaveout_indices(cv_1, 1) == [6:8] )
        Test.@test( PredictMD.get_leaveout_indices(cv_1, 2) == [9:10] )
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_indices(cv_1) )

        cv_2 = PredictMD.get_leavein_cv(cv_ranges, 2)
        Test.@test( !PredictMD.isleaf(cv_2) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_2) == 2 )
        Test.@test( PredictMD.get_leaveout_indices(cv_2, 1) ==  [1:3])
        Test.@test( PredictMD.get_leaveout_indices(cv_2, 2) ==  [4:5])
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_indices(cv_2) )

        cv_1_1 = PredictMD.get_leavein_cv(cv_1, 1)
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_indices(cv_1_1) )
        Test.@test( !PredictMD.isleaf(cv_1_1) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_1_1) == 2 )
        Test.@test( PredictMD.get_leaveout_indices(cv_1_1, 1) == [9:9] )
        Test.@test( PredictMD.get_leaveout_indices(cv_1_1, 2) == [10:10] )

        cv_1_2 = PredictMD.get_leavein_cv(cv_1, 2)
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_indices(cv_1_2) )
        Test.@test( !PredictMD.isleaf(cv_1_2) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_1_2) == 2 )
        Test.@test( PredictMD.get_leaveout_indices(cv_1_2, 1) == [6:7] )
        Test.@test( PredictMD.get_leaveout_indices(cv_1_2, 2) == [8:8] )

        cv_2_1 = PredictMD.get_leavein_cv(cv_2, 1)
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_indices(cv_2_1) )
        Test.@test( !PredictMD.isleaf(cv_2_1) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_2_1) == 2 )
        Test.@test( PredictMD.get_leaveout_indices(cv_2_1, 1) == [4:4] )
        Test.@test( PredictMD.get_leaveout_indices(cv_2_1, 2) == [5:5] )

        cv_2_2 = PredictMD.get_leavein_cv(cv_2, 2)
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_indices(cv_2_2) )
        Test.@test( !PredictMD.isleaf(cv_2_2) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_2_2) == 2 )
        Test.@test( PredictMD.get_leaveout_indices(cv_2_2, 1) == [1:2] )
        Test.@test( PredictMD.get_leaveout_indices(cv_2_2, 2) == [3:3] )

        cv_1_1_1 = PredictMD.get_leavein_cv(cv_1_1, 1)
        Test.@test( PredictMD.isleaf(cv_1_1_1) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_1_1_1) == 0 )
        Test.@test_throws( ArgumentError, PredictMD.get_leaveout_indices(cv_1_1_1, 1) )
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_cv(cv_1_1_1, 1) )
        Test.@test( PredictMD.get_leavein_indices(cv_1_1_1) == [10:10] )

        cv_1_1_2 = PredictMD.get_leavein_cv(cv_1_1, 2)
        Test.@test( PredictMD.isleaf(cv_1_1_2) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_1_1_2) == 0 )
        Test.@test_throws( ArgumentError, PredictMD.get_leaveout_indices(cv_1_1_2, 1) )
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_cv(cv_1_1_2, 1) )
        Test.@test( PredictMD.get_leavein_indices(cv_1_1_2) == [9:9] )

        cv_1_2_1 = PredictMD.get_leavein_cv(cv_1_2, 1)
        Test.@test( PredictMD.isleaf(cv_1_2_1) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_1_2_1) == 0 )
        Test.@test_throws( ArgumentError, PredictMD.get_leaveout_indices(cv_1_2_1, 1) )
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_cv(cv_1_2_1, 1) )
        Test.@test( PredictMD.get_leavein_indices(cv_1_2_1) == [8:8] )

        cv_1_2_2 = PredictMD.get_leavein_cv(cv_1_2, 2)
        Test.@test( PredictMD.isleaf(cv_1_2_2) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_1_2_2) == 0 )
        Test.@test_throws( ArgumentError, PredictMD.get_leaveout_indices(cv_1_2_2, 1) )
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_cv(cv_1_2_2, 1) )
        Test.@test( PredictMD.get_leavein_indices(cv_1_2_2) == [6:7] )

        cv_2_1_1 = PredictMD.get_leavein_cv(cv_2_1, 1)
        Test.@test( PredictMD.isleaf(cv_2_1_1) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_2_1_1) == 0 )
        Test.@test_throws( ArgumentError, PredictMD.get_leaveout_indices(cv_2_1_1, 1) )
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_cv(cv_2_1_1, 1) )
        Test.@test( PredictMD.get_leavein_indices(cv_2_1_1) == [5:5] )

        cv_2_1_2 = PredictMD.get_leavein_cv(cv_2_1, 2)
        Test.@test( PredictMD.isleaf(cv_2_1_2) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_2_1_2) == 0 )
        Test.@test_throws( ArgumentError, PredictMD.get_leaveout_indices(cv_2_1_2, 1) )
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_cv(cv_2_1_2, 1) )
        Test.@test( PredictMD.get_leavein_indices(cv_2_1_2) == [4:4] )

        cv_2_2_1 = PredictMD.get_leavein_cv(cv_2_2, 1)
        Test.@test( PredictMD.isleaf(cv_2_2_1) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_2_2_1) == 0 )
        Test.@test_throws( ArgumentError, PredictMD.get_leaveout_indices(cv_2_2_1, 1) )
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_cv(cv_2_2_1, 1) )
        Test.@test( PredictMD.get_leavein_indices(cv_2_2_1) == [3:3] )

        cv_2_2_2 = PredictMD.get_leavein_cv(cv_2_2, 2)
        Test.@test( PredictMD.isleaf(cv_2_2_2) )
        Test.@test( PredictMD.get_top_level_num_folds(cv_2_2_2) == 0 )
        Test.@test_throws( ArgumentError, PredictMD.get_leaveout_indices(cv_2_2_2, 1) )
        Test.@test_throws( ArgumentError, PredictMD.get_leavein_cv(cv_2_2_2, 1) )
        Test.@test( PredictMD.get_leavein_indices(cv_2_2_2) == [1:2] )
    end

    Test.@testset "nested cross validation, integer indices, large" begin
        @debug("nested cross validation, integer indices, large")
        num_samples = 1_000
        master_cv = PredictMD.CrossValidation{Int}(; all_indices = collect(1:num_samples), num_folds_per_level = (3,5,7))
        Test.@test( isa(master_cv, PredictMD.CrossValidation{Int}) )
        Test.@test( PredictMD.get_all_indices(master_cv) == 1:num_samples)
        Test.@test( !PredictMD.isleaf(master_cv) )
        Test.@test( PredictMD.get_top_level_num_folds(master_cv) == 3 )
        for i = 1:3
            @debug("i = $(i)")
            testing_indices = PredictMD.get_leaveout_indices(master_cv, i)
            supertuning_tuning_training_cv = PredictMD.get_leavein_cv(master_cv, i)
            Test.@test( length(intersect(testing_indices, PredictMD.get_all_indices(supertuning_tuning_training_cv))) == 0 )
            Test.@test( sort(unique(union(testing_indices, PredictMD.get_all_indices(supertuning_tuning_training_cv)))) == 1:num_samples )
            Test.@test( !PredictMD.isleaf(supertuning_tuning_training_cv) )
            Test.@test( PredictMD.get_top_level_num_folds(supertuning_tuning_training_cv) == 5 )
            for j = 1:5
                @debug("j = $(j)")
                supertuning_indices = PredictMD.get_leaveout_indices(supertuning_tuning_training_cv, j)
                tuning_training_cv = PredictMD.get_leavein_cv(supertuning_tuning_training_cv, j)
                Test.@test( !PredictMD.isleaf(tuning_training_cv) )
                Test.@test( PredictMD.get_top_level_num_folds(tuning_training_cv) == 7 )
                for k = 1:7
                    @debug("k = $(k)")
                    tuning_indices = PredictMD.get_leaveout_indices(tuning_training_cv, k)
                    training_cv = PredictMD.get_leavein_cv(tuning_training_cv, k)
                    Test.@test( PredictMD.isleaf(training_cv) )
                    Test.@test( PredictMD.get_top_level_num_folds(training_cv) == 0 )
                    training_indices = PredictMD.get_leavein_indices(training_cv)
                    Test.@test( length(intersect(testing_indices, supertuning_indices)) == 0 )
                    Test.@test( length(intersect(testing_indices, tuning_indices)) == 0 )
                    Test.@test( length(intersect(testing_indices, training_indices)) == 0 )
                    Test.@test( length(intersect(supertuning_indices, tuning_indices)) == 0 )
                    Test.@test( length(intersect(supertuning_indices, training_indices)) == 0 )
                    Test.@test( length(intersect(tuning_indices, training_indices)) == 0 )
                end
            end
        end
    end

    Test.@testset "nested cross validation, range indices, large" begin
        @debug("nested cross validation, range indices, large")
        num_samples = 1_000
        master_cv_integers = PredictMD.CrossValidation{Int}(; all_indices = collect(1:num_samples), num_folds_per_level = (3,5,7))
        Test.@test( isa(master_cv_integers, PredictMD.CrossValidation{Int}) )
        master_cv_ranges = PredictMD.CrossValidation{UnitRange{Int}}(master_cv_integers)
        Test.@test( isa(master_cv_ranges, PredictMD.CrossValidation{UnitRange{Int}}) )
        Test.@test( !PredictMD.isleaf(master_cv_ranges) )
        Test.@test( PredictMD.get_top_level_num_folds(master_cv_ranges) == 3 )
        for i = 1:3
            @debug("i = $(i)")
            testing_indices = PredictMD.get_leaveout_indices(master_cv_ranges, i)
            supertuning_tuning_training_cv = PredictMD.get_leavein_cv(master_cv_ranges, i)
            Test.@test( !PredictMD.isleaf(supertuning_tuning_training_cv) )
            Test.@test( PredictMD.get_top_level_num_folds(supertuning_tuning_training_cv) == 5 )
            for j = 1:5
                @debug("j = $(j)")
                supertuning_indices = PredictMD.get_leaveout_indices(supertuning_tuning_training_cv, j)
                tuning_training_cv = PredictMD.get_leavein_cv(supertuning_tuning_training_cv, j)
                Test.@test( !PredictMD.isleaf(tuning_training_cv) )
                Test.@test( PredictMD.get_top_level_num_folds(tuning_training_cv) == 7 )
                for k = 1:7
                    @debug("k = $(k)")
                    tuning_indices = PredictMD.get_leaveout_indices(tuning_training_cv, k)
                    training_cv = PredictMD.get_leavein_cv(tuning_training_cv, k)
                    Test.@test( PredictMD.isleaf(training_cv) )
                    Test.@test( PredictMD.get_top_level_num_folds(training_cv) == 0 )
                    training_indices = PredictMD.get_leavein_indices(training_cv)
                    Test.@test( length(intersect(testing_indices, supertuning_indices)) == 0 )
                    Test.@test( length(intersect(testing_indices, tuning_indices)) == 0 )
                    Test.@test( length(intersect(testing_indices, training_indices)) == 0 )
                    Test.@test( length(intersect(supertuning_indices, tuning_indices)) == 0 )
                    Test.@test( length(intersect(supertuning_indices, training_indices)) == 0 )
                    Test.@test( length(intersect(tuning_indices, training_indices)) == 0 )
                end
            end
        end
    end

    Test.@testset "roundtrip CV integer indices <-> CV range indices" begin
        @debug("roundtrip CV integer indices <-> CV range indices")
        num_samples = 1_000
        cv_integer_1 = PredictMD.CrossValidation{Int}(; all_indices = collect(1:num_samples), num_folds_per_level = (3,5,7))
        cv_ranges_2 = PredictMD.CrossValidation{UnitRange{Int}}(cv_integer_1)
        cv_integer_3 = PredictMD.CrossValidation{Int}(cv_ranges_2)
        cv_ranges_4 = PredictMD.CrossValidation{UnitRange{Int}}(cv_integer_3)
        cv_integer_5 = PredictMD.CrossValidation{Int}(cv_ranges_4)
        cv_ranges_6 = PredictMD.CrossValidation{UnitRange{Int}}(cv_integer_5)
        cv_integer_7 = PredictMD.CrossValidation{Int}(cv_ranges_6)
        cv_ranges_8 = PredictMD.CrossValidation{UnitRange{Int}}(cv_integer_7)
        cv_integer_9 = PredictMD.CrossValidation{Int}(cv_ranges_8)
        Test.@test( isa(cv_integer_1, PredictMD.CrossValidation{Int}) )
        Test.@test( isa(cv_integer_3, PredictMD.CrossValidation{Int}) )
        Test.@test( isa(cv_integer_5, PredictMD.CrossValidation{Int}) )
        Test.@test( isa(cv_integer_7, PredictMD.CrossValidation{Int}) )
        Test.@test( isa(cv_integer_9, PredictMD.CrossValidation{Int}) )
        Test.@test( isa(cv_ranges_2, PredictMD.CrossValidation{UnitRange{Int}}) )
        Test.@test( isa(cv_ranges_4, PredictMD.CrossValidation{UnitRange{Int}}) )
        Test.@test( isa(cv_ranges_6, PredictMD.CrossValidation{UnitRange{Int}}) )
        Test.@test( isa(cv_ranges_8, PredictMD.CrossValidation{UnitRange{Int}}) )
        Test.@test( cv_integer_1 == cv_integer_3 == cv_integer_5 == cv_integer_7 == cv_integer_9)
        Test.@test( cv_ranges_2 == cv_ranges_4 == cv_ranges_6 == cv_ranges_8)
    end

    Test.@testset "vectors_to_ranges" begin
        @debug("vectors_to_ranges")
        x = Int[
            12, 28, 27, 24, 21, 30, 6, 10, 4, 18,
            16, 36, 35, 29, 15, 9, 19, 34, 17, 5,
            23, 3, 26, 37, 20, 11, 7,
            ]
        y = UnitRange{Int}[
            3:7, 9:12, 15:21, 23:24, 26:30, 34:37,
            ]
        Test.@test y == PredictMD.vector_to_ranges(x)
        vector_1 = StatsBase.sample(1:100_000_000, 100_000)
        unique!(vector_1)
        sort!(vector_1)
        ranges_2 = PredictMD.vector_to_ranges(vector_1)
        vector_3 = PredictMD.ranges_to_vector(ranges_2)
        ranges_4 = PredictMD.vector_to_ranges(vector_3)
        vector_5 = PredictMD.ranges_to_vector(ranges_4)
        ranges_6 = PredictMD.vector_to_ranges(vector_5)
        vector_7 = PredictMD.ranges_to_vector(ranges_6)
        ranges_8 = PredictMD.vector_to_ranges(vector_7)
        vector_9 = PredictMD.ranges_to_vector(ranges_8)
        Test.@test vector_1 == vector_3 == vector_5 == vector_7 == vector_9
        Test.@test ranges_2 == ranges_4 == ranges_6 == ranges_8
    end

    Test.@testset "ranges_to_vectors" begin
        @debug("ranges_to_vectors")
        x = Int[
            12, 28, 27, 24, 21, 30, 6, 10, 4, 18,
            16, 36, 35, 29, 15, 9, 19, 34, 17, 5,
            23, 3, 26, 37, 20, 11, 7,
            ]
        y = UnitRange{Int}[
            3:7, 9:12, 15:21, 23:24, 26:30, 34:37,
            ]
        Test.@test sort(x) == PredictMD.ranges_to_vector(y)
        ranges_1 = Vector{UnitRange{Int}}(undef, 0)
        for i = 1:100
            a = StatsBase.sample((i)*(1_000_000):(i+1)*(1_000_000))
            b = StatsBase.sample((i)*(1_000_000):(i+1)*(1_000_000))
            push!(ranges_1, min(a,b):max(a,b))
        end
        unique!(ranges_1)
        sort!(ranges_1)
        vector_2 = PredictMD.ranges_to_vector(ranges_1)
        ranges_3 = PredictMD.vector_to_ranges(vector_2)
        vector_4 = PredictMD.ranges_to_vector(ranges_3)
        ranges_5 = PredictMD.vector_to_ranges(vector_4)
        vector_6 = PredictMD.ranges_to_vector(ranges_5)
        ranges_7 = PredictMD.vector_to_ranges(vector_6)
        vector_8 = PredictMD.ranges_to_vector(ranges_7)
        ranges_9 = PredictMD.vector_to_ranges(vector_8)
        Test.@test ranges_1 == ranges_3 == ranges_5 == ranges_7 == ranges_9
        Test.@test vector_2 == vector_4 == vector_6 == vector_8
    end
end


text:
#' #Pipeline with single subject (all conditions)
using Plots
gr()
using Measures
using MegTools
test = load_cont_epochs("condEpoch.mat")

test_average = average_across_trials(test)

#' Getting the averages of all conditions
test_average_auditory,auditory_left,auditory_right = select_channels(
    test_average,
    paradigm="auditoryN1m",
)

avs = plot(layout=(4,1), size=(1000,1000))
for (cond,val) in test_average_auditory
    plot!(avs, test_average_auditory[cond], legend=false)
end

avs

#' Filtering and baseline correcting
filtered = highlow_butterworth_filter(test_average_auditory, 1000)
baseline_corrected = baseline_correction(filtered)

peaks = find_peaks(baseline_corrected, auditory_left, auditory_right)

#' Loading trigger values from labels
cond_trigger_vals = load_trigger_values("regsoi")

#' For this specific analysis I need the peaks from all averaged and preprocecced ERFS
soi, left_amps, right_amps = collect_peaks(peaks)

#' Plotting peaks that we got
scatter(soi,left_amps, label="Left peak amplitudes");
scatter!(soi,right_amps, label="Right peak amplitudes")

ext: jl

max_stars_repo_name: ElectronicTeaCup/MegTools

max_issues_repo_name: ElectronicTeaCup/MegTools

max_forks_repo_name: ElectronicTeaCup/MegTools


text:
module Emporos

using MacroTools
using PyCall
using Unitful
using UnitfulAssets
using Dates

export IBKR

include("./util.jl")
include("./types.jl")

include("./IBKR/IBKR.jl")
using .IBKR

end


text:
module HexFFT

import Base: fft, ifft

export hfft2, ihfft2, fft, ifft, OffsetHexData

# HELPER FUNCTIONS
pad_zeros(arr) = hcat(arr, zeros(arr))
decimate(data, offset) = view(data, :, offset + 1:2:size(data, 2))
div2(n::Int) = n >> 1
function fold_half(data, op)
    n = size(data, 1)
    op.(view(data, 1:div2(n), :), view(data, div2(n) + 1:n, :))
end

# FORWARD TRANSFORM
function dft_nst1(data)
    out = fft(pad_zeros(data), 2)
    return decimate(out, 0), decimate(out, 1)
end

nst2(data) = repmat(fft(fold_half(data, (+)), 1), 2, 1)

coeff_inner(i, n) = exp(-2pi * im * i / n)
nst3_coeffs(n, m) = reshape(coeff_inner.(0:div2(n) - 1, n), div2(n), 1)
nst3(data) = repmat(fft(nst3_coeffs(size(data)...) .*
                        fold_half(data, (-)), 1), 2, 1)

w_inner(b, s, d, n, m) = exp(-pi * im * ((b + 2d) / 2m + (b + 2s) / n))
w(b, n, m) = [w_inner(b, s, d, n, m) for s in 0:n - 1, d in 0:m - 1]

function hfft2(data0, data1)
    g00, g01 = dft_nst1(data0)
    g10, g11 = dft_nst1(data1)
    return nst2(g00) .+ w(0, size(g10)...) .* nst2(g10),
           nst3(g01) .+ w(1, size(g11)...) .* nst3(g11)
end

# INVERSE TRANSFORM
function idft_inst1(data)
    out = 2ifft(pad_zeros(data), 2)
    return decimate(out, 0), decimate(out, 1)
end

inst2(data) = repmat(0.5ifft(fold_half(data, (+)), 1), 2, 1)

icoeff_inner(i, n) = exp(2pi * im * i / n)
inst3_coeffs(n, m) = reshape(icoeff_inner.(0:div2(n) - 1, n), div2(n), 1)
inst3(data) = repmat(0.5ifft(inst3_coeffs(size(data)...) .*
                             fold_half(data, (-)), 1), 2, 1)

iw_inner(a, r, c, n, m) = exp(pi * im * ((a + 2c) / 2m + (a + 2r) / n))
iw(a, n, m) = [iw_inner(a, r, c, n, m) for r in 0:n - 1, c in 0:m - 1]

function ihfft2(data0, data1)
    g00, g01 = idft_inst1(data0)
    g10, g11 = idft_inst1(data1)
    return 0.5(inst2(g00) .+ iw(0, size(g10)...) .* inst2(g10)),
           0.5(inst3(g01) .+ iw(1, size(g11)...) .* inst3(g11))
end

# Hex Array type
immutable OffsetHexData{T, S <: AbstractArray{T, 2}}
    odd_rows::S
    even_rows::S
end
function OffsetHexData{S}(odd_rows::S, even_rows::S)
    @assert size(odd_rows) == size(even_rows)
    OffsetHexData{eltype(odd_rows), S}(odd_rows, even_rows)
end

fft(H::OffsetHexData) = OffsetHexData(hfft2(H.odd_rows, H.even_rows)...)
ifft(H::OffsetHexData) = OffsetHexData(ihfft2(H.odd_rows, H.even_rows)...)

end # module


text:
module MLBStatsAPI

using HTTP, JSON
include("endpoints.jl")

export game, schedulemlb, attendance, awards, draft, people, standings, teams, meta

# --- private functions ---

function __formaturl(apiinfo, params)
    url = apiinfo["url"]
    pathparams = apiinfo["pathparams"]
    for param in pathparams
        paramval = nothing
        if haskey(params, param)
            paramval = params[param]
        else
            if haskey(apiinfo, param)
                paramval = apiinfo[param]
            else
                throw(ArgumentError("Required Path Parameters $pathparams not found."))
            end
        end
        url = replace(url, "$param" => "$paramval")
    end
    return url
end

function __genquery(apiinfo, params)
    paramkeys = intersect(apiinfo["params"], keys(params))
    return Dict(key => params[key] for key in paramkeys)
end

function __executeapi(apiname, params)
    apiinfo = ENDPOINTS[apiname]
    required = apiinfo["required"]
    if length(required) != 0 && length(intersect(required, keys(params))) == 0
        throw(ArgumentError("$required must be filled in."))
    end
    url = __formaturl(apiinfo, params)
    query = __genquery(apiinfo, params)
    res = HTTP.request("GET", url, query = query, status_exception = false)
    @info "HTTP Status Code: " * string(res.status)
    if res.status == 200
        return JSON.parse(String(res.body))
    else
        throw(ErrorException("No data."))
    end
end

# --- wrapper functions ---

"""
Wrapper function to execute game api.

e.g.
```
params = Dict("gamePk" => 000000)
gameresult = game(params)
```
"""
function game(params)
    return __executeapi("game", params)
end

"""
Wrapper function to execute game api.

e.g.
```
params = Dict("sportId" => 1, "date" => "2018-07-01")
scheduleresult = schedulemlb(params)
```
"""
function schedulemlb(params)
    return __executeapi("schedule", params)
end

"""
Wrapper function to execute attendance api.

e.g.
```
params = Dict("teamId" => 100)
attendanceresult = attendance(params)
```
"""
function attendance(params)
    return __executeapi("attendance", params)
end

"""
Wrapper function to execute awards api.

e.g.
```
params = Dict("awardId" => "MLBHOF")
hofresult = awards(params)
```
"""
function awards(params)
    return __executeapi("awards", params)
end

"""
Wrapper function to execute draft api.

e.g.
```
params = Dict("year" => 2019)
draftresult = draft(params)
```
"""
function draft(params)
    return __executeapi("draft", params)
end

"""
Wrapper function to execute people api.

e.g.
```
params = Dict("personIds" => 00000)
peopleresult = people(params)
```
"""
function people(params)
    return __executeapi("people", params)
end

"""
Wrapper function to execute standings api.

e.g.
```
params = Dict("leagueId" => 103)
standingsresult = people(params)
```
"""
function standings(params)
    return __executeapi("standings", params)
end

"""
Wrapper function to execute teams api.

e.g.
```
params = Dict()
teamsresult = teams(params)
```
"""
function teams(params)
    return __executeapi("teams", params)
end

"""
Wrapper function to execute meta api.

e.g.
```
params = Dict("type" => "awards")
metaresult = meta(params)
```
"""
function meta(params)
    return __executeapi("meta", params)
end

end


text:
# Make custom function for entry to AMIS so as to overload based on previous distribution
_amis_entry(log_f, previous::Distribution, S, logW, k::Int; nepochs=5, gmm_smps=1000, IS_tilt=1., terminate=0.75) =
    amis(log_f, previous, k; nepochs=nepochs, gmm_smps=gmm_smps, IS_tilt=IS_tilt, terminate=terminate)

_amis_entry(log_f, prev::GMMComp, S, logW, k::Int; nepochs=5, gmm_smps=1000, IS_tilt=1., terminate=0.75) =
    amis(log_f, prev.pis, prev.mus, prev.covs, S, logW; nepochs=nepochs, gmm_smps=gmm_smps, IS_tilt=IS_tilt, terminate=terminate)

"""
     seq_amis(target_logfs, init_dist::Union{Distribution, Int}, k::Int; nepochs=4, IS_tilt=1.3f0, gmm_smps=2000,
              terminate=0.6, max_retry=3, verbose=true)

Perform AMIS on a sequence of target log distributions. The initial distribution is specified with `init_dist`
and this function will perform AMIS using the previously tuned proposal as the initial distribution for the
next target.

Arguments:
- `target_logfs`: an iterable (e.g. vector) where element ``i`` corresponds to the ``i``th log target
    distribution. The distribution need not be normalized, and may often be a log joint.
- `init_dist`: a `Distribution` or `Int`, where the Int specifies the dimension of a standard multivariate
    Gaussian distribution. This is the initial proposal distribution for the first AMIS call.
- `k`: specifies the number of components to use for the proposal GMM.

Optional kwargs.
These are the sam as AMIS, except:
- `max_retry`: AMIS will sometimes fail to provide a good proposal (diagnosed by small ESS of the sample).
    This gives the maximum number of times AMIS will be retried (with increasing tilt and number of epochs).
    See also `min_ess`.
- `min_ess`: The minimum ESS required for each target distribution.
- `verbose`: Whether to display information regarding ESS (and retries) for each target.

Return value:
A tuple containing:
- A vector of the samples and corresponding log weights for each target. Each element can use the utility
    functions: `ess(x)`, `weights(x)` and `resample(x, n)`.
- A vector of the GMM components used for the final proposal for each target.
"""
function seq_amis(target_logfs, init_dist::Union{Distribution, Int}, k::Int; nepochs=4, IS_tilt=1.3f0,
                  gmm_smps=2000, terminate=0.6, max_retry=3, verbose=true, min_ess=min(100, gmm_smps/2))

    t0 = time()
    t_begin = copy(t0)
    if init_dist isa Int
        d = init_dist
        c_dist = MvNormal(zeros(Float32, d), convert(Matrix{Float32}, Matrix(I, d, d)))
    else
        c_dist = init_dist
    end

    smp_out = []
    dist_out = []

    pis, mus, covs = nothing, nothing, nothing    # local scope -- may not need since refactoring.
    S, logW, W = nothing, nothing, nothing

    for (i, lj) in enumerate(target_logfs)

        # reset quantities (in case changed due to non-convergence at prev iter)
        c_IS_tilt, c_nepochs, c_ess = IS_tilt, nepochs, 0

        for j = 1:max_retry
            S, logW, pis, mus, covs = _amis_entry(lj, c_dist, S, W, k; nepochs=c_nepochs, gmm_smps=gmm_smps,
                                                     IS_tilt=c_IS_tilt, terminate=terminate)
            W = softmax(logW)
            c_ess = eff_ss(W)
            (c_ess > 100) && break

            # increase computation / spread of proposal when ESS is low and retry
            # ==> This could be much smarter -- I'm currently throwing away the previous computation.
            # (Although note that not infrequently the previous computation is ~ uninformative.)
            c_IS_tilt, c_nepochs = c_IS_tilt*1.3, c_nepochs + 1
            verbose && printfmtln("Retrying... (ess = {:.1f})", c_ess)
        end
        c_dist = GMMComp(pis, mus, covs)

        push!(smp_out, MCPosterior(S, logW))
        push!(dist_out, GMMComp(pis, mus, covs))
        verbose && printfmtln("Posterior {:02d}, ess={:.2f}, cum. time taken = {:.2f}s ({:.2f})", i, c_ess,
            time()-t_begin, time() - t0); flush(stdout);
        t0 = time()
    end

    return smp_out, dist_out
end

max_stars_repo_name: ornithos/SeqAdaptiveIS

max_issues_repo_name: ornithos/SeqAdaptiveIS

max_forks_repo_name: ornithos/SeqAdaptiveIS


text:
module Normalizers

import StatsBase
using StatsBase: zscore, ZScoreTransform,UnitRangeTransform
using Dates
using DataFrames: DataFrame
using Statistics
import MultivariateStats
using Random

using ..AbsTypes
using ..Utils

import ..AbsTypes: fit, fit!, transform, transform!
export fit, fit!,transform, transform!
export Normalizer

const MV=MultivariateStats

const gmethods = [:zscore,:unitrange,:sqrt,:log,:pca,:ppca,:fa]

"""
    Normalizer(Dict(
       :method => :zscore
    ))

Transforms continuous features into normalized form such as zscore, unitrange, square-root, log, pca, ppca
with parameter:

- `:method` => `:zscore` or `:unitrange` or `:sqrt` or `:log` or `pca` or `ppca` or `fa`
- `:zscore` => standard z-score with centering and scaling
- `:unitrange` => unit range normalization with centering and scaling
- `:sqrt` => square-root transform
- `:pca` => principal component analysis transform
- `:ppca` => probabilistic pca
- `:fa` => factor analysis
- `:log` => log transform

Example:

    function generatedf()
        Random.seed!(123)
        gdate = DateTime(2014,1,1):Dates.Minute(15):DateTime(2016,1,1)
        gval1 = rand(length(gdate))
        gval2 = rand(length(gdate))
        gval3 = rand(length(gdate))
        X = DataFrame(Date=gdate,Value1=gval1,Value2=gval2,Value3=gval3)
        X
    end

    X = generatedf()
    norm = Normalizer(Dict(:method => :zscore))
    fit!(norm,X)
    res=transform!(norm,X)

Implements: `fit!`, `transform!`
"""
mutable struct Normalizer <: Transformer
   name::String
   model::Dict{Symbol,Any}

   function Normalizer(args=Dict())
      default_args = Dict(
         :name => "nrmlzr",
         :method => :zscore
      )
      cargs=nested_dict_merge(default_args,args)
      cargs[:method] in gmethods || throw(ArgmentError("invalid method"))
      cargs[:name] = cargs[:name]*"_"*randstring(3)
      new(cargs[:name],cargs)
   end
end

function Normalizer(st::Symbol)
   Normalizer(Dict(:method=>st))
end

"""
    fit!(st::Statifier, features::T, labels::Vector=[])

Validate argument features other than dates are continuous.
"""
function fit!(norm::Normalizer, features::DataFrame, labels::Vector=[])::Nothing
   # check features are in correct format and no categorical values
   (infer_eltype(features[:,1]) <: DateTime && infer_eltype(Matrix(features[:,2:end])) <: Real) ||
      (infer_eltype(Matrix(features)) <: Real) ||
      throw(ArgmentError("Normalizer.fit!: make sure features are purely float values or float values with Date on first column"))
      return nothing
end

function fit(norm::Normalizer, features::DataFrame, labels::Vector=[])::Normalizer
   fit!(norm,features,labels)
   return deepcopy(norm)
end

"""
    transform!(norm::Normalizer, features::T) where {T<:Union{Vector,Matrix,DataFrame}}

Compute statistics.
"""
function transform!(norm::Normalizer, pfeatures::DataFrame)::DataFrame
   pfeatures != DataFrame() || return DataFrame()
   res = Array{Float64,2}(undef,0,0)
   if (infer_eltype(pfeatures[:,1]) <: DateTime && infer_eltype(Matrix(pfeatures[:,2:end])) <: Real)
      res = processnumeric(norm,Matrix{Float64}(pfeatures[:,2:end])) |> x->DataFrame(x,:auto)
   elseif infer_eltype(Matrix(pfeatures)) <: Real
      features = pfeatures |> Array{Float64}
      res = processnumeric(norm,features) |> x->DataFrame(x,:auto)
   else
      error("Normalizer.transform!: make sure features are purely float values or float values with Date on first column")
   end
   res
end

function transform(norm::Normalizer, pfeatures::DataFrame)::DataFrame
   return transform!(norm,pfeatures)
end

function processnumeric(norm::Normalizer,features::Matrix)
  if norm.model[:method] == :zscore
    ztransform(features)
  elseif norm.model[:method] == :unitrange
    unitrtransform(features)
  elseif norm.model[:method] == :pca
    pca(features)
  elseif norm.model[:method] == :ppca
    ppca(features)
  elseif norm.model[:method] == :ica
    ica(features)
  elseif norm.model[:method] == :fa
    fa(features)
  elseif norm.model[:method] == :sqrt
    sqrtf(features)
  elseif norm.model[:method] == :log
    logf(features)
  else
    error("arg's :method is mapped to unknown keyword")
  end
end

# apply sqrt transform
function sqrtf(X)
   sqrt.(X)
end

# apply log transform
function logf(X)
   log.(X)
end

# apply z-score transform
function ztransform(X)
   xp = X' |> collect |> Matrix{Float64}
   StatsBase.fit(ZScoreTransform, xp,dims=2; center=true, scale=true) |> dt -> StatsBase.transform(dt,xp)' |> collect
end

# unit-range
function unitrtransform(X)
   xp = X' |> collect |> Matrix{Float64}
   StatsBase.fit(UnitRangeTransform,xp,dims=2) |> dt -> StatsBase.transform(dt,xp)' |> collect
end

# pca
function pca(X)
   xp = X' |> collect |> Matrix{Float64}
   m = MV.fit(MV.PCA,xp)
   MV.transform(m,xp)' |> collect
end

function ica(X,kk::Int=0)
   k = kk
   if k == 0
      k = size(X)[2]
   end
   xp = X' |> collect |> Matrix{Float64}
   m = MV.fit(MV.ICA,xp,k)
   MV.transform(m,xp)' |> collect
end

# ppca
function ppca(X)
   xp = X' |> collect |> Matrix{Float64}
   m = MV.fit(MV.PPCA,xp)
   MV.transform(m,xp)' |> collect
end

# fa
function fa(X)
   xp = X' |> collect |> Matrix{Float64}
   m = MV.fit(MV.FactorAnalysis,xp)
   MV.transform(m,xp)' |> collect
end

end


text:
# Data structures for beam search:

# A batch consists of multiple sentences and their beams.
# A Beam consists of a sentence, a beam array and a cand array.
# Both arrays consist of BeamStates sorted by cumulative score.
# We also track the mincost for training and early stop.

# Say we have a n sentence corpus.
# On return we give back n parsers.
# If training we give back x and y.  The size of x and y are uncertain.
# During processing we need a beam for each sentence of size nbeam.
# The beam stores BeamState's.  We allocate BeamState's dynamically.
# Internally we process sentences in batches of size nbatch.
# We need to keep whole of f, no reusing, and copy part of it to return x.
# We could use KUdense or KUsparse for f and use resize... that supports sparse as well.
# that would mean modifying features.jl and sfeatures.jl to use that too.
# calling convention for sparse and dense features.jl are not the same now.
# we could have both append a column at the end of a given array?
# or just have them return the feats vector and we do the appending here using ccat.

# beam: nbeam, beam, sent, mcs, cand, beam2; do we need?
# prealloc parsers but no prealloc of beamstates? parsers dynamic too.

type BeamState
    parent::BeamState   # previous state
    move::Move          # move executed to get to this state
    score::Float64      # cumulative score (logp) TODO: should be unnormalized logp!
    cost::Cost          # cumulative cost (gold arcs that have become impossible)
    fidx::Int           # column in feature matrix representing current parser state # TODO: what is this?
    parser::Parser 	# current parser state
    BeamState()=new()
    BeamState(p::Parser)=new(NullBeamState,0,0,0,0,p)
    BeamState(b::BeamState,m::Integer,s::Number,c::Integer)=
        new(b,convert(Move,m),convert(Float64,s),convert(Cost,c))
end

const NullBeamState = BeamState()
Base.isless(a::BeamState,b::BeamState)=(a.score < b.score)

type Beam
    beam::Vector{BeamState}
    cand::Vector{BeamState}
    sent::Sentence
    stop::Bool                  # TODO: what is this?
    function Beam{T<:Parser}(pt::Type{T}, s::Sentence)
        p = pt(s)
        b = BeamState(p)
        new(BeamState[b], BeamState[], s, !anyvalidmoves(p))
    end
end

typealias Batch Vector{Beam}

function bparse{T<:Parser}(pt::Type{T}, corpus::Corpus, feats::DFvec, net::Net, nbeam::Int;
                           nbatch::Int=1, returnxy::Bool=false, usepmap::Bool=false)
    usepmap ?
    bp_pmap(pt, corpus, feats, net, nbeam, nbatch, returnxy) : # TODO: test pmap
    bp_main(pt, corpus, feats, net, nbeam, nbatch, returnxy)
end

# Life cycle of a Beam:
# 1. bp_resize: beam[] (one element, no fidx), cand[] empty # TODO: what is this? bp_init?
# 2. bp_features: beam[i] gets fidx
# 3. predict: score matrix filled
# 4. bp_update_cand: cand (no parser/fidx) filled with children of beam, sorted
# 5. bp_earlystop: cand[1:nbeam] does not have a 0-cost element stop here
# 6. bp_update_beam: beam filled with top nbeam cand, parsers added
# 7. anyvalidmoves: if beam[1] has no valid moves stop here
# 8. goto step 2

function bp_main{T<:Parser}(pt::Type{T}, corpus::Corpus, feats::DFvec, net::Net, nbeam::Int, nbatch::Int, returnxy::Bool)
    #@date "bp_main Starting"
    if (nbatch < 1 || nbatch > length(corpus)); nbatch = length(corpus); end # TODO: sentences of different length?
    (parses, batch, fmatrix, score) = bp_init(pt, corpus, feats, returnxy) # TODO: replace net with predict::Function
    for s1=1:nbatch:length(corpus)
        s2=min(length(corpus), s1+nbatch-1)
        (batch, fmatrix, score) = bp_init(pt, corpus, s1:s2, nbeam, batch, fmatrix, score) # :86
        nf = 0
        while true              # TODO: how do we train the rnn model? this uses batch for whole corpus!
            nf1 = nf + 1
            nf2 = nf = bp_features(batch, feats, fmatrix, nf) # :35152
            nf2 < nf1 && break                                # TODO: replace sub -> view?
            predict(net, sub(fmatrix,:,nf1:nf2), sub(score,:,nf1:nf2); batch=0) # :26933
            bp_update(batch, score, nbeam, returnxy) # :10750 # TODO: is this where parsing happens?
        end
        # @show nf
        bp_append(parses, batch, fmatrix, score, returnxy) # :1
    end # for s1=1:nbatch:length(corpus)
    #@date "bp_main Finished"   # TODO: what is the input, who computes the features?
    return parses               # TODO: what is the output, what is the loss function? how do we backprop to blstm?
end # function bp_main

function bp_init{T<:Parser}(pt::Type{T}, corpus::Corpus, feats::DFvec, returnxy::Bool)
    # All returned arrays have correct ndims but 0 size
    batch = Array(Beam, 0)
    s0 = corpus[1]
    p0 = pt(s0)
    wt = wtype(corpus)
    fmatrix = Array(wt, (flen(p0,s0,feats), 0))
    score = Array(wt, (p0.nmove, 0))
    # If returnxy we prepare empty x and y matrices
    parses = (if returnxy; (pt[], Matrix{wt}[], Matrix{wt}[]); else; pt[]; end)
    return (parses, batch, fmatrix, score)
end

function bp_init{T<:Parser}(pt::Type{T}, corpus::Corpus, r::UnitRange{Int}, nbeam::Int, batch, fmatrix, score)
    resize!(batch, length(r))
    nword = j = 0
    for i in r
        batch[j+=1] = Beam(pt, corpus[i])
        nword += wcnt(corpus[i])
    end
    fcols = 2*nword*nbeam
    if (fcols > size(fmatrix, 2))
        fcols = round(Int, 1.3*fcols)
        fmatrix = similar(fmatrix, (size(fmatrix,1), fcols))
        score = similar(score, (size(score,1), fcols))
    end
    return (batch, fmatrix, score)
end

function bp_features(batch::Batch, feats::DFvec, fmatrix::Matrix, nf::Int)
    for s in batch
        s.stop && continue
        for bs in s.beam
            size(fmatrix,2) > nf || error("Out of bounds")
            f = features(bs.parser, s.sent, feats, fmatrix, (nf+=1)) # :35029
            bs.fidx = nf
        end
    end
    return nf
end

function bp_update(batch::Batch, score::Matrix, nbeam::Int, returnxy::Bool)
    for s in batch
        s.stop && continue
        bp_update_cand(s, score) # :5128
        returnxy && (s.stop = bp_earlystop(s,nbeam)) && continue
        bp_update_beam(s, nbeam) # :5604
        (s.stop = !anyvalidmoves(s.beam[1].parser)) && continue
    end
end

function bp_update_cand(s::Beam, score::Matrix)
    resize!(s.cand,0)
    cost = Array(Cost, size(score,1))
    for bs in s.beam
        movecosts(bs.parser, s.sent.head, s.sent.deprel, cost) # 162
        for j=1:length(cost)
            cost[j] == typemax(Cost) && continue
            push!(s.cand, BeamState(bs, j, bs.score + score[j,bs.fidx], bs.cost + cost[j])) # :4119
        end
    end
    (length(s.cand) > 0) || error("No candidates found")
    sort!(s.cand; rev=true) # 339
end

function bp_earlystop(s::Beam, nbeam::Int)
    for i=1:min(nbeam,length(s.cand))
        if s.cand[i].cost == 0
            return false
        end
    end
    return true
end

function bp_update_beam(s::Beam, nbeam::Int)
    resize!(s.beam,0)
    for i=1:min(nbeam, length(s.cand))
        bc = s.cand[i]
        bc.parser = copy(bc.parent.parser)     # :5480
        move!(bc.parser, bc.move)              # 20
        push!(s.beam, bc)                      # 3
    end
end

function bp_append(parses, batch::Batch, fmatrix::Matrix, score::Matrix, returnxy::Bool)
    returnxy ?
    bp_append(parses, batch, fmatrix, score) :
    bp_append(parses, batch)
end

function bp_append(parses, batch::Batch)
    for s in batch
        push!(parses, s.beam[1].parser)
    end
end

function bp_append(parses, batch::Batch, fmatrix::Matrix, score::Matrix)
    (p,x,y) = parses
    grad = zeros(Float64,0)
    mark = falses(size(fmatrix,2))
    for s in batch
        push!(p, s.beam[1].parser)
        bp_softmax(s, grad)     # grad <- softmax probabilities from s.cand
        foundgold = false       # in case there is more than one, only treat the highest scoring zero cost answer as gold
        for i=1:length(s.cand)
            bs = s.cand[i]
            bs.cost == 0 && !foundgold && (foundgold=true; grad[i] -= 1) # grad[gold]=p-1, grad[other]=p
            while bs.parent != NullBeamState
                fidx = bs.parent.fidx
                if !mark[fidx]
                    mark[fidx] = true
                    score[:,fidx] = zero(eltype(score))
                end
                score[bs.move,fidx] += grad[i]
                bs = bs.parent
            end
        end
    end
    push!(x, fmatrix[:,mark])
    push!(y, score[:,mark])
end

function bp_softmax(s::Beam, prob::Vector{Float64})
    smax = typemin(Float64)
    z = zero(Float64)
    nc = length(s.cand)
    resize!(prob, nc)
    for bc in s.cand; bc.score > smax && (smax = bc.score); end
    for i=1:nc; prob[i] = exp(s.cand[i].score - smax); z += prob[i]; end
    for i=1:nc; prob[i] /= z; end
end

function print_beam(s::Beam, y::Matrix, f2x::Vector{Int})
    for i = 1:length(s.beam)
        bs = s.cand[i]
        ba = Array(BeamState,0)
        while bs.parent != NullBeamState
            unshift!(ba, bs)
            bs = bs.parent
        end
        for bs in ba
            parent_id = isdefined(bs.parent,:parser) ? hash(bs.parent.parser) & 0xff : 0
            parser_id = isdefined(bs,:parser) ? hash(bs.parser) & 0xff : 0
            grad = y[bs.move, f2x[bs.parent.fidx]]
            @printf "%02x/%d/%.2f/%+.2f/%02x " parent_id bs.cost bs.score grad parser_id
            # print("$(bs.score)/$(bs.cost)/$(bs.grad) ")
            # @printf "%.2f/%d/%d/%d/%+f " bs.score bs.cost bs.fidx bs.parent.fidx y[bs.move, f2x[bs.parent.fidx]]
            # print("$(bs.fidx) ")
            # print("$(pointer(bs.parser.head)) ")
            # bs = bs.parent
        end
        println()
    end
end

function bp_pmap{T<:Parser}(pt::Type{T}, corpus::Corpus, feats::DFvec, net::Net, nbeam::Int, nbatch::Int, returnxy::Bool)
    d = distribute(corpus)
    net = testnet(net)
    #@date "pmap Starting"
    p = pmap(procs(d)) do x
        bp_main(pt, localpart(d), feats, gpucopy(net), nbeam, nbatch, returnxy)
    end
    #@date "pmap Finished"
    return pcat(p)
end

### DEAD CODE

#     batch = [Beam(pt, corpus[i], ndeps) for i in r]
#     s1 = batch[1].sent
#     p1 = batch[1].beam[1].parser
#     nmove = p1.nmove
#     ftype = wtype(s1)
#     frows = flen(p1, s1, feats)
#     nword = 0; for i in r; nword += wcnt(corpus[i]); end
#     fcols = 2 * nword * nbeam
#     fmatrix=Array(ftype, frows, fcols)
#     score=Array(ftype, nmove, fcols)
#     cost=Array(Cost, nmove)
#     return (batch, fmatrix, score, cost)
# end

# function bp_pathsum(s::Beam, fcols::Vector{Int8})
#     ncols = 0
#     for bs in s.beam
#         bs.cost == 0 && (bs.grad -= 1) # gradient is p for bad paths, p-1 for good paths
#         ps = bs.parent
#         while ps != NullBeamState
#             ps.grad += bs.grad
#             ps = ps.parent
#         end
#     end
#     return ncols
# end

        # first resize x/y and init f2x, zero y, copy to x, then do grad sum directly on y
        # dont count, use final shrink, just make sure we have enough space to fit all fmatrix columns

        # ncols = bp_pathsum(s, fcols)
        # print_beam(s); println()
        # we need to put bs.parent.fidx into x, y[bs.move]+=bs.grad, multiple children can update same y

    # if returnxy
    #     @show map(typeof, parses)
    #     @show map(size, parses)
    #     parses = (parses[1], parses[2].arr, parses[3].arr)
    # end

# Return values:
# single sentence: nbatch=1, corpus[1]
# bp_append is called for every batch.
# so nbatch sentences: nbatch parses, one x one y matrix?
# don't grow xy, gc fmatrix/score, copy?
# return parses followed by a bunch of x/y pairs?
# parse array, x array, y array?  grow using push?  or just merge?
# or just the parse array.
# problem with just merge is pmap does not have an easy merge.
# so we will go with p,x,y vectors grown with push.

#     return (returnxy ? (parses[1], copy(parses[2].arr), copy(parses[3].arr)) : parses)
#     (parses, batch, fmatrix, score) = bp_init(pt, corpus, ndeps, feats, returnxy)
# (returnxy ? (Any[], Any[], Any[]) : Any[])

# function bp_append(parses, batch::Batch, fmatrix::KUdense, returnxy::Bool) # returnxy version
#     returnxy || (return bp_append(parses,batch))
#     (p, x, y) = parses
#     (mx,nx) = size(x)
#     nf = size(fmatrix,2)
#     resize!(x, (size(x,1), nx+nf))
#     resize!(y, (size(y,1), nx+nf))
#     f2x = zeros(Int, nf) # f2x[fidx] => xidx
#     grad = zeros(Float64,0)
#     for s in batch
#         push!(p, s.beam[1].parser)
#         bp_softmax(s, grad)
#         foundgold = false       # in case there is more than one, only treat the highest scoring zero cost answer as gold
#         for i=1:length(s.cand)
#             bs = s.cand[i]
#             bs.cost == 0 && !foundgold && (foundgold=true; grad[i] -= 1)
#             while bs.parent != NullBeamState
#                 fidx = bs.parent.fidx
#                 xidx = f2x[fidx]
#                 if xidx == 0
#                     f2x[fidx] = xidx = (nx += 1)
#                     nx <= size(x,2) || error("Out of bounds")
#                     copy!(x, (xidx-1)*mx+1, fmatrix, (fidx-1)*mx+1, mx)
#                     y[:,xidx] = zero(eltype(y))
#                 end
#                 y[bs.move,xidx] += grad[i]
#                 bs = bs.parent
#             end
#         end
#         # s.sent.form[1]=="Influential" && print_beam(s, y.arr, f2x)
#     end
#     # @show nx
#     resize!(x, (size(x,1), nx))
#     resize!(y, (size(y,1), nx))
# end


text:
function parser_json_unit_tests()
    parse_from_json_tests()
end

function parse_from_json_tests()
    data = RVRP.parse_from_json_string("{\"id\":\"wololo\"}\n")
    @test typeof(data) == RVRP.RvrpInstance
    @test data.locations == RVRP.Location[]
    @test data.location_groups == RVRP.LocationGroup[]
    @test data.product_compatibility_classes == RVRP.ProductCompatibilityClass[]
    @test data.product_sharing_classes == RVRP.ProductSharingClass[]
    @test data.product_specification_classes == RVRP.ProductSpecificationClass[]
    @test data.requests == RVRP.Request[]
    @test data.vehicle_categories == RVRP.VehicleCategory[]
    @test data.vehicle_sets == RVRP.HomogeneousVehicleSet[]

    # data = RVRP.parse_cvrplib(dirname(@__FILE__)* "/../../data/cvrplib/P/P\\P-n16-k8.vrp")
    # @test isfile(dirname(@__FILE__)*"/../../dump/dummy.json") == false
    # RVRP.parse_to_json(data, dirname(@__FILE__)*"/../../dump/dummy.json")
    # @test isfile(dirname(@__FILE__)*"/../../dump/dummy.json") == true
    # data2 = RVRP.parse_from_json(dirname(@__FILE__)*"/../../dump/dummy.json")
    # @test data == data2
    # rm(dirname(@__FILE__)*"/../../dump/dummy.json")
    # @test isfile(dirname(@__FILE__)*"/../../dump/dummy.json") == false
end


text:
using Test
using Smoothers
import Smoothers: filter

@testset "filter" begin

    # Moving Average Filter
    w = 3;
    b = ones(w)/w;
    a = [1];
    x = [1,2,3,4,5];

    # Types' promotion
    fx = filter(Float16.(b),Float16.(a),Float16.(x));
    @test eltype(fx) == Float16
    fx = filter(Float32.(b),Float32.(a),Float32.(x));
    @test eltype(fx) == Float32
    fx = filter(BigFloat.(b),BigFloat.(a),BigFloat.(x));
    @test eltype(fx) == BigFloat
    fx = filter(Float16.(b),Int32.(a),Float64.(x))
    @test eltype(fx) == Float64

    # Standard
    fx = filter(b,a,x)
    @test eltype(fx) == Float64
    @test fx[1:2] ≈ [1/3,1]
    @test fx[3:5] ≈ sma(x,3)

    # Initial State
    si = [1/3,1]
    fx = filter(b,a,x,si)
    @test fx[1:2] ≈ [2/3,2]
    @test fx[3:5] ≈ sma(x,3)

end


text:
@testset "varinfo.jl" begin
    @testset "TypedVarInfo" begin
        @model gdemo(x, y) = begin
            s ~ InverseGamma(2, 3)
            m ~ truncated(Normal(0.0, sqrt(s)), 0.0, 2.0)
            x ~ Normal(m, sqrt(s))
            y ~ Normal(m, sqrt(s))
        end
        model = gdemo(1.0, 2.0)

        vi = VarInfo()
        model(vi, SampleFromUniform())
        tvi = TypedVarInfo(vi)

        meta = vi.metadata
        for f in fieldnames(typeof(tvi.metadata))
            fmeta = getfield(tvi.metadata, f)
            for vn in fmeta.vns
                @test tvi[vn] == vi[vn]
                ind = meta.idcs[vn]
                tind = fmeta.idcs[vn]
                @test meta.dists[ind] == fmeta.dists[tind]
                @test meta.orders[ind] == fmeta.orders[tind]
                @test meta.gids[ind] == fmeta.gids[tind]
                for flag in keys(meta.flags)
                    @test meta.flags[flag][ind] == fmeta.flags[flag][tind]
                end
                range = meta.ranges[ind]
                trange = fmeta.ranges[tind]
                @test all(meta.vals[range] .== fmeta.vals[trange])
            end
        end
    end
    @testset "Base" begin
        # Test Base functions:
        #   string, Symbol, ==, hash, in, keys, haskey, isempty, push!, empty!,
        #   getindex, setindex!, getproperty, setproperty!
        csym = gensym()
        vn1 = @varname x[1][2]
        @test string(vn1) == "x[1][2]"
        @test Symbol(vn1) == Symbol("x[1][2]")

        vn2 = @varname x[1][2]
        @test vn2 == vn1
        @test hash(vn2) == hash(vn1)
        @test inspace(vn1, (:x,))

        function test_base!(vi)
            empty!(vi)
            @test getlogp(vi) == 0
            @test get_num_produce(vi) == 0

            vn = @varname x
            dist = Normal(0, 1)
            r = rand(dist)
            gid = Selector()

            @test isempty(vi)
            @test ~haskey(vi, vn)
            @test !(vn in keys(vi))
            push!(vi, vn, r, dist, gid)
            @test ~isempty(vi)
            @test haskey(vi, vn)
            @test vn in keys(vi)

            @test length(vi[vn]) == 1
            @test length(vi[SampleFromPrior()]) == 1

            @test vi[vn] == r
            @test vi[SampleFromPrior()][1] == r
            vi[vn] = [2 * r]
            @test vi[vn] == 2 * r
            @test vi[SampleFromPrior()][1] == 2 * r
            vi[SampleFromPrior()] = [3 * r]
            @test vi[vn] == 3 * r
            @test vi[SampleFromPrior()][1] == 3 * r

            empty!(vi)
            @test isempty(vi)
            push!(vi, vn, r, dist, gid)

            function test_inspace()
                space = (:x, :y, @varname(z[1]), @varname(M[1:10, :]))

                @test inspace(@varname(x), space)
                @test inspace(@varname(y), space)
                @test inspace(@varname(x[1]), space)
                @test inspace(@varname(z[1][1]), space)
                @test inspace(@varname(z[1][:]), space)
                @test inspace(@varname(z[1][2:3:10]), space)
                @test inspace(@varname(M[[2, 3], 1]), space)
                @test inspace(@varname(M[:, 1:4]), space)
                @test inspace(@varname(M[1, [2, 4, 6]]), space)
                @test !inspace(@varname(z[2]), space)
                @test !inspace(@varname(z), space)
            end
            return test_inspace()
        end
        vi = VarInfo()
        test_base!(vi)
        test_base!(empty!(TypedVarInfo(vi)))
    end
    @testset "flags" begin
        # Test flag setting:
        #    is_flagged, set_flag!, unset_flag!
        function test_varinfo!(vi)
            vn_x = @varname x
            dist = Normal(0, 1)
            r = rand(dist)
            gid = Selector()

            push!(vi, vn_x, r, dist, gid)

            # del is set by default
            @test !is_flagged(vi, vn_x, "del")

            set_flag!(vi, vn_x, "del")
            @test is_flagged(vi, vn_x, "del")

            unset_flag!(vi, vn_x, "del")
            @test !is_flagged(vi, vn_x, "del")
        end
        vi = VarInfo()
        test_varinfo!(vi)
        test_varinfo!(empty!(TypedVarInfo(vi)))
    end
    @testset "setgid!" begin
        vi = VarInfo()
        meta = vi.metadata
        vn = @varname x
        dist = Normal(0, 1)
        r = rand(dist)
        gid1 = Selector()
        gid2 = Selector(2, :HMC)

        push!(vi, vn, r, dist, gid1)
        @test meta.gids[meta.idcs[vn]] == Set([gid1])
        setgid!(vi, gid2, vn)
        @test meta.gids[meta.idcs[vn]] == Set([gid1, gid2])

        vi = empty!(TypedVarInfo(vi))
        meta = vi.metadata
        push!(vi, vn, r, dist, gid1)
        @test meta.x.gids[meta.x.idcs[vn]] == Set([gid1])
        setgid!(vi, gid2, vn)
        @test meta.x.gids[meta.x.idcs[vn]] == Set([gid1, gid2])
    end
    @testset "setval! & setval_and_resample!" begin
        @model function testmodel(x)
            n = length(x)
            s ~ truncated(Normal(), 0, Inf)
            m ~ MvNormal(zeros(n), I)
            return x ~ MvNormal(m, s^2 * I)
        end

        @model function testmodel_univariate(x, ::Type{TV}=Vector{Float64}) where {TV}
            n = length(x)
            s ~ truncated(Normal(), 0, Inf)

            m = TV(undef, n)
            for i in eachindex(m)
                m[i] ~ Normal()
            end

            for i in eachindex(x)
                x[i] ~ Normal(m[i], s)
            end
        end

        x = randn(5)
        model_mv = testmodel(x)
        model_uv = testmodel_univariate(x)

        for model in [model_uv, model_mv]
            m_vns = model == model_uv ? [@varname(m[i]) for i in 1:5] : @varname(m)
            s_vns = @varname(s)

            vi_typed = VarInfo(model)
            vi_untyped = VarInfo()
            model(vi_untyped, SampleFromPrior())

            for vi in [vi_untyped, vi_typed]
                vicopy = deepcopy(vi)

                ### `setval` ###
                DynamicPPL.setval!(vicopy, (m=zeros(5),))
                # Setting `m` fails for univariate due to limitations of `setval!`
                # and `setval_and_resample!`. See docstring of `setval!` for more info.
                if model == model_uv
                    @test_broken vicopy[m_vns] == zeros(5)
                else
                    @test vicopy[m_vns] == zeros(5)
                end
                @test vicopy[s_vns] == vi[s_vns]

                # Ordering is NOT preserved => fails for multivariate model.
                DynamicPPL.setval!(
                    vicopy, (; (Symbol("m[$i]") => i for i in (1, 3, 5, 4, 2))...)
                )
                if model == model_uv
                    @test vicopy[m_vns] == 1:5
                else
                    @test vicopy[m_vns] == [1, 3, 5, 4, 2]
                end
                @test vicopy[s_vns] == vi[s_vns]

                DynamicPPL.setval!(
                    vicopy, (; (Symbol("m[$i]") => i for i in (1, 2, 3, 4, 5))...)
                )
                DynamicPPL.setval!(vicopy, (s=42,))
                @test vicopy[m_vns] == 1:5
                @test vicopy[s_vns] == 42

                ### `setval_and_resample!` ###
                if model == model_mv && vi == vi_untyped
                    # Trying to re-run model with `MvNormal` on `vi_untyped` will call
                    # `MvNormal(μ::Vector{Real}, Σ)` which causes `StackOverflowError`
                    # so we skip this particular case.
                    continue
                end

                vicopy = deepcopy(vi)
                DynamicPPL.setval_and_resample!(vicopy, (m=zeros(5),))
                model(vicopy)
                # Setting `m` fails for univariate due to limitations of `subsumes(::String, ::String)`
                if model == model_uv
                    @test_broken vicopy[m_vns] == zeros(5)
                else
                    @test vicopy[m_vns] == zeros(5)
                end
                @test vicopy[s_vns] != vi[s_vns]

                # Ordering is NOT preserved.
                DynamicPPL.setval_and_resample!(
                    vicopy, (; (Symbol("m[$i]") => i for i in (1, 3, 5, 4, 2))...)
                )
                model(vicopy)
                if model == model_uv
                    @test vicopy[m_vns] == 1:5
                else
                    @test vicopy[m_vns] == [1, 3, 5, 4, 2]
                end
                @test vicopy[s_vns] != vi[s_vns]

                # Correct ordering.
                DynamicPPL.setval_and_resample!(
                    vicopy, (; (Symbol("m[$i]") => i for i in (1, 2, 3, 4, 5))...)
                )
                model(vicopy)
                @test vicopy[m_vns] == 1:5
                @test vicopy[s_vns] != vi[s_vns]

                DynamicPPL.setval_and_resample!(vicopy, (s=42,))
                model(vicopy)
                @test vicopy[m_vns] != 1:5
                @test vicopy[s_vns] == 42
            end
        end

        # https://github.com/TuringLang/DynamicPPL.jl/issues/250
        @model function demo()
            return x ~ filldist(MvNormal([1, 100], I), 2)
        end

        vi = VarInfo(demo())
        vals_prev = vi.metadata.x.vals
        ks = [@varname(x[1, 1]), @varname(x[2, 1]), @varname(x[1, 2]), @varname(x[2, 2])]
        DynamicPPL.setval!(vi, vi.metadata.x.vals, ks)
        @test vals_prev == vi.metadata.x.vals

        DynamicPPL.setval_and_resample!(vi, vi.metadata.x.vals, ks)
        @test vals_prev == vi.metadata.x.vals
    end
end


text:
"""

    verify_jld2_data(path::AbstractString; extract_recent::Bool = false)

Verify if JLD2 files exist and match ZIP archives in `path`. Only the archives are version-controlled
and hosted on the remote repository.
If a JLD2 does not exist, it will be extracted from the archives and written to disk.

A warning will be returned if ZIP archives are more recent than their corresponding JLD2 files,
likely indicating upstream changes. Re-run with `extract_recent` set to `true` to extract the newer files.
"""
function verify_jld2_data(path::AbstractString; extract_recent::Bool=false)
    # List zip files
    files = readdir(path; join=true)
    zipfiles = filter(x -> occursin(".zip", x), files)
    # Check if corresponding .jld2 exist
    jldfiles = replace.(zipfiles, ".zip" => ".jld2")
    missing_files = map(!isfile, jldfiles)
    # Extract missing files from archive & write to disk
    if any(missing_files)
        @warn """\n
              One (or more) JLD2 file is missing.
              It will be extracted from the ZIP archives and written to disk.
              """
        for (i, j) in zip(findall(missing_files), 1:sum(missing_files))
            @info "Reading $(jldfiles[i]) from archive ($(j)/$(sum(missing_files)))"
            _unzip_jld2(zipfiles[i], jldfiles[i])
        end
    end
    # Warn if a ZIP archive is more recent than its corresponding JLD2 file
    more_recent = mtime.(zipfiles) .> mtime.(jldfiles)
    if any(more_recent)
        @warn """
              \n One (or more) ZIP archive is more recent than its corresponding JLD2 file.
              \n Run `verify_jld2_data("$(path)"; extract_recent = true)` to extract newer version.
              """
        # Extract JLD2 file from more recent archive, if specified
        if extract_recent
            for (i, j) in zip(findall(more_recent), 1:sum(more_recent))
                @info "Reading $(jldfiles[i]) from more recent archive ($(j)/$(sum(more_recent)))"
                _unzip_jld2(zipfiles[i], jldfiles[i])
            end
        end
    end
end
# @time verify_jld2_data(joinpath("data", "jld2"))
# @time verify_jld2_data(joinpath("data", "jld2"); extract_recent = true)

function _unzip_jld2(zip_path, jld_path)
    zf = ZipFile.Reader(zip_path)
    open(jld_path, "w") do io
        write(io, read(zf.files[1]))
    end
end

function _zip_jld2(zip_path, jld_path)
    w = ZipFile.Writer(zip_path) # create archive folder
    filename = split(jld_path, "/")[end]
    f = ZipFile.addfile(w, filename; method=ZipFile.Deflate) # create file for compression
    open(jld_path, "r") do io # file to compress
        write(f, io) # compress file, ~60sec
    end
    return close(w)
end
# @time _zip_jld2(zipfiles[2], jldfiles[2])

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
@for_libpetsc begin

    function incref(::Type{$PetscScalar}, obj)
        @chk ccall((:PetscObjectReference, $libpetsc), PetscErrorCode, (Ptr{Cvoid},), obj)
        return nothing
    end

    function decref(::Type{$PetscScalar}, obj)
        if !finalized($PetscScalar)
            @chk ccall((:PetscObjectDereference, $libpetsc), PetscErrorCode, (Ptr{Cvoid},), obj)
            obj.ptr = C_NULL
        end
        return nothing
    end

    function nrefs(::Type{$PetscScalar}, obj)
        r_n = Ref{$PetscInt}()
        @chk ccall((:PetscObjectGetReference, $libpetsc), PetscErrorCode, (Ptr{Cvoid},Ptr{$PetscInt}), obj, r_n)
        return r_n[]
    end

end

"""
    incref(obj)

Increment the reference counter fo `obj`.
This usually only needs to be called when accessing objects owned by other objects, e.g. via `KSPGetPC`.
"""
incref(obj) = incref(scalartype(obj), obj)

"""
    decref(obj)

Decrement the reference counter for `obj`.

In general we don't need to use this, as we can call `destroy` instead.
"""
decref(obj) = decref(scalartype(obj), obj)

"""
    nrefs(obj)

The current reference count for `obj`.
"""
nrefs(obj)  = nrefs(scalartype(obj), obj)


text:
# This example is used to generate plots in the paper, "Reflection from a multi-species material and its transmitted effective wavenumber." Proc. R. Soc. (2018): 20170864.

# Everything related to ploting has been commented so that this example can run with requiring the Plots package. That way, this example can be checked everytime the package is modified.

using EffectiveWaves

# using LaTeXStrings
# using Plots
# height=500
# unicodeplots()
 # pyplot(linewidth=3, size=(2.6*height,height), border=false)

 # Plots.scalefontsizes(1.8)

filename="fluid"
mediumname = "water"
## choose material

  # fluid
  medium = WaterDistilled
  # inclusion1 = SodiumNitrate # surfactant ?
  inclusion1 = Glycerol # surfactant ?
  inclusion2 = Hexadecane

  # medium = Hexadecane
  # # inclusion1 = SodiumNitrate # surfactant ?
  # inclusion1 = SodiumNitrate # surfactant ?
  # inclusion2 = Glycerol

  ωfactor = 4000;
  ωs = 2*ωfactor.*LinRange(real(medium.c/10000),real(medium.c),100) # k from 0 to 1
  # ωs = 2.0*pi*LinRange(1.0e1,1.0e7,100)

  volfrac = 0.22
  r1 = 0.1/ωfactor; vol1 = 0.11;
  # r1 = 0.1/ωfactor; vol1 = 0.06;
  r2 = 1.0/ωfactor; vol2 = volfrac - vol1

  sp1 = Specie(ρ=inclusion1.ρ, r=r1, c=inclusion1.c, volfrac = vol1)
  sp2 = Specie(ρ=inclusion2.ρ, r=r2, c=inclusion2.c, volfrac = volfrac-vol1)

  # True wavenumber
  kTs  = wavenumber_low_volumefraction(ωs, medium, [sp1,sp2]);
  # Approximate wavenumber
  kTLSs = two_species_approx_wavenumber(ωs, medium, [sp1,sp2]);
  # Approximate Challis wavenumber
  kTCs = wavenumber_challis(ωs, medium, [sp1,sp2]);

  speed_arr = [ ωs./real(kTs), ωs./real(kTLSs), ωs./real(kTCs), 0 .*ωs .+ real(medium.c)]
  atten_arr = imag([kTs,kTLSs,kTCs])

  styles = [:solid :dashdot :dashdot :dot]
  # labs = [L"Effective $k_{*}$" L"Approximate $k_{*LS}$" L"Approximate $k_{*C}$" mediumname];
  # xlabel = L"k a_S";
  xs = r1.*ωs ./ real(medium.c);
  ys_arr = speed_arr;
  m =5;
  # p1= plot(xs, ys_arr, xlabel=xlabel, ylabel="sound speed (m/s)", labels=labs, line = styles
  #           , ylims = (minimum(minimum.(ys_arr))*0.995, maximum(maximum.(ys_arr))*1.005));

  ys_arr = atten_arr;
  # labs = [L"Effective $k_{*}$" L"Approximate $k_{*LS}$" L"Approximate $k_{*C}$"];
  # p2 = plot(xs, ys_arr, labels=labs, xlabel=xlabel, ylabel="attenuation (1/m)", line = styles
  #             , ylims = (minimum(minimum.(ys_arr))*0.995, maximum(maximum.(ys_arr))*1.005));
  # plot(p1,p2)
  # try mkdir("media") end
  # savefig("media/compare_$(filename).png")
  # savefig("media/compare_$(filename).pdf")
  #
  # gui()

## How small does inclusion1 need to be ?
  ω = ωfactor * real(medium.c)/4.0
  r1s = (0.002:0.01:0.5) ./ ωfactor;

  sp2 = Specie(ρ=inclusion2.ρ ,r=r2, c=inclusion2.c, volfrac = volfrac-vol1)

  # True wavenumber
  kTs = map(r1s) do r1
    sp1 = Specie(ρ=inclusion1.ρ ,r=r1, c=inclusion1.c, volfrac = vol1)
    wavenumber_low_volumefraction(ω, medium, [sp1,sp2])
  end
  # Approximate wavenumber
  kTLSs = map(r1s) do r1
    sp1 = Specie(ρ=inclusion1.ρ ,r=r1, c=inclusion1.c, volfrac = vol1)
    two_species_approx_wavenumber(ω, medium, [sp1,sp2])
  end
  # Approximate Challis wavenumber
  kTCs = map(r1s) do r1
    sp1 = Specie(ρ=inclusion1.ρ ,r=r1, c=inclusion1.c, volfrac = vol1)
    wavenumber_challis(ω, medium, [sp1,sp2])
  end

  speed_arr = [ ω./real(kTs), ω./real(kTLSs), ω./real(kTCs), 0 .*real(kTLSs) .+ real(medium.c)]
  atten_arr = imag([kTs,kTLSs,kTCs])

  styles = [:solid :dashdot :dashdot :dot]
  # labs = [L"Effective $k_{*}$" L"Approximate $k_{*LS}$" L"Approximate $k_{*C}$" mediumname]
  ys_arr = speed_arr;
  xs = r1s.*ω ./ real(medium.c);
  # xlabel=L"k_0 a_S"
  # m =5;
  # p1 = plot(xs, ys_arr, xlabel=xlabel, ylabel="sound speed (m/s)", labels=labs , line = styles
            #   , ylims = (minimum(minimum.(ys_arr))*0.995, maximum(maximum.(ys_arr))*1.005));
  # y1 = min(ys_arr[1][1:m]..., ys_arr[2][1:m]...);
  # y2 = max(ys_arr[1][1:m]..., ys_arr[2][1:m]...);
  # p1 = gray_square!([xs[1],xs[m]],[y1,y2]);

  ys_arr = atten_arr;
  # labs = [L"Effective $k_{*}$" L"Approximate $k_{*LS}$" L"Approximate $k_{*C}$"]
  # p2 = plot(xs, ys_arr, labels=labs, xlabel=xlabel, ylabel="attenuation (1/m)", line=styles
            #   , ylims = (minimum(minimum.(ys_arr))*0.995, maximum(maximum.(ys_arr))*1.005));
  # plot(p1,p2)
  # gui()
  # savefig("media/compare_$(filename)_small.png")
  # savefig("media/compare_$(filename)_small.pdf")

# Plots.scalefontsizes(1/1.8)

max_stars_repo_name: arturgower/EffectiveWaves


text:
module ANBIMA

using HTTP
using Dates
using BusinessDays
using Evolutionary

include("types.jl")
include("IMA.jl")
include("IDA.jl")
include("ettj.jl")
include("credit_spread.jl")

function parse_value(T::Type, x::String)
    if x == "--" ||  x == ""
        return nothing
    end
    return parse(T, replace(replace(x, r"\." => ""), r"," => "."))
end

parse_value(T::Type, x::SubString{String}) = parse_value(T, String(x))

end # module


text:
const PopOrInput = Union{Population, InputPopulation}

struct NetworkEdge{WT<:AbstractMatrix{<:Real}, ST<:AbstractArray{<:AbstractSynapse, 2}, LT<:AbstractLearner}
    weights::WT
    synapses::ST
    learner::LT
end

struct Network <: AbstractDict{Symbol, PopOrInput}
    pops::Dict{Symbol, PopOrInput}
    fedgelist::Dict{Symbol, Vector{Symbol}}
    bedgelist::Dict{Symbol, Vector{Symbol}}
    connections::Dict{Tuple{Symbol, Symbol}, NetworkEdge}
end

"""
    size(net::Network)

Return the number of populations in a network.
"""
Base.length(net::Network) = length(net.pops)
Base.size(net::Network) = length(net)

Base.iterate(net::Network) = iterate(net.pops)
Base.getindex(net::Network, key) = net.pops[key]
Base.setindex!(net::Network, pop, key) = (net.pops[key] = pop)

Base.show(io::IO, net::Network) = print(io, "Network($(size(net)))")
function Base.show(io::IO, ::MIME"text/plain", net::Network)
    println(io, "Network($(size(net))):")
    for (name, pop) in net.pops
        println(io, "  $name => $pop")
    end
end

Network(pops::Dict = Dict{Symbol, PopOrInput}()) =
    Network(pops, Dict{Symbol, Vector{Symbol}}(),
                  Dict{Symbol, Vector{Symbol}}(),
                  Dict{Tuple{Symbol, Symbol}, NetworkEdge}())

function connect!(net::Network, src::Symbol, dst::Symbol; weights::AbstractMatrix{<:Real}, synapse = Synapse.Delta, learner = George())
    !(haskey(net.pops, src) && haskey(net.pops, dst)) && error("Cannot find populations called $src and/or $dst in network.")
    fedges = get!(net.fedgelist, src, Symbol[])
    push!(fedges, dst)
    bedges = get!(net.bedgelist, dst, Symbol[])
    push!(bedges, src)

    m = size(net.pops[src])
    n = size(net.pops[dst])
    synapses = StructArray(synapse() for i in 1:m, j in 1:n; unwrap = t -> t <: AbstractSynapse)

    net.connections[(src, dst)] = NetworkEdge(weights, synapses, learner)
end

function _process_spikes!(net::Network, t::Integer, spikes; dt::Real = 1.0)
    for (pop, spikevec) in spikes
        dsts = get!(net.fedgelist, pop, Symbol[])
        for dst in dsts
            edge = net.connections[(pop, dst)]
            weights = edge.weights
            synapses = edge.synapses

            # excite synapses
            map((row, s) -> (s > 0) && excite!(row, s + 1), eachrow(synapses), spikevec)

            # compute current
            current = vec(reduce(+, weights .* evaluate!(synapses, t + 1; dt = dt); dims = 1))
            excite!(net.pops[dst].somas, current)
        end

        srcs = get!(net.bedgelist, pop, Symbol[])
        for src in srcs
            # apply refactory period to synapses
            edge = net.connections[(src, pop)]
            map((col, s) -> (s > 0) && spike!(col, s; dt = dt), eachcol(edge.synapses), spikevec)
        end
    end
end

function _record_spikes!(net::Network, spikes; dt::Real = 1.0)
    for (pop, spikevec) in spikes
        dsts = get!(net.fedgelist, pop, Symbol[])
        for dst in dsts
            # record pre-synaptic spikes
            edge = net.connections[(pop, dst)]
            prespike!(edge.learner, edge.weights, spikevec; dt = dt)
        end

        srcs = get!(net.bedgelist, pop, Symbol[])
        for src in srcs
            # record post-synaptic spikes
            edge = net.connections[(src, pop)]
            postspike!(edge.learner, edge.weights, spikevec; dt = dt)
        end
    end
end

_evalnode(node::Population, t; dt, dense) = node(t; dt = dt, dense = dense)
_evalnode(node::InputPopulation, t; dt, dense) = node(t; dt = dt)

_resetnode!(node::Population) = reset!(node)
_resetnode!(node::InputPopulation) = nothing

function update!(net::Network, t::Integer, spikes; dt::Real = 1.0)
    @inbounds for (name, pop) in net.pops
        (pop isa Population) && update!(pop, t, spikes[name]; dt = dt)
    end

    _record_spikes!(net, spikes; dt = dt)

    @inbounds for (_, edge) in net.connections
        update!(edge.learner, edge.weights, t; dt = dt)
    end
end

function (net::Network)(t::Integer; dt::Real = 1.0, dense = false)
    spikes = Dict{Symbol, Union{Vector, CuVector}}()

    @inbounds for (name, pop) in net.pops
        spikes[name] = _evalnode(pop, t; dt = dt, dense = dense)
    end

    _process_spikes!(net, t, spikes; dt = dt)

    return spikes
end

function reset!(net::Network)
    _resetnode!.(values(net.pops))
    for (_, edge) in net.connections
        reset!(edge.synapses)
    end
end

function simulate!(net::Network, T::Integer; dt::Real = 1.0, cb = () -> (), dense = false)
    spiketimes = Dict{Symbol, Dict}()

    for t = 1:T
        cb()

        spikes = net(t; dt = dt, dense = dense)

        for (name, spikevec) in spikes
             # record spikes
            dict = get!(spiketimes, name, Dict([(i, Int[]) for i in 1:size(net.pops[name])]))
            _recordspikes!(dict, spikevec)
        end

        update!(net, t, spikes; dt = dt)
    end

    return spiketimes
end


text:
using Tullio, Test, ForwardDiff
# using Tracker; _gradient(x...) = Tracker.gradient(x...); GRAD = :Tracker

function gradtest(f, dims)
    x = randn(dims...)
    grad = ForwardDiff.gradient(x -> sum(sin, f(x)), x)
    grad ≈ _gradient(x -> sum(sin, f(x)), x)[1]
end

@testset "from TensorTrace" begin
    # These can all be handled using TensorOperations

    triv1(x) = Tullio.@tensor A[i,j] := 2 * x[i,j]
    @test gradtest(triv1, (2,3))

    r32 = randn(3,2);
    r312 = randn(3,1,2);

    ## trace!
    tr1(x) = Tullio.@tensor T[k] := 22 * x[i,i,k]
    @test gradtest(tr1, (3,3,4))

    tr2(x) = Tullio.@tensor T[k] := 22 * x[i,i,k,j,j]
    @test gradtest(tr2, (3,3,4,7,7))

    ## contract! A
    con1(x) = Tullio.@tensor C[i,j] := 5 * x[i,k] * r32[k,j]
    @test gradtest(con1, (2,3))

    r22 = rand(2,2);

    con3(x) = Tullio.@tensor C[i,j,m,n] := x[i,j,k] * r312[k,m,n]
    @test gradtest(con3, (1,2,3))

    con4(x) = Tullio.@tensor C[i,m] := x[i,kk,k] * r312[k,m,kk]
    @test gradtest(con4, (1,2,3))

    con5(x) = Tullio.@tensor C[j,i,n,m] := 44 * x[i,j,k] * r312[k,m,n]
    @test gradtest(con5, (1,2,3))

    r392 = randn(3,9,2);
    con6(x) = Tullio.@tensor C[n,i,m,j] := x[i,j,k] * r392[k,m,n]
    @test gradtest(con6, (9,2,3))

    con7(x) = Tullio.@tensor C[m,n,j,i] := 44 * x[i,j,k] * r392[k,m,n]
    @test gradtest(con7, (9,2,3))

    ## contract! B
    con8b(x) = Tullio.@tensor K[i,j] := 5 * r32[i,k] * x[k,j]
    @test gradtest(con8b, (2,3))

    con9b(x) = Tullio.@tensor K[i,j,m,n] := r312[i,j,k] * x[m,k,n]
    @test gradtest(con9b, (1,2,3))

    con10b(x) = Tullio.@tensor K[n,j,m,i] := r392[i,j,k] * x[m,k,n]
    @test gradtest(con10b, (9,2,3))

    r3399 = randn(3,3,9,9);

    con13(x) = Tullio.@tensor K[i,j] := r3399[s,s,j,k] * x[t,t,k,i]
    @test gradtest(con13, (3,3,9,9))

    r33 = rand(3,3);
    con14(x) = Tullio.@tensor K[i,j] := r3399[a,b,j,k] * x[b,c,k,i] * r33[a,c]
    @test gradtest(con14, (3,3,9,9))

    ## scalar -- one with :=, one without
    sc1(x) = Tullio.@tensor s = r22[b,β] * x[a,b,c] * r312[c,a,β]
    @test gradtest(sc1, (1,2,3))

    sc2(x) = Tullio.@tensor s := x[γ,c] * r3399[c,γ,i,i]
    @test gradtest(sc2, (3,3))

end

@testset "errors" begin
    @test_throws LoadError @eval Tullio.@tensor C[k] := A[i,i,k] + B[k]  # two terms
    @test_throws LoadError @eval Tullio.@tensor B[k] := conj(A[k])  # functions
    @test_throws LoadError @eval Tullio.@tensor C[k] := A[i, i+k]  # not a contraction
end


text:
using Documenter
using Tensorial

# Setup for doctests in docstrings
DocMeta.setdocmeta!(Tensorial, :DocTestSetup, recursive = true,
    quote
        using Tensorial
        using Random
        Random.seed!(1234)
    end
)

makedocs(;
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    modules = [Tensorial],
    sitename = "Tensorial.jl",
    pages=[
        "Home" => "index.md",
        "Getting Started" => [
            "Cheat Sheet.md",
        ],
        "Manual" => [
            "Tensor Type.md",
            # "Constructors.md",
            "Tensor Operations.md",
            "Continuum Mechanics.md",
            "Voigt form.md",
            "Tensor Inversion.md",
            "Broadcast.md",
            "Einstein summation.md",
            "Quaternion.md",
        ],
        "Benchmarks.md",
    ],
    doctest = true, # :fix
)

deploydocs(
    repo = "github.com/KeitaNakamura/Tensorial.jl.git",
    devbranch = "main",
)


text:
A = [2 3 0;
     1 2 4;
     0 2 1]

f_m(X) = X^3 + 2X^2 + 3X + 5*Matrix{Int8}(I, 3, 3,)
f_s(X) = X^3 + 2X^2 + 3X + 5

fX = f_m(A)

λ1, λ2, λ3 = eigvals(A)

println("eigvals(f(X))=", eigvals(fX))
println("f(λ1)=", f_s(λ1),
        "f(λ2)=", f_s(λ2),
        "f(λ3)=", f_s(λ3))


text:
"""
    MCIndices(sz::NTuple{N,AbstractVector{Int}}) -> R
    MCIndices(A::AbstractArray) -> R

A `CartesianIndices` like type defines mutable and disconnected region `R`.

# Examples

```jldoctest
julia> im = MCIndices(([1, 3], [2, 4]))
2×2 MCIndices{2}:
 (1, 2)  (1, 4)
 (3, 2)  (3, 4)

julia> im[1]
(1, 2)

julia> im[1, 2]
(1, 4)
```
"""
struct MCIndices{N} <: AbstractArray{Dims{N},N}
    indices::NTuple{N,Vector{Int}}
end
MCIndices(indices::NTuple{N,AbstractVector{Int}}) where {N} =
    MCIndices(map(collect, indices))
MCIndices(A::AbstractArray) = MCIndices(axes(A))

Base.size(im::MCIndices) = map(length, im.indices)

Base.@propagate_inbounds getdim(im::MCIndices, i::Int) = im.indices[i]
Base.@propagate_inbounds Base.getindex(im::MCIndices{N}, I::Vararg{Int,N}) where {N} =
    map(getindex, im.indices, I)::Dims{N}

# DOKArray
const DOK{T,N} = Dict{Dims{N},T}

"""
    DOKSparseArray{T,N}

A simple sparse array type storing as
[`DOK` format](https://en.wikipedia.org/wiki/Sparse_matrix#Dictionary_of_keys_(DOK)).
"""
struct DOKSparseArray{T,N} <: ResizingTools.AbstractRNArray{T,N}
    dok::DOK{T,N}
    sz::Size{N}
end

Base.parent(A::DOKSparseArray) = A.dok
Base.isassigned(A::DOKSparseArray, i::Integer...) =
    haskey(A.dok, Base._to_subscript_indices(A, i...))
ArrayInterface.parent_type(::Type{<:DOKSparseArray{T,N}}) where {T,N} = DOK{T,N}
ResizingTools.isresizable(::Type{<:DOKSparseArray}) = true
ResizingTools.has_resize_buffer(::Type{<:DOKSparseArray}) = true
ResizingTools.getsize(A::DOKSparseArray{T,N}) where {T,N} = A.sz
ResizingTools.size_type(::Type{<:DOKSparseArray{T,N}}) where {T,N} = Size{N}

function Base.getindex(A::DOKSparseArray{T,N}, I::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A, I...)
    return get(A.dok, I, zero(T))
end
# this methods not used in this package
function Base.setindex!(A::DOKSparseArray{T,N}, v, I::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A, I...)
    return A.dok[I] = v
end

# Tools
function Base.get!(A::DOKSparseArray{T,N}, I::Dims{N}, v::T) where {T,N}
    @boundscheck checkbounds(A, I...)
    return get!(A.dok, I, v)
end
Base.get(A::DOKSparseArray{T,N}, I::Dims{N}, default) where {T,N} =
    get(parent(A), I, default)
Base.get(A::DOKSparseArray{T,0}, I::Tuple{}, default) where {T} = get(parent(A), I, default) # to avoid ambiguities

Base.keys(A::DOKSparseArray) = CartesianIndex.(keys(parent(A)))
Base.values(A::DOKSparseArray) = values(parent(A))

# these two methods don't create elements, the elements are created by setindex!
function ResizingTools.resize_buffer!(A::DOKSparseArray{T,N}, sz::Vararg{Any,N}) where {T,N}
    sz′ = to_dims(sz)
    @boundscheck all(map(<=, size(A), sz′)) ||
                 throw(ArgumentError("new size must large than the old one"))
    setsize!(A, sz′)
    return A
end
function ResizingTools.resize_buffer_dim!(A::DOKSparseArray, d::Int, n)
    n′ = to_dims(n)
    @boundscheck size(A, d) <= n′ ||
                 throw(ArgumentError("new size must large than the old one"))
    setsize!(A, d, n′)
    return A
end


text:
using YAXArrays
using Test

include("Cubes/axes.jl")
include("Cubes/cubes.jl")
include("Cubes/transformedcubes.jl")

include("Datasets/datasets.jl")


text:
using PICDataStructures, Test
using PICDataStructures: dimensionality, unitname
using LinearAlgebra
using RecursiveArrayTools: recursive_bottom_eltype
using Unitful
using StaticArrays

@testset "Vector field interface" begin
    grids =
        SparseAxisGrid(0.001:0.01:1),
        SparseAxisGrid(0.001:0.01:1, 0:0.01:1),
        SparseAxisGrid(0.001:0.005:1, 0:0.01:1, 0:0.01:1),
        SparseAxisGrid(0.001u"m":0.01u"m":1u"m"),
        SparseAxisGrid(0.001u"m":0.01u"m":1u"m", 0u"m":0.01u"m":1u"m"),
        SparseAxisGrid(0.001u"m":0.005u"m":1u"m", 0u"m":0.01u"m":1u"m", 0u"m":0.01u"m":1u"m")

    scalar_fields = [
        scalarfield(x->inv(x...), grids[1]),
        scalarfield(grids[2]) do (x,y)
            1 / √(x^2 + y^2)
        end,
        scalarfield(grids[3]) do (x,y,z)
            1 / √(x^2 + y^2 + z^2)
        end,
        scalarfield(grids[4]) do (x,)
            1u"V" / x
        end,
        scalarfield(grids[5]) do (x,y)
            1u"V" / √(x^2 + y^2)
        end,
        scalarfield(grids[6]) do (x,y,z)
            1u"V" / √(x^2 + y^2 + z^2)
        end
    ]
    fields = [
        build_vector((scalar_fields[1], ), (:x,), name="1D"),
        build_vector((scalar_fields[2], scalar_fields[2]), (:x, :y), name="2D"),
        build_vector((scalar_fields[3], scalar_fields[3], scalar_fields[3]), (:x, :y, :z), name="3D"),
        build_vector((scalar_fields[4], ), (:x,), name="1D"),
        build_vector((scalar_fields[5], scalar_fields[5]), (:x, :y), name="2D"),
        build_vector((scalar_fields[6], scalar_fields[6], scalar_fields[6]), (:x, :y, :z), name="3D")
    ]

    @testset "$(dimensionality(f))D ($(unitname(f)))" for (grid,s,f) in zip(grids,scalar_fields,fields)
        T = typeof(f)
        N = dimensionality(f)
        @test isconcretetype(T)

        @testset "Traits" begin
            @test scalarness(T) === VectorQuantity()
            @test domain_discretization(T) === LatticeGrid{N}()
            @test domain_type(T) <: SparseAxisGrid
        end

        @test f.x == s

        @testset "Indexing" begin
            @test size(f) == size(grid)
            t = recursive_bottom_eltype(f)
            @test all(f[1] .== t(1000))
            @test f[1] isa SVector
            @test f[1,:,:] == getfield(f, :data)[1,:,:]
            @test f[end][1] .≈ oneunit(t) / √N atol=t(1e-2)
        end

        @testset "Broadcasting" begin
            f2 = f .* 2
            @test f .* 1 == f
            @test typeof(f2) == typeof(f)
            @test getdomain(f2) == getdomain(f)
            @test nameof(f2) == nameof(f)
        end

        @testset "Downsampling" begin
            if N == 3
                f_small = downsample(f, 51, 26, 26)
                @test size(f_small) == (51, 26, 26)
                f_small = downsample(f)
                @test all(size(f_small) .≤ 15)
            elseif N == 2
                f_small = downsample(f, 51, 51)
                @test size(f_small) == (51, 51)
                f_small = downsample(f)
                @test all(size(f_small) .≤ 25)
            else
                f_small = downsample(f, 5)
                @test size(f_small) == (5,)
            end
            @test nameof(f) == nameof(f_small) == "$(N)D"
        end

        @testset "Slicing" begin
            if N > 1
                𝟘 = zero(recursive_bottom_eltype(grid))
                f_slice = selectdim(f, :x, 𝟘)
                @test length(propertynames(f_slice)) == N - 1
                @test dimensionality(f_slice) == N - 1
                @test nameof(f_slice) == nameof(f)

                f_slice2 = f[x=𝟘]
                @test dimensionality(f_slice) == N - 1

                f_slice = f[x=𝟘]
                d = f_slice .⋅ f_slice
                @test scalarness(typeof(d)) === ScalarQuantity()
            end
        end

        @testset "LinearAlgebra" begin
            nvf = norm.(f)
            @test scalarness(typeof(nvf)) === ScalarQuantity()
            @test norm.(f)[2] == nvf[2]
            # When converting to scalars from vectors we drop the name
            @test nameof(nvf) == ""

            if N == 3
                c = f × f
                @test all(iszero.(c))
                nameof(c) == "3D × 3D"
            end
        end
    end

    @testset "Unit handling" begin
        @test ustrip(fields[4]) == fields[1]
        @test ustrip(fields[5]) == fields[2]
        @test ustrip(fields[6]) == fields[3]

        @test fields[1] .* u"V/m" == fields[4]
        @test fields[2] .* u"V/m" == fields[5]
        @test fields[3] .* u"V/m" == fields[6]

        @test ustrip(u"V/m", fields[4]) == fields[1]
        @test ustrip(u"V/m", fields[5]) == fields[2]
        @test ustrip(u"V/m", fields[6]) == fields[3]

        @test ustrip(u"V/m", u"m", fields[4]) == fields[1]
        @test ustrip(u"V/m", u"m", fields[5]) == fields[2]
        @test ustrip(u"V/m", u"m", fields[6]) == fields[3]

        @test ustrip(fields[1]) == fields[1]
        @test nameof(ustrip(fields[1])) == "1D"
        @test nameof(ustrip(u"V/m", fields[4])) == "1D"
        @test nameof(ustrip(u"V/m", u"m", fields[4])) == "1D"
    end
end

@testset "Vector variable interface" begin
    grids = [
        ParticlePositions(sin.(0:0.1:2π)),
        ParticlePositions(sin.(0:0.1:2π),cos.(0:0.1:2π)),
        ParticlePositions(sin.(0:0.1:2π),cos.(0:0.1:2π),sin.(0:0.1:2π)),
        ParticlePositions(sin.(0:0.1:2π).*u"m"),
        ParticlePositions(sin.(0:0.1:2π).*u"m",cos.(0:0.1:2π).*u"m"),
        ParticlePositions(sin.(0:0.1:2π).*u"m",cos.(0:0.1:2π).*u"m",sin.(0:0.1:2π).*u"m"),
    ]
    scalar_vars = [
        scalarvariable(grids[1]) do (x,)
            x^2
        end,
        scalarvariable(grids[2]) do (x,y)
            (x)/10
        end,
        scalarvariable(grids[3]) do (x,y,z)
            (x^2 + y^2 + z^2)/10
        end,
        scalarvariable(grids[4]) do (x,)
            x^2
        end,
        scalarvariable(grids[5]) do (x,y)
            (x)/10
        end,
        scalarvariable(grids[6]) do (x,y,z)
            (x^2 + y^2 + z^2)/10
        end,
    ]
    vars = [
        build_vector((scalar_vars[1], ), (:x,), name="1D"),
        build_vector((scalar_vars[2], scalar_vars[2]), (:x, :y), name="2D"),
        build_vector((scalar_vars[3], scalar_vars[3], scalar_vars[3]), (:x, :y, :z), name="3D"),
        build_vector((scalar_vars[4], ), (:x,), name="1D"),
        build_vector((scalar_vars[5], scalar_vars[5]), (:x, :y), name="2D"),
        build_vector((scalar_vars[6], scalar_vars[6], scalar_vars[6]), (:x, :y, :z), name="3D")
    ]

    @testset "$(dimensionality(v))D ($(unitname(v)))" for (grid,f,v) in zip(grids,scalar_vars,vars)
        @test v isa VectorVariable
        @test getdomain(v) == getdomain(f) == grid

        T = typeof(v)
        N = dimensionality(v)
        @test isconcretetype(T)

        @test v.x == f

        @testset "Traits" begin
            @test scalarness(T) === VectorQuantity()
            @test domain_discretization(T) === ParticleGrid()
            @test domain_type(T) <: ParticlePositions
        end

        @testset "Indexing" begin
            @test size(v) == size(grid)
            @test v[1] isa SVector
            @test v[1,:] == getfield(v, :data)[1,:]
        end

        @testset "Broadcasting" begin
            bc = Base.broadcasted(x->x*2, v)
            ibc = Base.Broadcast.instantiate(bc)
            ElType = Base.Broadcast.combine_eltypes(ibc.f, ibc.args)
            @test ElType == eltype(v)

            v2 = v .* 2
            @test v .* 1 == v
            @test typeof(v2) == typeof(v)
            @test getdomain(v2) == getdomain(v)
            @test nameof(v2) == nameof(v)
        end

        @testset "Downsampling" begin
            v_small = downsample(v, 3)
            @test size(v_small) == (3,)
            v_small = downsample(v)
            @test all(size(v_small) .≤ size(v))
            @test nameof(v) == nameof(v_small)
        end

        @testset "LinearAlgebra" begin
            nvf = norm.(v)
            @test scalarness(typeof(nvf)) === ScalarQuantity()
            @test norm.(v)[2] == nvf[2]
            # When converting to scalars from vectors we drop the name
            @test nameof(nvf) == ""

            if N == 3
                c = v × v
                @test all(iszero.(c))
                nameof(c) == "3D × 3D"
            end
        end
    end

    @testset "Sclicing" begin
        @testset "2D" begin
            v = vars[2]
            v_slice = selectdim(v, :x, 0.0, ϵ=0.1)
            @test dimensionality(v_slice) == 1
            @test propertynames(v_slice) == (:y,)
            @test size(v_slice) == size(getdomain(v_slice))
            @test length(v_slice) == 5
            @test all(isapprox.(only.(v_slice), 0.0, atol=1e-2))
            @test propertynames(getdomain(v_slice)) == (:y,)
            # TODO: make this less fragile
            @test all(isapprox.(getdomain(v_slice).y[[1,2,5]], 1.0, atol=1e-2))
            @test all(isapprox.(getdomain(v_slice).y[3:4], -1.0, atol=1e-2))

            @test nameof(v) == nameof(v_slice) == "2D"
        end
        @testset "3D" begin
            v = vars[3]
            v_slice = selectdim(v, :x, 0.0, ϵ=0.1)
            @test dimensionality(v_slice) == 2
            @test propertynames(v_slice) == (:y,:z)
            @test size(v_slice) == size(getdomain(v_slice))
            @test length(v_slice) == 5
            @test all(isapprox.(v_slice, ([0.1, 0.1],), atol=1e-2))
            @test propertynames(getdomain(v_slice)) == (:y,:z)
            # TODO: make this less fragile
            @test all(isapprox.(getdomain(v_slice).y[[1,2,5]], 1.0, atol=1e-2))
            @test all(isapprox.(getdomain(v_slice).y[3:4], -1.0, atol=1e-2))
            @test count(getdomain(v_slice).z .< 0.0) == 2
            @test iszero(first(getdomain(v_slice).z))

            @test nameof(v) == nameof(v_slice) == "3D"
        end
    end

    @testset "Units" begin
        @test ustrip(vars[4]) == vars[1]
        @test vars[1] .* 1u"m^2" == vars[4]

        @test ustrip(vars[5]) == vars[2]
        @test vars[2] .* 1u"m" == vars[5]

        @test ustrip(vars[6]) == vars[3]
        @test vars[3] .* 1u"m^2" == vars[6]

        @test nameof(ustrip(vars[1])) == "1D"
    end
end


text:
module HIP

import ..AMDGPU.libhip
using CEnum

include("libhip_common.jl")
include("error.jl")
include("libhip.jl")

end