text:
function list_points(mapping, grid, mesh, shift, ir_only)
    shift = shift ./ 2  # true / 2 = 0.5, false / 2 = 0
    if ir_only
        return map(unique(mapping)) do i
            (grid[:, i] .+ shift) ./ mesh
        end
    else
        return map(eachslice(grid; dims = 2)) do point
            (point .+ shift) ./ mesh  # Add 1 because `mapping` index starts from 0
        end
    end
end

# From https://github.com/unkcpz/LibSymspg.jl/blob/53d2f6d/test/test_api.jl#L89-L99
@testset "Test reciprocal mesh" begin
    lattice = [
        -2.0 2.0 2.0
        2.0 -2.0 2.0
        2.0 2.0 -2.0
    ]
    positions = [0.0 0.0 0.0]'
    types = [1]
    mesh = [4, 4, 4]
    is_shift = [0, 0, 0]
    cell = Cell(lattice, positions, types)
    nir, mapping, grid_address = get_ir_reciprocal_mesh(
        cell,
        mesh,
        is_shift;
        is_time_reversal = true,
        symprec = 1e-5,
    )
    @test nir == length(unique(mapping)) == 8
end

# See https://spglib.github.io/spglib/python-spglib.html
@testset "Python example" begin
    lattice = [[0.0, 0.5, 0.5], [0.5, 0.0, 0.5], [0.5, 0.5, 0.0]] * 5.4
    positions = [[0.875, 0.875, 0.875], [0.125, 0.125, 0.125]]
    types = [1, 1]
    cell = Cell(lattice, positions, types)
    mesh = [8, 8, 8]
    @test get_spacegroup_type(cell).international_short == "Fd-3m"
    @test get_spacegroup_type(cell).number == 227
    @testset "No shifts" begin
        shift = [0, 0, 0]
        nir, mapping, grid_address = get_ir_reciprocal_mesh(
            cell,
            mesh,
            shift;
            is_time_reversal = true,
            symprec = 1e-5,
        )
        python_mapping = [
            0,
            1,
            2,
            3,
            4,
            3,
            2,
            1,
            1,
            9,
            10,
            11,
            12,
            13,
            14,
            15,
            2,
            10,
            18,
            19,
            20,
            21,
            22,
            14,
            3,
            11,
            19,
            27,
            28,
            29,
            21,
            13,
            4,
            12,
            20,
            28,
            36,
            28,
            20,
            12,
            3,
            13,
            21,
            29,
            28,
            27,
            19,
            11,
            2,
            14,
            22,
            21,
            20,
            19,
            18,
            10,
            1,
            15,
            14,
            13,
            12,
            11,
            10,
            9,
            1,
            9,
            10,
            11,
            12,
            13,
            14,
            15,
            9,
            1,
            15,
            14,
            13,
            12,
            11,
            10,
            10,
            15,
            10,
            83,
            84,
            85,
            84,
            83,
            11,
            14,
            83,
            19,
            92,
            93,
            94,
            84,
            12,
            13,
            84,
            92,
            28,
            101,
            93,
            85,
            13,
            12,
            85,
            93,
            101,
            28,
            92,
            84,
            14,
            11,
            84,
            94,
            93,
            92,
            19,
            83,
            15,
            10,
            83,
            84,
            85,
            84,
            83,
            10,
            2,
            10,
            18,
            19,
            20,
            21,
            22,
            14,
            10,
            15,
            10,
            83,
            84,
            85,
            84,
            83,
            18,
            10,
            2,
            14,
            22,
            21,
            20,
            19,
            19,
            83,
            14,
            11,
            84,
            94,
            93,
            92,
            20,
            84,
            22,
            84,
            20,
            93,
            166,
            93,
            21,
            85,
            21,
            94,
            93,
            29,
            93,
            94,
            22,
            84,
            20,
            93,
            166,
            93,
            20,
            84,
            14,
            83,
            19,
            92,
            93,
            94,
            84,
            11,
            3,
            11,
            19,
            27,
            28,
            29,
            21,
            13,
            11,
            14,
            83,
            19,
            92,
            93,
            94,
            84,
            19,
            83,
            14,
            11,
            84,
            94,
            93,
            92,
            27,
            19,
            11,
            3,
            13,
            21,
            29,
            28,
            28,
            92,
            84,
            13,
            12,
            85,
            93,
            101,
            29,
            93,
            94,
            21,
            85,
            21,
            94,
            93,
            21,
            94,
            93,
            29,
            93,
            94,
            21,
            85,
            13,
            84,
            92,
            28,
            101,
            93,
            85,
            12,
            4,
            12,
            20,
            28,
            36,
            28,
            20,
            12,
            12,
            13,
            84,
            92,
            28,
            101,
            93,
            85,
            20,
            84,
            22,
            84,
            20,
            93,
            166,
            93,
            28,
            92,
            84,
            13,
            12,
            85,
            93,
            101,
            36,
            28,
            20,
            12,
            4,
            12,
            20,
            28,
            28,
            101,
            93,
            85,
            12,
            13,
            84,
            92,
            20,
            93,
            166,
            93,
            20,
            84,
            22,
            84,
            12,
            85,
            93,
            101,
            28,
            92,
            84,
            13,
            3,
            13,
            21,
            29,
            28,
            27,
            19,
            11,
            13,
            12,
            85,
            93,
            101,
            28,
            92,
            84,
            21,
            85,
            21,
            94,
            93,
            29,
            93,
            94,
            29,
            93,
            94,
            21,
            85,
            21,
            94,
            93,
            28,
            101,
            93,
            85,
            12,
            13,
            84,
            92,
            27,
            28,
            29,
            21,
            13,
            3,
            11,
            19,
            19,
            92,
            93,
            94,
            84,
            11,
            14,
            83,
            11,
            84,
            94,
            93,
            92,
            19,
            83,
            14,
            2,
            14,
            22,
            21,
            20,
            19,
            18,
            10,
            14,
            11,
            84,
            94,
            93,
            92,
            19,
            83,
            22,
            84,
            20,
            93,
            166,
            93,
            20,
            84,
            21,
            94,
            93,
            29,
            93,
            94,
            21,
            85,
            20,
            93,
            166,
            93,
            20,
            84,
            22,
            84,
            19,
            92,
            93,
            94,
            84,
            11,
            14,
            83,
            18,
            19,
            20,
            21,
            22,
            14,
            2,
            10,
            10,
            83,
            84,
            85,
            84,
            83,
            10,
            15,
            1,
            15,
            14,
            13,
            12,
            11,
            10,
            9,
            15,
            10,
            83,
            84,
            85,
            84,
            83,
            10,
            14,
            83,
            19,
            92,
            93,
            94,
            84,
            11,
            13,
            84,
            92,
            28,
            101,
            93,
            85,
            12,
            12,
            85,
            93,
            101,
            28,
            92,
            84,
            13,
            11,
            84,
            94,
            93,
            92,
            19,
            83,
            14,
            10,
            83,
            84,
            85,
            84,
            83,
            10,
            15,
            9,
            10,
            11,
            12,
            13,
            14,
            15,
            1,
        ]
        @test mapping == python_mapping .+ 1
        @test nir == length(unique(mapping)) == 29
        # Irreducible k-points
        python_results = [
            [0.0, 0.0, 0.0],
            [0.125, 0.0, 0.0],
            [0.25, 0.0, 0.0],
            [0.375, 0.0, 0.0],
            [0.5, 0.0, 0.0],
            [0.125, 0.125, 0.0],
            [0.25, 0.125, 0.0],
            [0.375, 0.125, 0.0],
            [0.5, 0.125, 0.0],
            [-0.375, 0.125, 0.0],
            [-0.25, 0.125, 0.0],
            [-0.125, 0.125, 0.0],
            [0.25, 0.25, 0.0],
            [0.375, 0.25, 0.0],
            [0.5, 0.25, 0.0],
            [-0.375, 0.25, 0.0],
            [-0.25, 0.25, 0.0],
            [0.375, 0.375, 0.0],
            [0.5, 0.375, 0.0],
            [-0.375, 0.375, 0.0],
            [0.5, 0.5, 0.0],
            [0.375, 0.25, 0.125],
            [0.5, 0.25, 0.125],
            [-0.375, 0.25, 0.125],
            [0.5, 0.375, 0.125],
            [-0.375, 0.375, 0.125],
            [-0.25, 0.375, 0.125],
            [-0.375, 0.5, 0.125],
            [-0.25, 0.5, 0.25],
        ]
        @test list_points(mapping, grid_address, mesh, shift, true) == python_results
        @test length(list_points(mapping, grid_address, mesh, shift, false)) == 512
    end
    @testset "With shifts" begin
        shift = [1, 1, 1]
        nir, mapping, grid_address = get_ir_reciprocal_mesh(
            cell,
            mesh,
            shift;
            is_time_reversal = true,
            symprec = 1e-5,
        )
        python_mapping = [
            0,
            1,
            2,
            3,
            4,
            5,
            6,
            7,
            1,
            9,
            10,
            11,
            12,
            13,
            14,
            15,
            2,
            10,
            18,
            19,
            20,
            21,
            22,
            23,
            3,
            11,
            19,
            27,
            28,
            29,
            30,
            31,
            4,
            12,
            20,
            28,
            36,
            37,
            38,
            31,
            5,
            13,
            21,
            29,
            37,
            45,
            46,
            23,
            6,
            14,
            22,
            30,
            38,
            46,
            54,
            15,
            7,
            15,
            23,
            31,
            31,
            23,
            15,
            7,
            1,
            9,
            10,
            11,
            12,
            13,
            14,
            15,
            9,
            73,
            74,
            75,
            76,
            77,
            78,
            54,
            10,
            74,
            82,
            83,
            84,
            85,
            86,
            46,
            11,
            75,
            83,
            91,
            92,
            93,
            94,
            38,
            12,
            76,
            84,
            92,
            100,
            101,
            94,
            30,
            13,
            77,
            85,
            93,
            101,
            109,
            86,
            22,
            14,
            78,
            86,
            94,
            94,
            86,
            78,
            14,
            15,
            54,
            46,
            38,
            30,
            22,
            14,
            6,
            2,
            10,
            18,
            19,
            20,
            21,
            22,
            23,
            10,
            74,
            82,
            83,
            84,
            85,
            86,
            46,
            18,
            82,
            146,
            147,
            148,
            149,
            109,
            45,
            19,
            83,
            147,
            155,
            156,
            157,
            101,
            37,
            20,
            84,
            148,
            156,
            164,
            157,
            93,
            29,
            21,
            85,
            149,
            157,
            157,
            149,
            85,
            21,
            22,
            86,
            109,
            101,
            93,
            85,
            77,
            13,
            23,
            46,
            45,
            37,
            29,
            21,
            13,
            5,
            3,
            11,
            19,
            27,
            28,
            29,
            30,
            31,
            11,
            75,
            83,
            91,
            92,
            93,
            94,
            38,
            19,
            83,
            147,
            155,
            156,
            157,
            101,
            37,
            27,
            91,
            155,
            219,
            220,
            164,
            100,
            36,
            28,
            92,
            156,
            220,
            220,
            156,
            92,
            28,
            29,
            93,
            157,
            164,
            156,
            148,
            84,
            20,
            30,
            94,
            101,
            100,
            92,
            84,
            76,
            12,
            31,
            38,
            37,
            36,
            28,
            20,
            12,
            4,
            4,
            12,
            20,
            28,
            36,
            37,
            38,
            31,
            12,
            76,
            84,
            92,
            100,
            101,
            94,
            30,
            20,
            84,
            148,
            156,
            164,
            157,
            93,
            29,
            28,
            92,
            156,
            220,
            220,
            156,
            92,
            28,
            36,
            100,
            164,
            220,
            219,
            155,
            91,
            27,
            37,
            101,
            157,
            156,
            155,
            147,
            83,
            19,
            38,
            94,
            93,
            92,
            91,
            83,
            75,
            11,
            31,
            30,
            29,
            28,
            27,
            19,
            11,
            3,
            5,
            13,
            21,
            29,
            37,
            45,
            46,
            23,
            13,
            77,
            85,
            93,
            101,
            109,
            86,
            22,
            21,
            85,
            149,
            157,
            157,
            149,
            85,
            21,
            29,
            93,
            157,
            164,
            156,
            148,
            84,
            20,
            37,
            101,
            157,
            156,
            155,
            147,
            83,
            19,
            45,
            109,
            149,
            148,
            147,
            146,
            82,
            18,
            46,
            86,
            85,
            84,
            83,
            82,
            74,
            10,
            23,
            22,
            21,
            20,
            19,
            18,
            10,
            2,
            6,
            14,
            22,
            30,
            38,
            46,
            54,
            15,
            14,
            78,
            86,
            94,
            94,
            86,
            78,
            14,
            22,
            86,
            109,
            101,
            93,
            85,
            77,
            13,
            30,
            94,
            101,
            100,
            92,
            84,
            76,
            12,
            38,
            94,
            93,
            92,
            91,
            83,
            75,
            11,
            46,
            86,
            85,
            84,
            83,
            82,
            74,
            10,
            54,
            78,
            77,
            76,
            75,
            74,
            73,
            9,
            15,
            14,
            13,
            12,
            11,
            10,
            9,
            1,
            7,
            15,
            23,
            31,
            31,
            23,
            15,
            7,
            15,
            54,
            46,
            38,
            30,
            22,
            14,
            6,
            23,
            46,
            45,
            37,
            29,
            21,
            13,
            5,
            31,
            38,
            37,
            36,
            28,
            20,
            12,
            4,
            31,
            30,
            29,
            28,
            27,
            19,
            11,
            3,
            23,
            22,
            21,
            20,
            19,
            18,
            10,
            2,
            15,
            14,
            13,
            12,
            11,
            10,
            9,
            1,
            7,
            6,
            5,
            4,
            3,
            2,
            1,
            0,
        ]
        @test mapping == python_mapping .+ 1
        @test nir == length(unique(mapping)) == 60
        # Irreducible k-points
        python_results = [
            [0.0625, 0.0625, 0.0625],
            [0.1875, 0.0625, 0.0625],
            [0.3125, 0.0625, 0.0625],
            [0.4375, 0.0625, 0.0625],
            [0.5625, 0.0625, 0.0625],
            [-0.3125, 0.0625, 0.0625],
            [-0.1875, 0.0625, 0.0625],
            [-0.0625, 0.0625, 0.0625],
            [0.1875, 0.1875, 0.0625],
            [0.3125, 0.1875, 0.0625],
            [0.4375, 0.1875, 0.0625],
            [0.5625, 0.1875, 0.0625],
            [-0.3125, 0.1875, 0.0625],
            [-0.1875, 0.1875, 0.0625],
            [-0.0625, 0.1875, 0.0625],
            [0.3125, 0.3125, 0.0625],
            [0.4375, 0.3125, 0.0625],
            [0.5625, 0.3125, 0.0625],
            [-0.3125, 0.3125, 0.0625],
            [-0.1875, 0.3125, 0.0625],
            [-0.0625, 0.3125, 0.0625],
            [0.4375, 0.4375, 0.0625],
            [0.5625, 0.4375, 0.0625],
            [-0.3125, 0.4375, 0.0625],
            [-0.1875, 0.4375, 0.0625],
            [-0.0625, 0.4375, 0.0625],
            [0.5625, 0.5625, 0.0625],
            [-0.3125, 0.5625, 0.0625],
            [-0.1875, 0.5625, 0.0625],
            [-0.3125, -0.3125, 0.0625],
            [-0.1875, -0.3125, 0.0625],
            [-0.1875, -0.1875, 0.0625],
            [0.1875, 0.1875, 0.1875],
            [0.3125, 0.1875, 0.1875],
            [0.4375, 0.1875, 0.1875],
            [0.5625, 0.1875, 0.1875],
            [-0.3125, 0.1875, 0.1875],
            [-0.1875, 0.1875, 0.1875],
            [0.3125, 0.3125, 0.1875],
            [0.4375, 0.3125, 0.1875],
            [0.5625, 0.3125, 0.1875],
            [-0.3125, 0.3125, 0.1875],
            [-0.1875, 0.3125, 0.1875],
            [0.4375, 0.4375, 0.1875],
            [0.5625, 0.4375, 0.1875],
            [-0.3125, 0.4375, 0.1875],
            [-0.1875, 0.4375, 0.1875],
            [0.5625, 0.5625, 0.1875],
            [-0.3125, 0.5625, 0.1875],
            [-0.3125, -0.3125, 0.1875],
            [0.3125, 0.3125, 0.3125],
            [0.4375, 0.3125, 0.3125],
            [0.5625, 0.3125, 0.3125],
            [-0.3125, 0.3125, 0.3125],
            [0.4375, 0.4375, 0.3125],
            [0.5625, 0.4375, 0.3125],
            [-0.3125, 0.4375, 0.3125],
            [0.5625, 0.5625, 0.3125],
            [0.4375, 0.4375, 0.4375],
            [0.5625, 0.4375, 0.4375],
        ]
        @test list_points(mapping, grid_address, mesh, shift, true) == python_results
        @test length(list_points(mapping, grid_address, mesh, shift, false)) == 512
    end
end

@testset "Test rutile structure" begin
    lattice = [
        4 0 0
        0 4 0
        0 0 3
    ]
    positions = [
        0.0 0.5 0.3 0.7 0.2 0.8
        0.0 0.5 0.3 0.7 0.8 0.2
        0.0 0.5 0.0 0.0 0.5 0.5
    ]
    types = [14, 14, 8, 8, 8, 8]
    rutile = Cell(lattice, positions, types)
    @testset "No shift" begin
        shift = [0, 0, 0]
        mesh = [6, 6, 6]
        nir, mapping, grid_address =
            get_ir_reciprocal_mesh(rutile, mesh; is_time_reversal = true, symprec = 1e-5)
        @test nir == length(unique(mapping)) == 40
        @test list_points(mapping, grid_address, mesh, shift, true) ≈ [
            [0.0, 0.0, 0.0],
            [0.16666667, 0.0, 0.0],
            [0.33333333, 0.0, 0.0],
            [0.5, 0.0, 0.0],
            [0.16666667, 0.16666667, 0.0],
            [0.33333333, 0.16666667, 0.0],
            [0.5, 0.16666667, 0.0],
            [0.33333333, 0.33333333, 0.0],
            [0.5, 0.33333333, 0.0],
            [0.5, 0.5, 0.0],
            [0.0, 0.0, 0.16666667],
            [0.16666667, 0.0, 0.16666667],
            [0.33333333, 0.0, 0.16666667],
            [0.5, 0.0, 0.16666667],
            [0.16666667, 0.16666667, 0.16666667],
            [0.33333333, 0.16666667, 0.16666667],
            [0.5, 0.16666667, 0.16666667],
            [0.33333333, 0.33333333, 0.16666667],
            [0.5, 0.33333333, 0.16666667],
            [0.5, 0.5, 0.16666667],
            [0.0, 0.0, 0.33333333],
            [0.16666667, 0.0, 0.33333333],
            [0.33333333, 0.0, 0.33333333],
            [0.5, 0.0, 0.33333333],
            [0.16666667, 0.16666667, 0.33333333],
            [0.33333333, 0.16666667, 0.33333333],
            [0.5, 0.16666667, 0.33333333],
            [0.33333333, 0.33333333, 0.33333333],
            [0.5, 0.33333333, 0.33333333],
            [0.5, 0.5, 0.33333333],
            [0.0, 0.0, 0.5],
            [0.16666667, 0.0, 0.5],
            [0.33333333, 0.0, 0.5],
            [0.5, 0.0, 0.5],
            [0.16666667, 0.16666667, 0.5],
            [0.33333333, 0.16666667, 0.5],
            [0.5, 0.16666667, 0.5],
            [0.33333333, 0.33333333, 0.5],
            [0.5, 0.33333333, 0.5],
            [0.5, 0.5, 0.5],
        ]
        @test list_points(mapping, grid_address, mesh, shift, false) ≈ [
            [0.0, 0.0, 0.0],
            [0.16666667, 0.0, 0.0],
            [0.33333333, 0.0, 0.0],
            [0.5, 0.0, 0.0],
            [-0.33333333, 0.0, 0.0],
            [-0.16666667, 0.0, 0.0],
            [0.0, 0.16666667, 0.0],
            [0.16666667, 0.16666667, 0.0],
            [0.33333333, 0.16666667, 0.0],
            [0.5, 0.16666667, 0.0],
            [-0.33333333, 0.16666667, 0.0],
            [-0.16666667, 0.16666667, 0.0],
            [0.0, 0.33333333, 0.0],
            [0.16666667, 0.33333333, 0.0],
            [0.33333333, 0.33333333, 0.0],
            [0.5, 0.33333333, 0.0],
            [-0.33333333, 0.33333333, 0.0],
            [-0.16666667, 0.33333333, 0.0],
            [0.0, 0.5, 0.0],
            [0.16666667, 0.5, 0.0],
            [0.33333333, 0.5, 0.0],
            [0.5, 0.5, 0.0],
            [-0.33333333, 0.5, 0.0],
            [-0.16666667, 0.5, 0.0],
            [0.0, -0.33333333, 0.0],
            [0.16666667, -0.33333333, 0.0],
            [0.33333333, -0.33333333, 0.0],
            [0.5, -0.33333333, 0.0],
            [-0.33333333, -0.33333333, 0.0],
            [-0.16666667, -0.33333333, 0.0],
            [0.0, -0.16666667, 0.0],
            [0.16666667, -0.16666667, 0.0],
            [0.33333333, -0.16666667, 0.0],
            [0.5, -0.16666667, 0.0],
            [-0.33333333, -0.16666667, 0.0],
            [-0.16666667, -0.16666667, 0.0],
            [0.0, 0.0, 0.16666667],
            [0.16666667, 0.0, 0.16666667],
            [0.33333333, 0.0, 0.16666667],
            [0.5, 0.0, 0.16666667],
            [-0.33333333, 0.0, 0.16666667],
            [-0.16666667, 0.0, 0.16666667],
            [0.0, 0.16666667, 0.16666667],
            [0.16666667, 0.16666667, 0.16666667],
            [0.33333333, 0.16666667, 0.16666667],
            [0.5, 0.16666667, 0.16666667],
            [-0.33333333, 0.16666667, 0.16666667],
            [-0.16666667, 0.16666667, 0.16666667],
            [0.0, 0.33333333, 0.16666667],
            [0.16666667, 0.33333333, 0.16666667],
            [0.33333333, 0.33333333, 0.16666667],
            [0.5, 0.33333333, 0.16666667],
            [-0.33333333, 0.33333333, 0.16666667],
            [-0.16666667, 0.33333333, 0.16666667],
            [0.0, 0.5, 0.16666667],
            [0.16666667, 0.5, 0.16666667],
            [0.33333333, 0.5, 0.16666667],
            [0.5, 0.5, 0.16666667],
            [-0.33333333, 0.5, 0.16666667],
            [-0.16666667, 0.5, 0.16666667],
            [0.0, -0.33333333, 0.16666667],
            [0.16666667, -0.33333333, 0.16666667],
            [0.33333333, -0.33333333, 0.16666667],
            [0.5, -0.33333333, 0.16666667],
            [-0.33333333, -0.33333333, 0.16666667],
            [-0.16666667, -0.33333333, 0.16666667],
            [0.0, -0.16666667, 0.16666667],
            [0.16666667, -0.16666667, 0.16666667],
            [0.33333333, -0.16666667, 0.16666667],
            [0.5, -0.16666667, 0.16666667],
            [-0.33333333, -0.16666667, 0.16666667],
            [-0.16666667, -0.16666667, 0.16666667],
            [0.0, 0.0, 0.33333333],
            [0.16666667, 0.0, 0.33333333],
            [0.33333333, 0.0, 0.33333333],
            [0.5, 0.0, 0.33333333],
            [-0.33333333, 0.0, 0.33333333],
            [-0.16666667, 0.0, 0.33333333],
            [0.0, 0.16666667, 0.33333333],
            [0.16666667, 0.16666667, 0.33333333],
            [0.33333333, 0.16666667, 0.33333333],
            [0.5, 0.16666667, 0.33333333],
            [-0.33333333, 0.16666667, 0.33333333],
            [-0.16666667, 0.16666667, 0.33333333],
            [0.0, 0.33333333, 0.33333333],
            [0.16666667, 0.33333333, 0.33333333],
            [0.33333333, 0.33333333, 0.33333333],
            [0.5, 0.33333333, 0.33333333],
            [-0.33333333, 0.33333333, 0.33333333],
            [-0.16666667, 0.33333333, 0.33333333],
            [0.0, 0.5, 0.33333333],
            [0.16666667, 0.5, 0.33333333],
            [0.33333333, 0.5, 0.33333333],
            [0.5, 0.5, 0.33333333],
            [-0.33333333, 0.5, 0.33333333],
            [-0.16666667, 0.5, 0.33333333],
            [0.0, -0.33333333, 0.33333333],
            [0.16666667, -0.33333333, 0.33333333],
            [0.33333333, -0.33333333, 0.33333333],
            [0.5, -0.33333333, 0.33333333],
            [-0.33333333, -0.33333333, 0.33333333],
            [-0.16666667, -0.33333333, 0.33333333],
            [0.0, -0.16666667, 0.33333333],
            [0.16666667, -0.16666667, 0.33333333],
            [0.33333333, -0.16666667, 0.33333333],
            [0.5, -0.16666667, 0.33333333],
            [-0.33333333, -0.16666667, 0.33333333],
            [-0.16666667, -0.16666667, 0.33333333],
            [0.0, 0.0, 0.5],
            [0.16666667, 0.0, 0.5],
            [0.33333333, 0.0, 0.5],
            [0.5, 0.0, 0.5],
            [-0.33333333, 0.0, 0.5],
            [-0.16666667, 0.0, 0.5],
            [0.0, 0.16666667, 0.5],
            [0.16666667, 0.16666667, 0.5],
            [0.33333333, 0.16666667, 0.5],
            [0.5, 0.16666667, 0.5],
            [-0.33333333, 0.16666667, 0.5],
            [-0.16666667, 0.16666667, 0.5],
            [0.0, 0.33333333, 0.5],
            [0.16666667, 0.33333333, 0.5],
            [0.33333333, 0.33333333, 0.5],
            [0.5, 0.33333333, 0.5],
            [-0.33333333, 0.33333333, 0.5],
            [-0.16666667, 0.33333333, 0.5],
            [0.0, 0.5, 0.5],
            [0.16666667, 0.5, 0.5],
            [0.33333333, 0.5, 0.5],
            [0.5, 0.5, 0.5],
            [-0.33333333, 0.5, 0.5],
            [-0.16666667, 0.5, 0.5],
            [0.0, -0.33333333, 0.5],
            [0.16666667, -0.33333333, 0.5],
            [0.33333333, -0.33333333, 0.5],
            [0.5, -0.33333333, 0.5],
            [-0.33333333, -0.33333333, 0.5],
            [-0.16666667, -0.33333333, 0.5],
            [0.0, -0.16666667, 0.5],
            [0.16666667, -0.16666667, 0.5],
            [0.33333333, -0.16666667, 0.5],
            [0.5, -0.16666667, 0.5],
            [-0.33333333, -0.16666667, 0.5],
            [-0.16666667, -0.16666667, 0.5],
            [0.0, 0.0, -0.33333333],
            [0.16666667, 0.0, -0.33333333],
            [0.33333333, 0.0, -0.33333333],
            [0.5, 0.0, -0.33333333],
            [-0.33333333, 0.0, -0.33333333],
            [-0.16666667, 0.0, -0.33333333],
            [0.0, 0.16666667, -0.33333333],
            [0.16666667, 0.16666667, -0.33333333],
            [0.33333333, 0.16666667, -0.33333333],
            [0.5, 0.16666667, -0.33333333],
            [-0.33333333, 0.16666667, -0.33333333],
            [-0.16666667, 0.16666667, -0.33333333],
            [0.0, 0.33333333, -0.33333333],
            [0.16666667, 0.33333333, -0.33333333],
            [0.33333333, 0.33333333, -0.33333333],
            [0.5, 0.33333333, -0.33333333],
            [-0.33333333, 0.33333333, -0.33333333],
            [-0.16666667, 0.33333333, -0.33333333],
            [0.0, 0.5, -0.33333333],
            [0.16666667, 0.5, -0.33333333],
            [0.33333333, 0.5, -0.33333333],
            [0.5, 0.5, -0.33333333],
            [-0.33333333, 0.5, -0.33333333],
            [-0.16666667, 0.5, -0.33333333],
            [0.0, -0.33333333, -0.33333333],
            [0.16666667, -0.33333333, -0.33333333],
            [0.33333333, -0.33333333, -0.33333333],
            [0.5, -0.33333333, -0.33333333],
            [-0.33333333, -0.33333333, -0.33333333],
            [-0.16666667, -0.33333333, -0.33333333],
            [0.0, -0.16666667, -0.33333333],
            [0.16666667, -0.16666667, -0.33333333],
            [0.33333333, -0.16666667, -0.33333333],
            [0.5, -0.16666667, -0.33333333],
            [-0.33333333, -0.16666667, -0.33333333],
            [-0.16666667, -0.16666667, -0.33333333],
            [0.0, 0.0, -0.16666667],
            [0.16666667, 0.0, -0.16666667],
            [0.33333333, 0.0, -0.16666667],
            [0.5, 0.0, -0.16666667],
            [-0.33333333, 0.0, -0.16666667],
            [-0.16666667, 0.0, -0.16666667],
            [0.0, 0.16666667, -0.16666667],
            [0.16666667, 0.16666667, -0.16666667],
            [0.33333333, 0.16666667, -0.16666667],
            [0.5, 0.16666667, -0.16666667],
            [-0.33333333, 0.16666667, -0.16666667],
            [-0.16666667, 0.16666667, -0.16666667],
            [0.0, 0.33333333, -0.16666667],
            [0.16666667, 0.33333333, -0.16666667],
            [0.33333333, 0.33333333, -0.16666667],
            [0.5, 0.33333333, -0.16666667],
            [-0.33333333, 0.33333333, -0.16666667],
            [-0.16666667, 0.33333333, -0.16666667],
            [0.0, 0.5, -0.16666667],
            [0.16666667, 0.5, -0.16666667],
            [0.33333333, 0.5, -0.16666667],
            [0.5, 0.5, -0.16666667],
            [-0.33333333, 0.5, -0.16666667],
            [-0.16666667, 0.5, -0.16666667],
            [0.0, -0.33333333, -0.16666667],
            [0.16666667, -0.33333333, -0.16666667],
            [0.33333333, -0.33333333, -0.16666667],
            [0.5, -0.33333333, -0.16666667],
            [-0.33333333, -0.33333333, -0.16666667],
            [-0.16666667, -0.33333333, -0.16666667],
            [0.0, -0.16666667, -0.16666667],
            [0.16666667, -0.16666667, -0.16666667],
            [0.33333333, -0.16666667, -0.16666667],
            [0.5, -0.16666667, -0.16666667],
            [-0.33333333, -0.16666667, -0.16666667],
            [-0.16666667, -0.16666667, -0.16666667],
        ]
    end
    @testset "With shift" begin
        shift = trues(3)
        @testset "6×6×6" begin
            mesh = [6, 6, 6]
            nir, mapping, grid_address = get_ir_reciprocal_mesh(
                rutile,
                mesh,
                shift;
                is_time_reversal = true,
                symprec = 1e-5,
            )
            @test nir == length(unique(mapping)) == 18
            @test list_points(mapping, grid_address, mesh, shift, true) ≈ [
                [0.08333333, 0.08333333, 0.08333333],
                [0.25, 0.08333333, 0.08333333],
                [0.41666667, 0.08333333, 0.08333333],
                [0.25, 0.25, 0.08333333],
                [0.41666667, 0.25, 0.08333333],
                [0.41666667, 0.41666667, 0.08333333],
                [0.08333333, 0.08333333, 0.25],
                [0.25, 0.08333333, 0.25],
                [0.41666667, 0.08333333, 0.25],
                [0.25, 0.25, 0.25],
                [0.41666667, 0.25, 0.25],
                [0.41666667, 0.41666667, 0.25],
                [0.08333333, 0.08333333, 0.41666667],
                [0.25, 0.08333333, 0.41666667],
                [0.41666667, 0.08333333, 0.41666667],
                [0.25, 0.25, 0.41666667],
                [0.41666667, 0.25, 0.41666667],
                [0.41666667, 0.41666667, 0.41666667],
            ]
        end
        @testset "4×4×4" begin
            mesh = [4, 4, 4]
            nir, mapping, grid_address = get_ir_reciprocal_mesh(
                rutile,
                mesh,
                shift;
                is_time_reversal = true,
                symprec = 1e-5,
            )
            @test nir == length(unique(mapping)) == 6
            @test list_points(mapping, grid_address, mesh, shift, true) == [
                [0.125, 0.125, 0.125],
                [0.375, 0.125, 0.125],
                [0.375, 0.375, 0.125],
                [0.125, 0.125, 0.375],
                [0.375, 0.125, 0.375],
                [0.375, 0.375, 0.375],
            ]
            @test list_points(mapping, grid_address, mesh, shift, false) == [
                [0.125, 0.125, 0.125],
                [0.375, 0.125, 0.125],
                [0.625, 0.125, 0.125],
                [-0.125, 0.125, 0.125],
                [0.125, 0.375, 0.125],
                [0.375, 0.375, 0.125],
                [0.625, 0.375, 0.125],
                [-0.125, 0.375, 0.125],
                [0.125, 0.625, 0.125],
                [0.375, 0.625, 0.125],
                [0.625, 0.625, 0.125],
                [-0.125, 0.625, 0.125],
                [0.125, -0.125, 0.125],
                [0.375, -0.125, 0.125],
                [0.625, -0.125, 0.125],
                [-0.125, -0.125, 0.125],
                [0.125, 0.125, 0.375],
                [0.375, 0.125, 0.375],
                [0.625, 0.125, 0.375],
                [-0.125, 0.125, 0.375],
                [0.125, 0.375, 0.375],
                [0.375, 0.375, 0.375],
                [0.625, 0.375, 0.375],
                [-0.125, 0.375, 0.375],
                [0.125, 0.625, 0.375],
                [0.375, 0.625, 0.375],
                [0.625, 0.625, 0.375],
                [-0.125, 0.625, 0.375],
                [0.125, -0.125, 0.375],
                [0.375, -0.125, 0.375],
                [0.625, -0.125, 0.375],
                [-0.125, -0.125, 0.375],
                [0.125, 0.125, 0.625],
                [0.375, 0.125, 0.625],
                [0.625, 0.125, 0.625],
                [-0.125, 0.125, 0.625],
                [0.125, 0.375, 0.625],
                [0.375, 0.375, 0.625],
                [0.625, 0.375, 0.625],
                [-0.125, 0.375, 0.625],
                [0.125, 0.625, 0.625],
                [0.375, 0.625, 0.625],
                [0.625, 0.625, 0.625],
                [-0.125, 0.625, 0.625],
                [0.125, -0.125, 0.625],
                [0.375, -0.125, 0.625],
                [0.625, -0.125, 0.625],
                [-0.125, -0.125, 0.625],
                [0.125, 0.125, -0.125],
                [0.375, 0.125, -0.125],
                [0.625, 0.125, -0.125],
                [-0.125, 0.125, -0.125],
                [0.125, 0.375, -0.125],
                [0.375, 0.375, -0.125],
                [0.625, 0.375, -0.125],
                [-0.125, 0.375, -0.125],
                [0.125, 0.625, -0.125],
                [0.375, 0.625, -0.125],
                [0.625, 0.625, -0.125],
                [-0.125, 0.625, -0.125],
                [0.125, -0.125, -0.125],
                [0.375, -0.125, -0.125],
                [0.625, -0.125, -0.125],
                [-0.125, -0.125, -0.125],
            ]
        end
        @testset "5×5×5" begin
            mesh = [5, 5, 5]
            nir, mapping, grid_address = get_ir_reciprocal_mesh(
                rutile,
                mesh,
                shift;
                is_time_reversal = true,
                symprec = 1e-5,
            )
            @test nir == length(unique(mapping)) == 18
            @test list_points(mapping, grid_address, mesh, shift, true) == [
                [0.1, 0.1, 0.1],
                [0.3, 0.1, 0.1],
                [0.5, 0.1, 0.1],
                [0.3, 0.3, 0.1],
                [0.5, 0.3, 0.1],
                [0.5, 0.5, 0.1],
                [0.1, 0.1, 0.3],
                [0.3, 0.1, 0.3],
                [0.5, 0.1, 0.3],
                [0.3, 0.3, 0.3],
                [0.5, 0.3, 0.3],
                [0.5, 0.5, 0.3],
                [0.1, 0.1, 0.5],
                [0.3, 0.1, 0.5],
                [0.5, 0.1, 0.5],
                [0.3, 0.3, 0.5],
                [0.5, 0.3, 0.5],
                [0.5, 0.5, 0.5],
            ]
            @test list_points(mapping, grid_address, mesh, shift, false) == [
                [0.1, 0.1, 0.1],
                [0.3, 0.1, 0.1],
                [0.5, 0.1, 0.1],
                [-0.3, 0.1, 0.1],
                [-0.1, 0.1, 0.1],
                [0.1, 0.3, 0.1],
                [0.3, 0.3, 0.1],
                [0.5, 0.3, 0.1],
                [-0.3, 0.3, 0.1],
                [-0.1, 0.3, 0.1],
                [0.1, 0.5, 0.1],
                [0.3, 0.5, 0.1],
                [0.5, 0.5, 0.1],
                [-0.3, 0.5, 0.1],
                [-0.1, 0.5, 0.1],
                [0.1, -0.3, 0.1],
                [0.3, -0.3, 0.1],
                [0.5, -0.3, 0.1],
                [-0.3, -0.3, 0.1],
                [-0.1, -0.3, 0.1],
                [0.1, -0.1, 0.1],
                [0.3, -0.1, 0.1],
                [0.5, -0.1, 0.1],
                [-0.3, -0.1, 0.1],
                [-0.1, -0.1, 0.1],
                [0.1, 0.1, 0.3],
                [0.3, 0.1, 0.3],
                [0.5, 0.1, 0.3],
                [-0.3, 0.1, 0.3],
                [-0.1, 0.1, 0.3],
                [0.1, 0.3, 0.3],
                [0.3, 0.3, 0.3],
                [0.5, 0.3, 0.3],
                [-0.3, 0.3, 0.3],
                [-0.1, 0.3, 0.3],
                [0.1, 0.5, 0.3],
                [0.3, 0.5, 0.3],
                [0.5, 0.5, 0.3],
                [-0.3, 0.5, 0.3],
                [-0.1, 0.5, 0.3],
                [0.1, -0.3, 0.3],
                [0.3, -0.3, 0.3],
                [0.5, -0.3, 0.3],
                [-0.3, -0.3, 0.3],
                [-0.1, -0.3, 0.3],
                [0.1, -0.1, 0.3],
                [0.3, -0.1, 0.3],
                [0.5, -0.1, 0.3],
                [-0.3, -0.1, 0.3],
                [-0.1, -0.1, 0.3],
                [0.1, 0.1, 0.5],
                [0.3, 0.1, 0.5],
                [0.5, 0.1, 0.5],
                [-0.3, 0.1, 0.5],
                [-0.1, 0.1, 0.5],
                [0.1, 0.3, 0.5],
                [0.3, 0.3, 0.5],
                [0.5, 0.3, 0.5],
                [-0.3, 0.3, 0.5],
                [-0.1, 0.3, 0.5],
                [0.1, 0.5, 0.5],
                [0.3, 0.5, 0.5],
                [0.5, 0.5, 0.5],
                [-0.3, 0.5, 0.5],
                [-0.1, 0.5, 0.5],
                [0.1, -0.3, 0.5],
                [0.3, -0.3, 0.5],
                [0.5, -0.3, 0.5],
                [-0.3, -0.3, 0.5],
                [-0.1, -0.3, 0.5],
                [0.1, -0.1, 0.5],
                [0.3, -0.1, 0.5],
                [0.5, -0.1, 0.5],
                [-0.3, -0.1, 0.5],
                [-0.1, -0.1, 0.5],
                [0.1, 0.1, -0.3],
                [0.3, 0.1, -0.3],
                [0.5, 0.1, -0.3],
                [-0.3, 0.1, -0.3],
                [-0.1, 0.1, -0.3],
                [0.1, 0.3, -0.3],
                [0.3, 0.3, -0.3],
                [0.5, 0.3, -0.3],
                [-0.3, 0.3, -0.3],
                [-0.1, 0.3, -0.3],
                [0.1, 0.5, -0.3],
                [0.3, 0.5, -0.3],
                [0.5, 0.5, -0.3],
                [-0.3, 0.5, -0.3],
                [-0.1, 0.5, -0.3],
                [0.1, -0.3, -0.3],
                [0.3, -0.3, -0.3],
                [0.5, -0.3, -0.3],
                [-0.3, -0.3, -0.3],
                [-0.1, -0.3, -0.3],
                [0.1, -0.1, -0.3],
                [0.3, -0.1, -0.3],
                [0.5, -0.1, -0.3],
                [-0.3, -0.1, -0.3],
                [-0.1, -0.1, -0.3],
                [0.1, 0.1, -0.1],
                [0.3, 0.1, -0.1],
                [0.5, 0.1, -0.1],
                [-0.3, 0.1, -0.1],
                [-0.1, 0.1, -0.1],
                [0.1, 0.3, -0.1],
                [0.3, 0.3, -0.1],
                [0.5, 0.3, -0.1],
                [-0.3, 0.3, -0.1],
                [-0.1, 0.3, -0.1],
                [0.1, 0.5, -0.1],
                [0.3, 0.5, -0.1],
                [0.5, 0.5, -0.1],
                [-0.3, 0.5, -0.1],
                [-0.1, 0.5, -0.1],
                [0.1, -0.3, -0.1],
                [0.3, -0.3, -0.1],
                [0.5, -0.3, -0.1],
                [-0.3, -0.3, -0.1],
                [-0.1, -0.3, -0.1],
                [0.1, -0.1, -0.1],
                [0.3, -0.1, -0.1],
                [0.5, -0.1, -0.1],
                [-0.3, -0.1, -0.1],
                [-0.1, -0.1, -0.1],
            ]
        end
    end
end

@testset "Test distorted silicon structure" begin
    lattice = [
        4.01 0 0
        0 4 0
        0 0 3.99
    ]
    positions = [
        0.001 0 0
        0 0.5 0.5
        0.5 0 0.5
        0.5 0.5 0
        0.25 0.25 0.251
        0.25 0.75 0.75
        0.75 0.25 0.75
        0.75 0.75 0.25
    ]
    types = [14, 14, 14, 14, 14, 14, 14, 14]
    silicon_dist = Cell(lattice, positions, types)
    mesh = [6, 6, 6]
    nir, mapping, grid_address =
        get_ir_reciprocal_mesh(silicon_dist, mesh; is_time_reversal = true, symprec = 1e-5)
    @test nir == length(unique(mapping)) == 112
    @test list_points(mapping, grid_address, mesh, [0, 0, 0], true) ≈ [
        [0.0, 0.0, 0.0],
        [0.16666667, 0.0, 0.0],
        [0.33333333, 0.0, 0.0],
        [0.5, 0.0, 0.0],
        [0.0, 0.16666667, 0.0],
        [0.16666667, 0.16666667, 0.0],
        [0.33333333, 0.16666667, 0.0],
        [0.5, 0.16666667, 0.0],
        [-0.33333333, 0.16666667, 0.0],
        [-0.16666667, 0.16666667, 0.0],
        [0.0, 0.33333333, 0.0],
        [0.16666667, 0.33333333, 0.0],
        [0.33333333, 0.33333333, 0.0],
        [0.5, 0.33333333, 0.0],
        [-0.33333333, 0.33333333, 0.0],
        [-0.16666667, 0.33333333, 0.0],
        [0.0, 0.5, 0.0],
        [0.16666667, 0.5, 0.0],
        [0.33333333, 0.5, 0.0],
        [0.5, 0.5, 0.0],
        [0.0, 0.0, 0.16666667],
        [0.16666667, 0.0, 0.16666667],
        [0.33333333, 0.0, 0.16666667],
        [0.5, 0.0, 0.16666667],
        [-0.33333333, 0.0, 0.16666667],
        [-0.16666667, 0.0, 0.16666667],
        [0.0, 0.16666667, 0.16666667],
        [0.16666667, 0.16666667, 0.16666667],
        [0.33333333, 0.16666667, 0.16666667],
        [0.5, 0.16666667, 0.16666667],
        [-0.33333333, 0.16666667, 0.16666667],
        [-0.16666667, 0.16666667, 0.16666667],
        [0.0, 0.33333333, 0.16666667],
        [0.16666667, 0.33333333, 0.16666667],
        [0.33333333, 0.33333333, 0.16666667],
        [0.5, 0.33333333, 0.16666667],
        [-0.33333333, 0.33333333, 0.16666667],
        [-0.16666667, 0.33333333, 0.16666667],
        [0.0, 0.5, 0.16666667],
        [0.16666667, 0.5, 0.16666667],
        [0.33333333, 0.5, 0.16666667],
        [0.5, 0.5, 0.16666667],
        [-0.33333333, 0.5, 0.16666667],
        [-0.16666667, 0.5, 0.16666667],
        [0.0, -0.33333333, 0.16666667],
        [0.16666667, -0.33333333, 0.16666667],
        [0.33333333, -0.33333333, 0.16666667],
        [0.5, -0.33333333, 0.16666667],
        [-0.33333333, -0.33333333, 0.16666667],
        [-0.16666667, -0.33333333, 0.16666667],
        [0.0, -0.16666667, 0.16666667],
        [0.16666667, -0.16666667, 0.16666667],
        [0.33333333, -0.16666667, 0.16666667],
        [0.5, -0.16666667, 0.16666667],
        [-0.33333333, -0.16666667, 0.16666667],
        [-0.16666667, -0.16666667, 0.16666667],
        [0.0, 0.0, 0.33333333],
        [0.16666667, 0.0, 0.33333333],
        [0.33333333, 0.0, 0.33333333],
        [0.5, 0.0, 0.33333333],
        [-0.33333333, 0.0, 0.33333333],
        [-0.16666667, 0.0, 0.33333333],
        [0.0, 0.16666667, 0.33333333],
        [0.16666667, 0.16666667, 0.33333333],
        [0.33333333, 0.16666667, 0.33333333],
        [0.5, 0.16666667, 0.33333333],
        [-0.33333333, 0.16666667, 0.33333333],
        [-0.16666667, 0.16666667, 0.33333333],
        [0.0, 0.33333333, 0.33333333],
        [0.16666667, 0.33333333, 0.33333333],
        [0.33333333, 0.33333333, 0.33333333],
        [0.5, 0.33333333, 0.33333333],
        [-0.33333333, 0.33333333, 0.33333333],
        [-0.16666667, 0.33333333, 0.33333333],
        [0.0, 0.5, 0.33333333],
        [0.16666667, 0.5, 0.33333333],
        [0.33333333, 0.5, 0.33333333],
        [0.5, 0.5, 0.33333333],
        [-0.33333333, 0.5, 0.33333333],
        [-0.16666667, 0.5, 0.33333333],
        [0.0, -0.33333333, 0.33333333],
        [0.16666667, -0.33333333, 0.33333333],
        [0.33333333, -0.33333333, 0.33333333],
        [0.5, -0.33333333, 0.33333333],
        [-0.33333333, -0.33333333, 0.33333333],
        [-0.16666667, -0.33333333, 0.33333333],
        [0.0, -0.16666667, 0.33333333],
        [0.16666667, -0.16666667, 0.33333333],
        [0.33333333, -0.16666667, 0.33333333],
        [0.5, -0.16666667, 0.33333333],
        [-0.33333333, -0.16666667, 0.33333333],
        [-0.16666667, -0.16666667, 0.33333333],
        [0.0, 0.0, 0.5],
        [0.16666667, 0.0, 0.5],
        [0.33333333, 0.0, 0.5],
        [0.5, 0.0, 0.5],
        [0.0, 0.16666667, 0.5],
        [0.16666667, 0.16666667, 0.5],
        [0.33333333, 0.16666667, 0.5],
        [0.5, 0.16666667, 0.5],
        [-0.33333333, 0.16666667, 0.5],
        [-0.16666667, 0.16666667, 0.5],
        [0.0, 0.33333333, 0.5],
        [0.16666667, 0.33333333, 0.5],
        [0.33333333, 0.33333333, 0.5],
        [0.5, 0.33333333, 0.5],
        [-0.33333333, 0.33333333, 0.5],
        [-0.16666667, 0.33333333, 0.5],
        [0.0, 0.5, 0.5],
        [0.16666667, 0.5, 0.5],
        [0.33333333, 0.5, 0.5],
        [0.5, 0.5, 0.5],
    ]
end

@testset "Test MgB₂ structure" begin
    lattice = [[3.07, 0.0, 0.0], [-1.535, 2.65869799, 0.0], [0.0, 0.0, 3.52]]
    positions = [[0, 0, 0], [0.33333333, 0.66666667, 0.5], [0.66666667, 0.33333333, 0.5]]
    types = [12, 5, 5]
    MgB₂ = Cell(lattice, positions, types)
    mesh = [7, 7, 7]
    nir, mapping, grid_address =
        get_ir_reciprocal_mesh(MgB₂, mesh; is_time_reversal = true, symprec = 1e-5)
    @test nir == length(unique(mapping)) == 32
    @test list_points(mapping, grid_address, mesh, [0, 0, 0], true) ≈ [
        [0.0, 0.0, 0.0],
        [0.14285714, 0.0, 0.0],
        [0.28571429, 0.0, 0.0],
        [0.42857143, 0.0, 0.0],
        [0.14285714, 0.14285714, 0.0],
        [0.28571429, 0.14285714, 0.0],
        [0.42857143, 0.14285714, 0.0],
        [0.28571429, 0.28571429, 0.0],
        [0.0, 0.0, 0.14285714],
        [0.14285714, 0.0, 0.14285714],
        [0.28571429, 0.0, 0.14285714],
        [0.42857143, 0.0, 0.14285714],
        [0.14285714, 0.14285714, 0.14285714],
        [0.28571429, 0.14285714, 0.14285714],
        [0.42857143, 0.14285714, 0.14285714],
        [0.28571429, 0.28571429, 0.14285714],
        [0.0, 0.0, 0.28571429],
        [0.14285714, 0.0, 0.28571429],
        [0.28571429, 0.0, 0.28571429],
        [0.42857143, 0.0, 0.28571429],
        [0.14285714, 0.14285714, 0.28571429],
        [0.28571429, 0.14285714, 0.28571429],
        [0.42857143, 0.14285714, 0.28571429],
        [0.28571429, 0.28571429, 0.28571429],
        [0.0, 0.0, 0.42857143],
        [0.14285714, 0.0, 0.42857143],
        [0.28571429, 0.0, 0.42857143],
        [0.42857143, 0.0, 0.42857143],
        [0.14285714, 0.14285714, 0.42857143],
        [0.28571429, 0.14285714, 0.42857143],
        [0.42857143, 0.14285714, 0.42857143],
        [0.28571429, 0.28571429, 0.42857143],
    ]
    @test length(list_points(mapping, grid_address, mesh, [0, 0, 0], false)) == 343
end


text:
random_seed = nothing

data_path = ""
data_prefix = "2d1ksample"

#Model Parameters
iterations = 99
hard_clustering = false
argmax_sample_stop = 10 #Change to hard assignment from soft at iterations - argmax_sample_stop
split_stop = 10 #Stop split/merge moves at  iterations - split_stop

total_dim = 2

α = 100000.0

initial_clusters = 1

use_dict_for_global = false

hyper_params = niw_hyperparams(1.0, [0, 0], 5.0, Matrix{Float64}(I, 2, 2) * 1.0)

#Saving specifics:
enable_saving = true
model_save_interval = 50
save_path = ""
overwrite_prec = false
save_file_prefix = "checkpoint_"

ext: jl


text:
using Test

include("common.jl")

@testset "TeaFiles.jl" begin
    @testset "Header" begin
        include("header_core.jl")
        include("header_read.jl")
        include("header_write.jl")
        include("header_create.jl")
    end
    @testset "Body" begin
        include("body_read.jl")
    end
    @testset "API" begin
        include("api.jl")
    end
    @testset "Example" begin
        include("example.jl")
    end
end


text:
using DifferentialEquations
using LinearAlgebra
using Flux
using Interpolations
# using SparseArrays

import Flux.Optimise: apply!

###### global constants ###########
const varnames = [:S,:I,:R,:D]
const parnames = [:β,:γ,:δ]

###### structures #################

mutable struct Baseline{T<:Real,D<:TimeType}
    time::StepRange{D}
    data::TimeArray{T,2,D}
    C::Matrix{T}
    σ::TimeArray{T,2,D}
    μ::TimeArray{T,2,D}
    u::TimeArray{T,2,D}
    g::TimeArray{T,2,D}
    p::TimeArray{T,2,D}

    function Baseline(time::StepRange{D},
                      data::TimeArray{T,2,D},
                      C::Matrix{T},
                      σ::TimeArray{T,2,D},
                      μ::TimeArray{T,2,D},
                      u::TimeArray{T,2,D},
                      g::TimeArray{T,2,D},
                      p::TimeArray{T,2,D},
                      ) where {T,D}

       @assert timestamp(data) == timestamp(σ) == timestamp(μ) ==
            timestamp(u) == timestamp(g) == timestamp(p) == time
            "input arguments timestamps are difform"

       size_C = (size(data,2),length(varnames))
       @assert size(C) == size_C "C must have size $size_C"

       @assert size(σ,2) == size(data,2) "σ must have $(size(data,2)) columns"
       @assert colnames(σ) == colnames(data) "columns of σ must have names $(colnames(data))"

       @assert size(μ,2) == length(parnames) "μ must have $(length(parnames)) columns"
       @assert colnames(μ) == parnames "columns of μ must have names $parnames"

       @assert size(u,2) == length(varnames) "u must have $(length(varnames)) columns"
       @assert colnames(u) == varnames "columns of u must have names $varnames"

       @assert size(p,2) == length(parnames) "u must have $(length(parnames)) columns"
       @assert colnames(p) == parnames "columns of p must have names $parnames"

       new{T,D}(time, data, C, σ, μ, u, g, p)
    end

    function Baseline(
                data::TimeArray{<:Real,2,<:TimeType},
                u₀::Dict{Symbol,<:Real},
                p::Union{Dict{Symbol,<:Real},TimeArray{<:Real,2,<:TimeType}};
                C=missing::Union{Missing,Matrix{<:Real}},
                start=timestamp(data)[1]::TimeType,
                stop=start+Day(30)::Period,
                step=Day(1)::Period,
                σ=missing::Union{Missing,TimeArray{<:Real,2,<:TimeType}},
                μ=missing::Union{Missing,TimeArray{<:Real,2,<:TimeType}},
                interptype=BSpline(Linear())::Union{Interpolations.InterpolationType,Tuple{Vararg{Interpolations.InterpolationType}}}
                )

        @assert length(u₀) == length(varnames) "u₀ must have length $(length(varnames))"

        if ismissing(C)
            if size(data,2) == length(varnames)
                _C = Diagonal(ones(length(varnames)))
            else
                error("C must be always specified when the number of observed and model parameters differ")
            end
        else
            if size(C) != (size(data,2),length(varnames))
                error("datamap must have size $(size(data,2))×$(length(varnames))")
            else
                _C = float(C)
            end
        end

        @assert start >= timestamp(data)[1] "start date must be after the last date available in data"
        @assert start < timestamp(data)[end] "start date must be before the last date available in data"
        @assert stop > start "stop date must be after the start date"
        @assert step > Day(0) "step must be a positive time period"

        @assert size(σ,2) == size(data,2) "σ must have $(size(data,2)) columns"
        @assert colnames(σ) == colnames(data) "columns of σ must have names $(colnames(data))"

        @assert size(μ,2) == length(parnames) "μ must have $(length(parnames)) columns"
        @assert colnames(μ) == parnames "columns of μ must have names $parnames"

        if interptype isa Tuple
            @assert length(interptype) == size(data,2) "interptype must be a tuple with $(size(data,2)) elements"
        else
            interptype = ntuple(i -> interptype, size(data,2))
        end

        # time
        if step < Day(1)
            t = DateTime(start):step:DateTime(stop)
        else
            t = start:step:stop
        end

        # Copy data to working data array ux
        _data = TimeArray(t,zeros(length(t),size(data,2)),colnames(data),Dict{String,Any}())
        for j=1:size(data,2)
            values(_data)[:,j] = _interpolation(timerange(timestamp(data)),
                                    values(data)[:,j], t + Millisecond(step)/2,
                                    #BSpline(Constant()), Flat())
                                    interptype[j], Flat())
        end
        data_t = convert.(eltype(t), timestamp(data))
        meta(_data)["lastdate"] = convert(eltype(t), timestamp(data)[end])

        # Initialize forcing window
        _σ = TimeArray(t,zeros(size(_data)),colnames(data))
        if ismissing(σ)
            σ = TimeArray(timestamp(data),ones(size(data,2)),varnames)
        end
        for j=1:size(σ,2)
            values(_σ)[:,j] = _interpolation(timerange(timestamp(σ)),
                                values(σ)[:,j], t, BSpline(Linear()), 0.0)
        end

        # Initialize control window
        _μ = TimeArray(t,zeros(length(t),length(parnames)),parnames)
        if ismissing(μ)
            μ = TimeArray(timestamp(data),ones(size(data,1),length(parnames)),parnames)
        end
        for j=1:size(μ,2)
            values(_μ)[:,j] = _interpolation(timerange(timestamp(μ)),
                                values(μ)[:,j], t, BSpline(Linear()), 0.0)
        end

        # Initialize model parameters
        _p = TimeArray(t,zeros(length(t),length(parnames)),parnames)
        if p isa Dict
            for (k,v) in enumerate(parnames)
                haskey(p,v) || error("p does not have key $v")
                values(_p)[:,k] .= p[v]
            end
        else
            for j=1:size(p,2)
                values(_p)[:,j] = _interpolation(timerange(timestamp(p)),
                                    values(p)[:,j], t, BSpline(Linear()), Flat())
            end
        end

        # Initialize baseline data array
        _u₀ = OrderedDict{Symbol,Float64}()
        for v in varnames
            haskey(u₀,v) || error("u₀ does not have key $v")
            _u₀[v] = u₀[v]
        end
        _u = forward(t, values(_p), _u₀)

        # Initialize forcing
        _g = forcing(_data,_C,_σ,_u)

        new{Float64,eltype(t)}(t, _data, _C, _σ, _μ, _u, _g, _p)
    end
end

struct BaselineInterpolation{T<:Real}
    data::AbstractInterpolation{T,2}
    u::AbstractInterpolation{T,2}
    g::AbstractInterpolation{T,2}
    p::AbstractInterpolation{T,2}

    function BaselineInterpolation(base::Baseline{T}, trange::AbstractRange{<:Real}) where {T}
        interptype = (BSpline(Linear()), NoInterp())
        data = scale(interpolate(values(base.data), interptype), trange, 1:size(base.data,2))
        u = scale(interpolate(values(base.u), interptype), trange, 1:size(base.u,2))
        g = scale(interpolate(values(base.g), interptype), trange, 1:size(base.g,2))
        p = scale(interpolate(values(base.p), interptype), trange, 1:size(base.p,2))
        new{T}(data, u, g, p)
    end
end

mutable struct DA{T<:Real,D<:TimeType}
    time::StepRange{D}
    data::TimeArray{T,2,D}
    C::Matrix{T}
    σ::TimeArray{T,2,D}
    μ::TimeArray{T,2,D}
    u::TimeArray{T,2,D}
    g::TimeArray{T,2,D}
    p::TimeArray{T,2,D}
    v::TimeArray{T,2,D}
    δp::TimeArray{T,2,D}

    function DA(base::Baseline{T,D}) where {T,D}
        v = backward(base)
        δp = gradient(base,v)
        new{T,D}(base.time, base.data, base.C, base.σ, base.μ,
                    base.u, base.g,base.p, v, δp)
    end

    function DA(da::DA, p::TimeArray{<:Real,2})
        u₀ = OrderedDict(colnames(da.u) .=> values(da.u)[1,:])
        u = forward(da.time, values(p), u₀)
        g = forcing(da.data, da.C ,da.σ ,u)
        base = Baseline(da.time, da.data, da.C, da.σ, da.μ, u, g, p)
        DA(base)
    end
end

function baseline(da::DA)
    Baseline(da.time, da.data, da.C, da.σ, da.μ, da.u, da.g, da.p)
end

function update(da::DA, α::Real)
    update(da, α.*da.δp)
end

function update(da::DA, δp::TimeArray{<:Real,2})
    update(da, values(δp))
end

function update(da::DA, δp::Array{<:Real,2})
    p = max.(da.p .- δp, 0.0)
    DA(da, p)
end

function residual(base::Baseline; relative=false::Bool)
    J = 0.5 * norm(values(base.g), 2)
    if relative
        J = J/datanorm(base)
    end
    J
end

function datanorm(base::Baseline)
    0.5 * norm(values(base.data) .*values(base.σ) * base.C, 2)
end

function extend_solution!(da::DA)
    lastdate = meta(da.data)["lastdate"] - Day(1)
    idxs = findfirst(lastdate .== timestamp(da.data))
    values(da.p)[idxs+1:end,:] .= values(da.p)[idxs:idxs,:]
end

function timerange(d::AbstractArray{<:TimeType})
    drange = range(d[1], d[end], step = d[2] - d[1])
    @assert d == drange "d must be equispaced in time"
    drange
end

function timestep(trange::StepRange{<:TimeType,<:Period})
    Dates.toms(step(trange))/Dates.toms(Day(1))
end

import Interpolations.interpolate
interpolate(range::AbstractRange, vs::AbstractVector, interptype::BSpline) = scale(interpolate(vs, interptype), range)

function _interpolation(t::AbstractRange{<:TimeType}, u::AbstractVector{<:Real},
                    tq::AbstractVector{<:TimeType}, interptype, extrapscheme)
        @assert step(t) > Millisecond(0) "step of t must be positive"
        x = datetime2unix(DateTime(t.start)):Second(t.step).value:datetime2unix(DateTime(t.stop))
        xq = datetime2unix.(DateTime.(tq))
        itp = interpolate(x, u, interptype)
        etp = extrapolate(itp, extrapscheme)
        etp.(xq)
end

"Classic Epidemic Model  Hethcote (2000), added deaths"
function sir!(du::Vector{Float64}, u::Vector{Float64},
              p::AbstractInterpolation{Float64,2}, t::Float64)
    β,γ,δ = p(t,1:3)
    S,I,R,D = u[:]
    N = sum(u)
    du[1]=-β*I*S/N
    du[2]= β*I*S/N - (γ+δ)*I
    du[3]= γ*I
    du[4]= δ*I
    return du
end

function dfdu(u::Vector{Float64}, p::Vector{Float64})
    S,I,R,D = u[:]
    β,γ,δ = p[:]

    N=S+I+R+D
    A=zeros(4,4)
    A[1,1]=-β*I*(N+S)/N^2
    A[1,2]=-β*S*(N+I)/N^2
    A[1,3]=+β*I*S/N^2
    A[1,4]=+β*I*S/N^2

    A[2,1]= β*I*(N+S)/N^2
    A[2,2]= β*S*(N+I)/N^2 - (γ+δ)
    A[2,3]=-β*I*S/N^2
    A[2,4]=-β*I*S/N^2

    A[3,2]=γ
    A[4,2]=δ
    return A
end

function sir_adj(v::Vector{Float64}, base::BaselineInterpolation, t::Float64)
    u = base.u(t,1:4)
    p = base.p(t,1:3)
    g = base.g(t,1:4)
    A = dfdu(u,p)
    dv = -A'*v - g
    return dv
end

function dfda(u::Matrix{Float64})
    S=u[:,1]
    I=u[:,2]
    R=u[:,3]
    D=u[:,4]

    N=S+I+R+D
    f=zeros(length(S),4,3)
    f[:,1,1]=-I.*S./N
    f[:,1,2].= 0.0
    f[:,1,3].= 0.0

    f[:,2,1]= I.*S./N
    f[:,2,2]=-I
    f[:,2,3]=-I

    f[:,3,1].=0.0
    f[:,3,2]= I
    f[:,3,3].=0.0

    f[:,4,1].=0.0
    f[:,4,2].=0.0
    f[:,4,3]= I

    return f
end

function forward(t::StepRange{<:TimeType}, p::Array{Float64,2}, u₀::OrderedDict{Symbol,Float64})
    Δt = timestep(t)
    tspan = (0.0, (size(p,1) - 1)*Δt)
    trange = range(tspan..., length=size(p,1))
    _u₀ = collect(values(u₀))
    interptype = (BSpline(Linear()), NoInterp())
    itp = scale(interpolate(p, interptype), trange, 1:size(p,2))
    problem = ODEProblem(sir!, _u₀, tspan, itp)
    solution = solve(problem, RK4(), dt=1.0, adaptive=false)
    u = collect(solution(trange)')
    TimeArray(t, u, collect(keys(u₀)))
end

function backward(base::Baseline)
    Δt = timestep(base.time)
    tspan = ((size(base.u,1) - 1)*Δt, 0.0)
    trange = reverse(range(tspan...,length=size(base.u,1)))
    v₀ = zeros(size(base.u,2))
    base_itp = BaselineInterpolation(base, trange)
    adj = ODEProblem(sir_adj, v₀, tspan, base_itp)
    adj_sol = solve(adj, tstops=trange)
    adj_sol = solve(adj, RK4(), dt=1.0, adaptive=false)
    v = collect(adj_sol(trange)')
    TimeArray(base.time, v, colnames(base.u))
end

function gradient(base::Baseline, v::TimeArray)
    u = values(base.u)
    p = values(base.p)
    f = dfda(u)
    δp = zeros(size(p))
    for i=1:size(u,2)
        δp = δp + values(v)[:,i] .* f[:,i,:]
    end
    δp ./= maximum(abs.(δp))
    TimeArray(timestamp(base.p), δp, colnames(base.p))
end

function forcing(data::TimeArray{Float64,2}, C::Matrix{Float64},
                 σ::TimeArray{Float64,2}, u::TimeArray{Float64,2})
    g = (values(u) * C' - values(data)) .* values(σ) * C
    TimeArray(timestamp(u), g, colnames(u))
end

function apply!(opt, da::DA)
    δp = values(da.δp) .* values(da.μ)
    Δp = apply!(opt, values(da.p), δp)
    update(da, Δp)
end


text:
"""   `calc_basis_scores_scalpels(rvs, ccfs; σ_rvs, num_basis )`
Compute the CCF basis functions and scores for a Scalpels reconstruction of CCFs
Inputs:
- rvs:  vector of estimated radial velocities
- ccfs: 2d array of CCFS of size (num_velocities, num_spectra)
Optional Inputs:
- σ_rvs:  vector of measurement uncertainties for estimated radial velocities (default: ones)
- num_basis:  number of basis vectors to use for SVD reconstruction of CCFs
- sort_by_responce:  set true to sort basis vectors by RV responce (default: true)
Output:
(scores, basis): a NamedTuple

Notes:
- Currently Scalpels weights all velocity pixels equally and uses the σ_rvs to weight each observation.
- First element of output starts with RMS with zero basis vectors (i.e., the input RVs)
"""
function calc_basis_scores_scalpels(rvs::AbstractVector{T1}, ccfs::AbstractArray{T2,2}
                ; σ_rvs::AbstractVector{T3} = ones(length(rvs)),
                num_basis::Integer = 3, assume_centered::Bool = false, sort_by_responce::Bool = true ) where { T1<:Real, T2<:Real, T3<:Real }
    @assert length(rvs) == length(σ_rvs)
    @assert length(rvs) == size(ccfs,2)
    @assert 1 <= num_basis < length(rvs)
    if assume_centered
        rvs_centered = rvs
    else
        mean_rv = mean(rvs, weights(1.0 ./ σ_rvs.^2 ))
        rvs_centered = rvs .- mean_rv
    end

    acfs = autocor(ccfs,0:size(ccfs,1)-1)
    if all(σ_rvs.==first(σ_rvs))
        acfs_minus_mean = acfs .- mean(acfs,dims=2)
    else
        acfs_minus_mean = acfs .- mean(acfs,weights(1.0 ./ σ_rvs.^2),dims=2)
    end
    #Δv_grid = convert(Float64,v_grid.step).*(0:size(acfs,1)-1)
    svd_acfs = svd(acfs_minus_mean')
    #return (basis=svd_acfs.U, scores=svd_acfs.S)  For ACF's SVD eigenvectors and eigenvalues

    alpha = svd_acfs.U'*rvs_centered
    if sort_by_responce
        idx = sortperm(abs.(alpha),rev=true)
        #println("# Ordered by alphas: ", idx[1:num_basis])
    else
        # To throw away sorting by projection onto RVs
        idx = 1:length(alpha)
    end
    U_keep = view(svd_acfs.U,:,idx[1:num_basis])
    P_keep = view(svd_acfs.Vt',:,idx[1:num_basis])
    return (scores=U_keep, basis=P_keep) # , responce=alpha[idx[1:num_basis]])

    #U_keep = view(svd_acfs.U,:,1:num_basis)
    #return (basis=U_keep, scores = view(alpha,1:num_basis) )
    #=  To clean RVs
    Δrv_shape = U_keep*U_keep'*rvs_centered
    rvs_clean = rvs .- Δrv_shape
    return rvs_clean
    =#
end

"""   `clean_rvs_scalpels(rvs, ccfs; σ_rvs, num_basis )`
Inputs:
- rvs:  vector of estimated radial velocities
- ccfs: 2d array of CCFS of size (num_velocities, num_spectra)
Optional Inputs:
- σ_rvs:  vector of measurement uncertainties for estimated radial velocities (default: ones)
- num_basis:  number of basis vectors to use for SVD reconstruction of CCFs
- sort_by_responce:  set true to sort basis vectors by RV responce (default: true)
Output:
rvs_clean: vector of estimated RVs after cleaning by scalpels

Notes:
- Currently Scalpels weights all observations equally and doesn't use the σ_rvs.
- First element of output starts with RMS with zero basis vectors (i.e., the input RVs)
"""
function clean_rvs_scalpels(rvs::AbstractVector{T1}, ccfs::AbstractArray{T2,2}
                ; σ_rvs::AbstractVector{T3} = ones(length(rvs)),
                num_basis::Integer = 3, sort_by_responce::Bool = true ) where { T1<:Real, T2<:Real, T3<:Real }
    @assert length(rvs) == length(σ_rvs)
    @assert length(rvs) == size(ccfs,2)
    if num_basis == 0   return rvs   end
    @assert 0 <= num_basis < length(rvs)
    mean_rv = mean(rvs, weights(1.0 ./ σ_rvs.^2 ))
    rvs_centered = rvs .- mean_rv

    #=
    acfs = autocor(ccfs,0:size(ccfs,1)-1)
    acfs_minus_mean = acfs .- mean(acfs,dims=2)
    #Δv_grid = convert(Float64,v_grid.step).*(0:size(acfs,1)-1)
    svd_acfs = svd(acfs_minus_mean')
    alpha = svd_acfs.U'*rvs_centered

    idx = sortperm(abs.(alpha),rev=true)
    U_keep = view(svd_acfs.U,:,idx[1:num_basis])

    #(U, scores ) = calc_basis_scores_scalpels(rvs_centered, ccfs, σ_rvs=σ_rvs, num_basis=num_basis, assume_centered=true)
    idx = sortperm(abs.(alpha),rev=true)
    alpha = U'*rvs_centered
    U_keep = view(U,:,idx[1:num_basis])
    =#

    (U_keep, ) = calc_basis_scores_scalpels(rvs_centered, ccfs, σ_rvs=σ_rvs, num_basis=num_basis, assume_centered=true, sort_by_responce=sort_by_responce)
    Δrv_shape = U_keep*U_keep'*rvs_centered
    rvs_clean = rvs .- Δrv_shape
    return rvs_clean
end

"""   `rms_clean_rvs_vs_num_basis_scalpels(rvs, ccfs; σ_rvs, max_num_basis )`
Compute RMS of estimated RVs after cleaning raw RVS with Scalpels algorithm (based on CCFs)

Inputs:
- rvs:  vector of estimated radial velocities
- ccfs: 2d array of CCFS of size (num_velocities, num_spectra)
Optional Inputs:
- σ_rvs:  vector of measurement uncertainties for estimated radial velocities (default: ones)
- max_num_basis:  maximum number of basis vectors to use for SVD reconstruction of CCFs
- sort_by_responce:  set true to sort basis vectors by RV responce (default: true)
Output:
rms_scalpels: vector of RMS estimated RVs after cleaning by scalpels as a function of the number of basis vectors

Notes:
- Currently Scalpels weights all observations equally and doesn't use the σ_rvs.
- First element of output starts with RMS with zero basis vectors (i.e., the input RVs)
"""
function rms_clean_rvs_vs_num_basis_scalpels(rvs::AbstractVector{T1}, ccfs::AbstractArray{T2,2}
                ; σ_rvs::AbstractVector{T3} = ones(length(rvs)),
                  max_num_basis::Integer = min(length(rvs)-1,default_max_num_basis),
                  sort_by_responce::Bool = true ) where { T1<:Real, T2<:Real, T3<:Real }
    rms_scalpels = map(b->std(clean_rvs_scalpels(rvs,ccfs,num_basis=b, sort_by_responce=sort_by_responce)), 0:max_num_basis)
end

"""   `calc_clean_rvs_scores_basis_scalpels(rvs, ccfs; σ_rvs, num_basis, sort_by_responce )`
Computes cleaned rvs as well as the CCF basis functions and scores for a Scalpels reconstruction of CCFs
Inputs:
- rvs:  vector of estimated radial velocities
- ccfs: 2d array of CCFS of size (num_velocities, num_spectra)
Optional Inputs:
- σ_rvs:  vector of measurement uncertainties for estimated radial velocities (default: ones)
- num_basis:  number of basis vectors to use for SVD reconstruction of CCFs
- sort_by_responce:  set true to sort basis vectors by RV responce (default: true)
Output:
(rvs, scores, basis): a NamedTuple

Notes:
- Currently Scalpels weights all velocity pixels equally and uses the σ_rvs to weight each observation.
- First element of output starts with RMS with zero basis vectors (i.e., the input RVs)
"""
function calc_clean_rvs_scores_basis_scalpels(rvs::AbstractVector{T1}, ccfs::AbstractArray{T2,2}
                ; σ_rvs::AbstractVector{T3} = ones(length(rvs)),
                num_basis::Integer = 3,
                sort_by_responce::Bool = true ) where { T1<:Real, T2<:Real, T3<:Real }
    @assert length(rvs) == length(σ_rvs)
    @assert length(rvs) == size(ccfs,2)
    if num_basis == 0   return rvs   end
    @assert 0 <= num_basis < length(rvs)
    mean_rv = mean(rvs, weights(1.0 ./ σ_rvs.^2 ))
    rvs_centered = rvs .- mean_rv

    (U_keep, basis ) = calc_basis_scores_scalpels(rvs_centered, ccfs, σ_rvs=σ_rvs, num_basis=num_basis, assume_centered=true, sort_by_responce=sort_by_responce)
    Δrv_shape = U_keep*U_keep'*rvs_centered
    rvs_clean = rvs .- Δrv_shape

    return (rvs=rvs_clean, scores=U_keep, basis=basis )
end


text:
@safetestset "isip" begin
    using Validator

    # IPV4
    @test Validator.isip("127.0.0.1", 4) === true
    @test Validator.isip("0.0.0.0", 4) === true
    @test Validator.isip("255.255.255.255", 4) === true
    @test Validator.isip("1.2.3.4", 4) === true
    @test Validator.isip("255.0.0.1", 4) === true
    @test Validator.isip("0.0.1.1", 4) === true

    @test Validator.isip("::1", 4) === false
    @test Validator.isip("2001:db8:0000:1:1:1:1:1", 4) === false
    @test Validator.isip("::ffff:127.0.0.1", 4) === false
    @test Validator.isip("137.132.10.01", 4) === false
    @test Validator.isip("0.256.0.256", 4) === false
    @test Validator.isip("255.256.255.256", 4) === false

    # IPV6
    @test Validator.isip("::1", 6) === true
    @test Validator.isip("2001:db8:0000:1:1:1:1:1", 6) === true
    @test Validator.isip("::ffff:127.0.0.1", 6) === true
    @test Validator.isip("fe80::1234%1", 6) === true
    @test Validator.isip("ff08::9abc%10", 6) === true
    @test Validator.isip("ff08::9abc%interface10", 6) === true
    @test Validator.isip("ff02::5678%pvc1.3", 6) === true

    @test Validator.isip("127.0.0.1", 6) === false
    @test Validator.isip("0.0.0.0'", 6) === false
    @test Validator.isip("255.255.255.255", 6) === false
    @test Validator.isip("1.2.3.4", 6) === false
    @test Validator.isip("::ffff:287.0.0.1", 6) === false
    @test Validator.isip("%", 6) === false
    @test Validator.isip("fe80::1234%", 6) === false
    @test Validator.isip("fe80::1234%1%3%4", 6) === false
    @test Validator.isip("fe80%fe80%", 6) === false

    # No type specified
    @test Validator.isip("127.0.0.1") === true
    @test Validator.isip("0.0.0.0") === true
    @test Validator.isip("255.255.255.255") === true
    @test Validator.isip("1.2.3.4") === true
    @test Validator.isip("::1") === true
    @test Validator.isip("2001:db8:0000:1:1:1:1:1") === true
    @test Validator.isip("2001:db8:3:4::192.0.2.33") === true
    @test Validator.isip("2001:41d0:2:a141::1") === true
    @test Validator.isip("::ffff:127.0.0.1") === true
    @test Validator.isip("::0000") === true
    @test Validator.isip("0000::") === true
    @test Validator.isip("1::") === true
    @test Validator.isip("1111:1:1:1:1:1:1:1") === true
    @test Validator.isip("fe80::a6db:30ff:fe98:e946") === true
    @test Validator.isip("::") === true
    @test Validator.isip("::8") === true
    @test Validator.isip("::ffff:127.0.0.1") === true
    @test Validator.isip("::ffff:255.255.255.255") === true
    @test Validator.isip("::ffff:0:255.255.255.255") === true
    @test Validator.isip("::2:3:4:5:6:7:8") === true
    @test Validator.isip("::255.255.255.255") === true
    @test Validator.isip("0:0:0:0:0:ffff:127.0.0.1") === true
    @test Validator.isip("1:2:3:4:5:6:7::") === true
    @test Validator.isip("1:2:3:4:5:6::8") === true
    @test Validator.isip("1::7:8") === true
    @test Validator.isip("1:2:3:4:5::7:8") === true
    @test Validator.isip("1:2:3:4:5::8") === true
    @test Validator.isip("1::6:7:8") === true
    @test Validator.isip("1:2:3:4::6:7:8") === true
    @test Validator.isip("1:2:3:4::8") === true
    @test Validator.isip("1::5:6:7:8") === true
    @test Validator.isip("1:2:3::5:6:7:8") === true
    @test Validator.isip("1:2:3::8") === true
    @test Validator.isip("1::4:5:6:7:8") === true
    @test Validator.isip("1:2::4:5:6:7:8") === true
    @test Validator.isip("1:2::8") === true
    @test Validator.isip("1::3:4:5:6:7:8") === true
    @test Validator.isip("1::8") === true
    @test Validator.isip("fe80::7:8%eth0") === true
    @test Validator.isip("fe80::7:8%1") === true
    @test Validator.isip("64:ff9b::192.0.2.33") === true
    @test Validator.isip("0:0:0:0:0:0:10.0.0.1") === true

    @test Validator.isip("abc") === false
    @test Validator.isip("256.0.0.0") === false
    @test Validator.isip("0.0.0.256") === false
    @test Validator.isip("26.0.0.256") === false
    @test Validator.isip("0200.200.200.200") === false
    @test Validator.isip("200.0200.200.200") === false
    @test Validator.isip("200.200.0200.200") === false
    @test Validator.isip("200.200.200.0200") === false
    @test Validator.isip("::banana") === false
    @test Validator.isip("banana::") === false
    @test Validator.isip("::1banana") === false
    @test Validator.isip("::1::") === false
    @test Validator.isip("1:") === false
    @test Validator.isip(":1") === false
    @test Validator.isip(":1:1:1::2") === false
    @test Validator.isip("1:1:1:1:1:1:1:1:1:1:1:1:1:1:1:1") === false
    @test Validator.isip("::11111") === false
    @test Validator.isip("11111:1:1:1:1:1:1:1") === false
    @test Validator.isip("2001:db8:0000:1:1:1:1::1") === false
    @test Validator.isip("0:0:0:0:0:0:ffff:127.0.0.1") === false
    @test Validator.isip("0:0:0:0:ffff:127.0.0.1") === false
end


text:
function run(s)
    # Your code here
    joltages = readlines(IOBuffer(s))
    joltages = [parse(Int, j) for j in joltages]
    sort!(joltages)
    pushfirst!(joltages, 0)
    diffs = joltages[2:end] - joltages[1:end - 1]
    return sum(diffs .== 1) * (sum(diffs .== 3) + 1)
end

#########################################

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time * 1000)")
    println(res)
end

main()


text:
function construct1S(Z,center::Point3D{T},cell::Array{Point3D{T}},dims,supercells) where T
  n_a,n_b,n_c = dims
  origin = -(div(supercells-1,2)+1)*sum(cell)
  a_vec,b_vec,c_vec = cell.*supercells
  points = Array{WfPoint3D{T},3}(n_a,n_b,n_c)
  R,Y = zero(T),zero(T)
  for (ia,a) in enumerate(linspace(0,1,n_a)),(ib,b) in enumerate(linspace(0,1,n_b)),(ic,c) in enumerate(linspace(0,1,n_c))
    point = origin+(a*a_vec+b*b_vec+c*c_vec)
    R = 2*Z^1.5 * exp(-Z*norm(point-center))
    Y = sqrt(1/4pi)
    points[ia,ib,ic] = WfPoint3D{T}(R*Y+0.0im,origin+point)
  end
  return Wfc3D{T}(points,cell,PhysAtom(T))
end

function construct2S(Z,center::Point3D{T},cell::Array{Point3D{T}},dims,supercells) where T
  n_a,n_b,n_c = dims
  origin = -(div(supercells-1,2)+1)*sum(cell)
  a_vec,b_vec,c_vec = cell.*supercells
  points = Array{WfPoint3D{T},3}(n_a,n_b,n_c)
  R,Y = zero(T),zero(T)
  for (ia,a) in enumerate(linspace(0,1,n_a)),(ib,b) in enumerate(linspace(0,1,n_b)),(ic,c) in enumerate(linspace(0,1,n_c))
    point = origin+(a*a_vec+b*b_vec+c*c_vec)
    rho = 2*Z*norm(point-center)/2
    R = 1/(2*sqrt(2)) *(2 - rho) * Z^1.5 *exp(-rho/2)
    Y = sqrt(1/4pi)
    points[ia,ib,ic] = WfPoint3D{T}(R*Y+0.0im,origin+point)
  end
  return Wfc3D{T}(points,cell,PhysAtom(T))
end

function construct3S(Z,center::Point3D{T},cell::Array{Point3D{T}},dims,supercells) where T
  n_a,n_b,n_c = dims
  origin = -(div(supercells-1,2)+1)*sum(cell)
  a_vec,b_vec,c_vec = cell.*supercells
  points = Array{WfPoint3D{T},3}(n_a,n_b,n_c)
  R,Y = zero(T),zero(T)
  for (ia,a) in enumerate(linspace(0,1,n_a)),(ib,b) in enumerate(linspace(0,1,n_b)),(ic,c) in enumerate(linspace(0,1,n_c))
    point = origin+(a*a_vec+b*b_vec+c*c_vec)
    rho = 2*Z*norm(point-center)/3
    R = 1/(9*sqrt(3)) *(6-6rho+rho^2) * Z^1.5 *exp(-rho/2)
    Y = sqrt(1/4pi)
    points[ia,ib,ic] = WfPoint3D{T}(R*Y+0.0im,origin+point)
  end
  return Wfc3D{T}(points,cell,PhysAtom(T))
end

function construct4S(Z,center::Point3D{T},cell::Array{Point3D{T}},dims,supercells) where T
  n_a,n_b,n_c = dims
  origin = -(div(supercells-1,2)+1)*sum(cell)
  a_vec,b_vec,c_vec = cell.*supercells
  points = Array{WfPoint3D{T},3}(n_a,n_b,n_c)
  R,Y = zero(T),zero(T)
  for (ia,a) in enumerate(linspace(0,1,n_a)),(ib,b) in enumerate(linspace(0,1,n_b)),(ic,c) in enumerate(linspace(0,1,n_c))
    point = origin+(a*a_vec+b*b_vec+c*c_vec)
    rho = 2*Z*norm(point-center)/4
    R = 1/(96) *(24-36rho+12rho^2-rho^3) * Z^1.5 *exp(-rho/2)
    Y = sqrt(1/4pi)
    points[ia,ib,ic] = WfPoint3D{T}(R*Y+0.0im,origin+point)
  end
  return Wfc3D{T}(points,cell,PhysAtom(T))
end

function construct5S(Z,center::Point3D{T},cell::Array{Point3D{T}},dims,supercells) where T
  n_a,n_b,n_c = dims
  origin = -(div(supercells-1,2)+1)*sum(cell)
  a_vec,b_vec,c_vec = cell.*supercells
  points = Array{WfPoint3D{T},3}(n_a,n_b,n_c)
  R,Y = zero(T),zero(T)
  for (ia,a) in enumerate(linspace(0,1,n_a)),(ib,b) in enumerate(linspace(0,1,n_b)),(ic,c) in enumerate(linspace(0,1,n_c))
    point = origin+(a*a_vec+b*b_vec+c*c_vec)
    rho = 2*Z*norm(point-center)/5
    R = 1/(300sqrt(5)) *(120-240rho+120rho^2-20rho^3+rho^4) * Z^1.5 *exp(-rho/2)
    Y = sqrt(1/4pi)
    points[ia,ib,ic] = WfPoint3D{T}(R*Y+0.0im,origin+point)
  end
  return Wfc3D{T}(points,cell,PhysAtom(T))
end

function construct2P(direction,Z,center::Point3D{T},cell::Array{Point3D{T}},dims,supercells) where T
  n_a,n_b,n_c = dims
  origin = -(div(supercells-1,2)+1)*sum(cell)
  a_vec,b_vec,c_vec = cell.*supercells
  points = Array{WfPoint3D{T},3}(n_a,n_b,n_c)
  R,Y = zero(T),zero(T)
  for (ia,a) in enumerate(linspace(0,1,n_a)),(ib,b) in enumerate(linspace(0,1,n_b)),(ic,c) in enumerate(linspace(0,1,n_c))
    point = origin+(a*a_vec+b*b_vec+c*c_vec)
    rho = 2*Z*norm(point-center)/2
    R = 1/2sqrt(6) * rho * Z^1.5 *exp(-rho/2)
    Y = sqrt(3)*getfield(point-center,direction)/norm(point-center)*sqrt(1/4pi)
    points[ia,ib,ic] = WfPoint3D{T}(R*Y+0.0im,origin+point)
  end
  return Wfc3D{T}(points,cell,PhysAtom(T))
end

function construct3P(direction,Z,center::Point3D{T},cell::Array{Point3D{T}},dims,supercells) where T
  n_a,n_b,n_c = dims
  origin = -(div(supercells-1,2)+1)*sum(cell)
  a_vec,b_vec,c_vec = cell.*supercells
  points = Array{WfPoint3D{T},3}(n_a,n_b,n_c)
  R,Y = zero(T),zero(T)
  for (ia,a) in enumerate(linspace(0,1,n_a)),(ib,b) in enumerate(linspace(0,1,n_b)),(ic,c) in enumerate(linspace(0,1,n_c))
    point = origin+(a*a_vec+b*b_vec+c*c_vec)
    rho = 2*Z*norm(point-center)/3
    R = 1/9sqrt(6) * rho*(4-rho) * Z^1.5 *exp(-rho/2)
    Y = sqrt(3)*getfield(point-center,direction)/norm(point-center)*sqrt(1/4pi)
    points[ia,ib,ic] = WfPoint3D{T}(R*Y+0.0im,origin+point)
  end
  return Wfc3D{T}(points,cell,PhysAtom(T))
end

function construct4P(direction,Z,center::Point3D{T},cell::Array{Point3D{T}},dims,supercells) where T
  n_a,n_b,n_c = dims
  origin = -(div(supercells-1,2)+1)*sum(cell)
  a_vec,b_vec,c_vec = cell.*supercells
  points = Array{WfPoint3D{T},3}(n_a,n_b,n_c)
  R,Y = zero(T),zero(T)
  for (ia,a) in enumerate(linspace(0,1,n_a)),(ib,b) in enumerate(linspace(0,1,n_b)),(ic,c) in enumerate(linspace(0,1,n_c))
    point = origin+(a*a_vec+b*b_vec+c*c_vec)
    rho = 2*Z*norm(point-center)/4
    R = 1/(32*sqrt(15)) * rho*(20-10rho+rho^2) * Z^1.5 *exp(-rho/2)
    Y = sqrt(3)*getfield(point-center,direction)/norm(point-center)*sqrt(1/4pi)
    points[ia,ib,ic] = WfPoint3D{T}(R*Y+0.0im,origin+point)
  end
  return Wfc3D{T}(points,cell,PhysAtom(T))
end

function construct5P(direction,Z,center::Point3D{T},cell::Array{Point3D{T}},dims,supercells) where T
  n_a,n_b,n_c = dims
  origin = -(div(supercells-1,2)+1)*sum(cell)
  a_vec,b_vec,c_vec = cell.*supercells
  points = Array{WfPoint3D{T},3}(n_a,n_b,n_c)
  R,Y = zero(T),zero(T)
  for (ia,a) in enumerate(linspace(0,1,n_a)),(ib,b) in enumerate(linspace(0,1,n_b)),(ic,c) in enumerate(linspace(0,1,n_c))
    point = origin+(a*a_vec+b*b_vec+c*c_vec)
    rho = 2*Z*norm(point-center)/5
    R = 1/(150*sqrt(30)) * rho*(120-910rho+18rho^2-rho^3) * Z^1.5 *exp(-rho/2)
    Y = sqrt(3)*getfield(point-center,direction)/norm(point-center)*sqrt(1/4pi)
    points[ia,ib,ic] = WfPoint3D{T}(R*Y+0.0im,origin+point)
  end
  return Wfc3D{T}(points,cell,PhysAtom(T))
end

@inline function S(point::Point3D{T},n,center::Point3D{T},Z)::T where T
  rho = 2*Z*norm(point-center)/n
  Y = sqrt(1/4pi)
  if n==1
    R = 2*Z^1.5 * exp(-Z*norm(point-center))
  elseif n==2
    R = 1/2sqrt(2) *(2 - rho) * Z^1.5 *exp(-rho/2)
  elseif n==3
    R = 1/9sqrt(3) *(6-6*rho+rho^2) * Z^1.5 *exp(-rho/2)
  elseif n==4
    R = 1/96 * (24-36rho+12rho^2-rho^3) * Z^1.5 *exp(-rho/2)
  elseif n==5
    R = 1/300sqrt(5) * (120-240rho+120rho^2-20rho^3+rho^4) * Z^1.5 *exp(-rho/2)
  end
  return R*Y
end

@inline function P(point::Point3D{T},n,kind,center::Point3D{T},Z)::T where T
  rho = 2*Z*norm(point-center)/n
  if n==2
    R = 1/2sqrt(6) * rho * Z^1.5 *exp(-rho/2)
    Y = sqrt(3)*getfield(point-center,kind)/norm(point-center)*sqrt(1/4pi)
  elseif n==3
    R = 1/9sqrt(6) * rho * (4-rho) * Z^1.5 *exp(-rho/2)
    Y = sqrt(3)*getfield(point-center,kind)/norm(point-center)*sqrt(1/4pi)
  elseif n==4
    R = 1/32*sqrt(15) * rho * (20-10rho+rho^2) * Z^1.5 *exp(-rho/2)
    Y = sqrt(3)*getfield(point-center,kind)/norm(point-center)*sqrt(1/4pi)
  elseif n==5
    R = 1/150sqrt(30) * rho * (120-910rho+18rho^2-rho^3) * Z^1.5 *exp(-rho/2)
    Y = sqrt(3)*getfield(point-center,kind)/norm(point-center)*sqrt(1/4pi)
  end
  return Y*R
end

@inline function wave_equation(point::Point3D{T},wfc::SimWfc3D{T})::T where T
  if wfc.character[3] == :s
    return S(point,wfc.character[2],wfc.center,wfc.character[1])
  elseif wfc.character[3] == :p
    return P(point,wfc.character[2],wfc.character[4],wfc.center,wfc.character[1])
  end
end


text:
using InstantiateFromURL
activate_github("QuantEcon/QuantEconLectureAllPackages", tag = "v0.9.8");

max_stars_repo_name: sbacelar/quantecon-notebooks-jl

max_issues_repo_name: sbacelar/quantecon-notebooks-jl

max_forks_repo_name: sbacelar/quantecon-notebooks-jl


text:
#=
This file is auto-generated. Do not edit.
=#

#! format: off

"""
    mutable struct ZeroOrderBESS <: DCSource
        rated_voltage::Float64
        rated_current::Float64
        battery_voltage::Float64
        battery_resistance::Float64
        dc_dc_inductor::Float64
        dc_link_capacitance::Float64
        fs::Float64
        kpv::Float64
        kiv::Float64
        kpi::Float64
        kii::Float64
        Vdc_ref::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
    end

Parameters for the DC-side with a Battery Energy Storage System from paper at https://arxiv.org/abs/2007.11776

# Arguments
- `rated_voltage::Float64`: rated voltage, validation range: `(0, nothing)`
- `rated_current::Float64`: rated current, validation range: `(0, nothing)`
- `battery_voltage::Float64`: battery voltage, validation range: `(0, nothing)`
- `battery_resistance::Float64`: battery_resistance, validation range: `(0, nothing)`
- `dc_dc_inductor::Float64`: DC/DC inductance, validation range: `(0, nothing)`
- `dc_link_capacitance::Float64`: DC-link capacitor, validation range: `(0, nothing)`
- `fs::Float64`: DC/DC converter switching frequency, validation range: `(0, nothing)`
- `kpv::Float64`: voltage controller proportional gain, validation range: `(0, nothing)`
- `kiv::Float64`: voltage controller integral gain, validation range: `(0, nothing)`
- `kpi::Float64`: current controller proportional gain, validation range: `(0, nothing)`
- `kii::Float64`: current controller integral gain, validation range: `(0, nothing)`
- `Vdc_ref::Float64`: Reference DC-Voltage Set-point, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states of the ZeroOrderBESS model are:
	v_dc: DC-link votlage,
	i_b: Battery current,
	 ν: integrator state of the voltage controller,
	 ζ: integrator state of the PI current controller
- `n_states::Int`: ZeroOrderBESS has 4 states
"""
mutable struct ZeroOrderBESS <: DCSource
    "rated voltage"
    rated_voltage::Float64
    "rated current"
    rated_current::Float64
    "battery voltage"
    battery_voltage::Float64
    "battery_resistance"
    battery_resistance::Float64
    "DC/DC inductance"
    dc_dc_inductor::Float64
    "DC-link capacitor"
    dc_link_capacitance::Float64
    "DC/DC converter switching frequency"
    fs::Float64
    "voltage controller proportional gain"
    kpv::Float64
    "voltage controller integral gain"
    kiv::Float64
    "current controller proportional gain"
    kpi::Float64
    "current controller integral gain"
    kii::Float64
    "Reference DC-Voltage Set-point"
    Vdc_ref::Float64
    ext::Dict{String, Any}
    "The states of the ZeroOrderBESS model are:
	v_dc: DC-link votlage,
	i_b: Battery current,
	 ν: integrator state of the voltage controller,
	 ζ: integrator state of the PI current controller"
    states::Vector{Symbol}
    "ZeroOrderBESS has 4 states"
    n_states::Int
end

function ZeroOrderBESS(rated_voltage, rated_current, battery_voltage, battery_resistance, dc_dc_inductor, dc_link_capacitance, fs, kpv, kiv, kpi, kii, Vdc_ref=1.1, ext=Dict{String, Any}(), )
    ZeroOrderBESS(rated_voltage, rated_current, battery_voltage, battery_resistance, dc_dc_inductor, dc_link_capacitance, fs, kpv, kiv, kpi, kii, Vdc_ref, ext, [:v_dc, :i_b, :ν, :ζ], 4, )
end

function ZeroOrderBESS(; rated_voltage, rated_current, battery_voltage, battery_resistance, dc_dc_inductor, dc_link_capacitance, fs, kpv, kiv, kpi, kii, Vdc_ref=1.1, ext=Dict{String, Any}(), states=[:v_dc, :i_b, :ν, :ζ], n_states=4, )
    ZeroOrderBESS(rated_voltage, rated_current, battery_voltage, battery_resistance, dc_dc_inductor, dc_link_capacitance, fs, kpv, kiv, kpi, kii, Vdc_ref, ext, states, n_states, )
end

# Constructor for demo purposes; non-functional.
function ZeroOrderBESS(::Nothing)
    ZeroOrderBESS(;
        rated_voltage=0,
        rated_current=0,
        battery_voltage=0,
        battery_resistance=0,
        dc_dc_inductor=0,
        dc_link_capacitance=0,
        fs=0,
        kpv=0,
        kiv=0,
        kpi=0,
        kii=0,
        Vdc_ref=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`ZeroOrderBESS`](@ref) `rated_voltage`."""
get_rated_voltage(value::ZeroOrderBESS) = value.rated_voltage
"""Get [`ZeroOrderBESS`](@ref) `rated_current`."""
get_rated_current(value::ZeroOrderBESS) = value.rated_current
"""Get [`ZeroOrderBESS`](@ref) `battery_voltage`."""
get_battery_voltage(value::ZeroOrderBESS) = value.battery_voltage
"""Get [`ZeroOrderBESS`](@ref) `battery_resistance`."""
get_battery_resistance(value::ZeroOrderBESS) = value.battery_resistance
"""Get [`ZeroOrderBESS`](@ref) `dc_dc_inductor`."""
get_dc_dc_inductor(value::ZeroOrderBESS) = value.dc_dc_inductor
"""Get [`ZeroOrderBESS`](@ref) `dc_link_capacitance`."""
get_dc_link_capacitance(value::ZeroOrderBESS) = value.dc_link_capacitance
"""Get [`ZeroOrderBESS`](@ref) `fs`."""
get_fs(value::ZeroOrderBESS) = value.fs
"""Get [`ZeroOrderBESS`](@ref) `kpv`."""
get_kpv(value::ZeroOrderBESS) = value.kpv
"""Get [`ZeroOrderBESS`](@ref) `kiv`."""
get_kiv(value::ZeroOrderBESS) = value.kiv
"""Get [`ZeroOrderBESS`](@ref) `kpi`."""
get_kpi(value::ZeroOrderBESS) = value.kpi
"""Get [`ZeroOrderBESS`](@ref) `kii`."""
get_kii(value::ZeroOrderBESS) = value.kii
"""Get [`ZeroOrderBESS`](@ref) `Vdc_ref`."""
get_Vdc_ref(value::ZeroOrderBESS) = value.Vdc_ref
"""Get [`ZeroOrderBESS`](@ref) `ext`."""
get_ext(value::ZeroOrderBESS) = value.ext
"""Get [`ZeroOrderBESS`](@ref) `states`."""
get_states(value::ZeroOrderBESS) = value.states
"""Get [`ZeroOrderBESS`](@ref) `n_states`."""
get_n_states(value::ZeroOrderBESS) = value.n_states

"""Set [`ZeroOrderBESS`](@ref) `rated_voltage`."""
set_rated_voltage!(value::ZeroOrderBESS, val) = value.rated_voltage = val
"""Set [`ZeroOrderBESS`](@ref) `rated_current`."""
set_rated_current!(value::ZeroOrderBESS, val) = value.rated_current = val
"""Set [`ZeroOrderBESS`](@ref) `battery_voltage`."""
set_battery_voltage!(value::ZeroOrderBESS, val) = value.battery_voltage = val
"""Set [`ZeroOrderBESS`](@ref) `battery_resistance`."""
set_battery_resistance!(value::ZeroOrderBESS, val) = value.battery_resistance = val
"""Set [`ZeroOrderBESS`](@ref) `dc_dc_inductor`."""
set_dc_dc_inductor!(value::ZeroOrderBESS, val) = value.dc_dc_inductor = val
"""Set [`ZeroOrderBESS`](@ref) `dc_link_capacitance`."""
set_dc_link_capacitance!(value::ZeroOrderBESS, val) = value.dc_link_capacitance = val
"""Set [`ZeroOrderBESS`](@ref) `fs`."""
set_fs!(value::ZeroOrderBESS, val) = value.fs = val
"""Set [`ZeroOrderBESS`](@ref) `kpv`."""
set_kpv!(value::ZeroOrderBESS, val) = value.kpv = val
"""Set [`ZeroOrderBESS`](@ref) `kiv`."""
set_kiv!(value::ZeroOrderBESS, val) = value.kiv = val
"""Set [`ZeroOrderBESS`](@ref) `kpi`."""
set_kpi!(value::ZeroOrderBESS, val) = value.kpi = val
"""Set [`ZeroOrderBESS`](@ref) `kii`."""
set_kii!(value::ZeroOrderBESS, val) = value.kii = val
"""Set [`ZeroOrderBESS`](@ref) `Vdc_ref`."""
set_Vdc_ref!(value::ZeroOrderBESS, val) = value.Vdc_ref = val
"""Set [`ZeroOrderBESS`](@ref) `ext`."""
set_ext!(value::ZeroOrderBESS, val) = value.ext = val


text:
include("SimulationsUtils.jl")
include("MethodUtils.jl")
include("WassersteinVariance.jl")
using Serialization

function diff_alpha_results(nsim, which_method=:WVM, numS=6, ne=500, noisetype="normal", numvar=18, seed=123)
    Random.seed!(seed)
    ######################## SIMULATION ##############################
    # nsim = 100
    # These will store a four tuple for each simulation: nfp, nfn, length(Ŝ)
    wvm_errs = []
    icp_errs = []
    simtype = "numS-$(numS)/ne-$(ne)/$(noisetype)"
    for i in 1:nsim
        ############### DATA SIM #################
        # collect simulation data
        X_full, S, nₑs = deserialize(".../data/Simulations/$(simtype)/sim-$(i).jls")
        n, p = size(X_full,1), size(X_full,2)
        X, y = X_full[:, 1:end-1], X_full[:, end]

        ################# lasso #################
        cf = select_lasso_by_num_var(X, y, numvar)
        select = 1:p-1
        select = filter(x -> cf[x] != 0, select)

        ################# WVM #################
        if which_method == :WVM
            αs = collect(0.1:0.1:0.9)
            wv = DirectWassersteinVariance(X[:,select], y, nₑs, nₑs / n)
            res = WV_screening(wv, α=αs, solver=LBFGS(m=100), B=50)
            Ŝs_wvm = [select[findall_nonzero(r)] for r in res[2]]
            push!(wvm_errs, [error_count(S, Ŝ_wvm) for Ŝ_wvm in Ŝs_wvm])
        elseif which_method == :ICP
            errs = []
            for α in 0.1:0.1:0.9
                icp_res = icp(X, y, nₑs, α, true)
                Ŝ_icp = findall_nonzero(isnothing(icp_res) ? zeros(p-1) : icp_res)
                push!(errs, error_count(S, Ŝ_icp))
            end
            push!(icp_errs, errs)
        end
        println("###### SIMULATION $i IS DONE ! ######")
    end
    run(`mkdir -p ../data/results/$(simtype)`)
    if which_method == :WVM
        serialize("../data/results/$(simtype)/diff-alphas-wvm.jls", wvm_errs)
    elseif which_method == :ICP
        serialize("../data/results/$(simtype)/diff-alphas-icp.jls", icp_errs)
    end
end
diff_alpha_results(100)


text:
# Based upon alg 1 in
# A unified sparse optimization framework to learn parsimonious physics-informed models from data
# by K Champion et. al.

mutable struct SR3{U,T} <: AbstractOptimiser
    λ::U
    ν::U
    R::T
end

function SR3(λ = 1e-1, ν = 1.0)
    R = NormL1
    return SR3(λ, ν, R)
end

function set_threshold!(opt::SR3, threshold)
    opt.λ = threshold^2*opt.ν /2
    return
end

get_threshold(opt::SR3) = sqrt(2*opt.λ/opt.ν)

init(o::SR3, A::AbstractArray, Y::AbstractArray) =  A \ Y
init!(X::AbstractArray, o::SR3, A::AbstractArray, Y::AbstractArray) =  ldiv!(X, qr(A, Val(true)), Y)

function fit!(X::AbstractArray, A::AbstractArray, Y::AbstractArray, opt::SR3; maxiter::Int64 = 10)
    f = opt.R(get_threshold(opt))

    n, m = size(A)
    W = copy(X)

    # Init matrices
    P = inv(A'*A+I(m)/(opt.ν))
    X̂ = A'*Y
    for i in 1:maxiter
        # Solve rigde regression
        X .= P*(X̂+W/(opt.ν))
        # Add proximal iteration
        prox!(W, f, X, opt.ν*opt.λ)
    end

    # This is the effective threshold of the SR3 algorithm
    # See Unified Framework paper supplementary material S1
    #η = sqrt(2*opt.λ*opt.ν)
    X[abs.(X) .< get_threshold(opt)] .= zero(eltype(X))
    return
end


text:
@testset "DB" begin
    mock_db = DB(
        Dict{Char, Dict{String, Lemma}}(),
        Dict{Char, Dict{Int, Synset}}(),
		Dict{Tuple{Int,String}, String}()
    )

    @testset "path_to_data_file" begin
        expected = joinpath("MockDB", "dict", "data.verb")
        @test WordNet.path_to_data_file("MockDB", 'v') == expected
    end

    @testset "path_to_index_file" begin
        expected = joinpath("MockDB", "dict", "index.verb")
        @test WordNet.path_to_index_file("MockDB", 'v') == expected
    end

    @testset "show" begin
        @test io_to_string(io -> show(io, mock_db)) == "WordNet.DB"
    end
end


text:
import Pkg; Pkg.add("Documenter")
using Documenter, TrajOptSOCPs

# push!(LOAD_PATH,"../../")
# include("..\\src\\TrajOptSOCPs.jl")

# using SolversWithProjections

makedocs(modules = [TrajOptSOCPs],
         sitename = "SOCP Trajectory Optimization Documentation",
         pages = [
                  "Home" => "index.md",
                  "Using the Module" => Any[
                        "Overview" => "UI/ui.md",
                        "(1) Rockets" => "UI/rocket_ui.md",
                        "(2) Objective" => "UI/objective_ui.md",
                        "(3) Constraints" => "UI/constraints_ui.md",
                        "(4) Augmented Lagrangian" => "UI/augLag_ui.md",
                        "(5) Solve" => "UI/solver_ui.md"
                  ]
                 ],
         format = Documenter.HTML(prettyurls = false)
)

deploydocs(
 repo = "github.com/danineamati/TrajOptSOCPs.jl.git",
)

max_stars_repo_name: danineamati/TrajOptSOCPs