text:
using Oceananigans.CubedSpheres
using Oceananigans.BoundaryConditions: fill_halo_regions!
using Oceananigans.CubedSpheres: west_halo, east_halo, south_halo, north_halo, fill_horizontal_velocity_halos!, get_face

# Opposite of the `Base.digits` function
# Source: https://stackoverflow.com/a/55529778
function undigits(d; base=10)
    (s, b) = promote(zero(eltype(d)), base)
    mult = one(s)
    for val in d
        s += val * mult
        mult *= b
    end
    return s
end

cs32_filepath = datadep"cubed_sphere_32_grid/cubed_sphere_32_grid.jld2"

for arch in archs
    @testset "Cubed sphere tracer halo exchange [$(typeof(arch))]" begin
        @info "  Testing cubed sphere tracer halo exchange [$(typeof(arch))]..."

        grid = ConformalCubedSphereGrid(cs32_filepath, arch, Nz=1, z=(-1, 0))
        field = CenterField(arch, grid)

        ## We will fill each grid point with a 5-digit integer "fiijj" where
        ## the f digit is the face number, the ii digits are the i index, and
        ## the jj digits are the j index. We then check that the halo exchange
        ## happened correctly.

        face_digit(n) = digits(abs(Int(n)))[5]
        i_digits(n) = digits(abs(Int(n)))[3:4] |> undigits
        j_digits(n) = digits(abs(Int(n)))[1:2] |> undigits

        for (face_number, field_face) in enumerate(faces(field))
            for i in 1:field_face.grid.Nx, j in 1:field_face.grid.Ny
                CUDA.@allowscalar field_face[i, j, 1] = parse(Int, @sprintf("%d%02d%02d", face_number, i, j))
            end
        end

        fill_halo_regions!(field, arch)

        CUDA.allowscalar(true)

        @testset "Source and destination faces are correct" begin
            for (face_number, field_face) in enumerate(faces(field))
                west_halo_vals =  west_halo(field_face, include_corners=false) |> Array
                east_halo_vals =  east_halo(field_face, include_corners=false) |> Array
                south_halo_vals = south_halo(field_face, include_corners=false) |> Array
                north_halo_vals = north_halo(field_face, include_corners=false) |> Array

                @test all(face_digit.(west_halo_vals)  .== grid.face_connectivity[face_number].west.face)
                @test all(face_digit.(east_halo_vals)  .== grid.face_connectivity[face_number].east.face)
                @test all(face_digit.(south_halo_vals) .== grid.face_connectivity[face_number].south.face)
                @test all(face_digit.(north_halo_vals) .== grid.face_connectivity[face_number].north.face)
            end
        end

        @testset "1W halo <- 5N boundary halo exchange" begin
            # Grid point [i, j] = [0, 1] in 1W halo should be from [i, j] = [32, 32] in 5N boundary.
            west_halo_south_value = get_face(field, 1)[0, 1, 1]
            @test face_digit(west_halo_south_value) == 5
            @test i_digits(west_halo_south_value) == 32
            @test j_digits(west_halo_south_value) == 32

            # Grid point [i, j] = [0, 32] in 1W halo should be from [i, j] = [1, 32] in 5N boundary.
            west_halo_north_value = get_face(field, 1)[0, 32, 1]
            @test face_digit(west_halo_north_value) == 5
            @test i_digits(west_halo_north_value) == 1
            @test j_digits(west_halo_north_value) == 32

            west_halo_values = west_halo(get_face(field, 1), include_corners=false)[:] |> Array
            @test all(face_digit.(west_halo_values) .== 5)
            @test all(i_digits.(west_halo_values) .== reverse(1:32))
            @test all(j_digits.(west_halo_values) .== 32)
        end

        @testset "1E halo <- 2W boundary halo exchange" begin
            # Grid point [i, j] = [33, 1] in 1E halo should be from [i, j] = [1, 1] in 2W boundary.
            east_halo_south_value = get_face(field, 1)[33, 1, 1]
            @test face_digit(east_halo_south_value) == 2
            @test i_digits(east_halo_south_value) == 1
            @test j_digits(east_halo_south_value) == 1

            # Grid point [i, j] = [33, 32] in 1E halo should be from [i, j] = [1, 32] in 2W boundary.
            east_halo_north_value = get_face(field, 1)[33, 32, 1]
            @test face_digit(east_halo_north_value) == 2
            @test i_digits(east_halo_north_value) == 1
            @test j_digits(east_halo_north_value) == 32

            east_halo_values = east_halo(get_face(field, 1), include_corners=false)[:] |> Array
            @test all(face_digit.(east_halo_values) .== 2)
            @test all(i_digits.(east_halo_values) .== 1)
            @test all(j_digits.(east_halo_values) .== 1:32)
        end

        @testset "1S halo <- 6N boundary halo exchange" begin
            # Grid point [i, j] = [1, 0] in 1S halo should be from [i, j] = [1, 32] in 6N boundary.
            south_halo_west_value = get_face(field, 1)[1, 0, 1]
            @test face_digit(south_halo_west_value) == 6
            @test i_digits(south_halo_west_value) == 1
            @test j_digits(south_halo_west_value) == 32

            # Grid point [i, j] = [32, 0] in 1S halo should be from [i, j] = [32, 32] in 6N boundary.
            south_halo_east_value = get_face(field, 1)[32, 0, 1]
            @test face_digit(south_halo_east_value) == 6
            @test i_digits(south_halo_east_value) == 32
            @test j_digits(south_halo_east_value) == 32

            south_halo_values = south_halo(get_face(field, 1), include_corners=false)[:] |> Array
            @test all(face_digit.(south_halo_values) .== 6)
            @test all(i_digits.(south_halo_values) .== 1:32)
            @test all(j_digits.(south_halo_values) .== 32)
        end

        @testset "1N halo <- 3W boundary halo exchange" begin
            # Grid point [i, j] = [1, 33] in 1N halo should be from [i, j] = [1, 32] in 3W boundary.
            north_halo_west_value = get_face(field, 1)[1, 33, 1]
            @test face_digit(north_halo_west_value) == 3
            @test i_digits(north_halo_west_value) == 1
            @test j_digits(north_halo_west_value) == 32

            # Grid point [i, j] = [32, 33] in 1N halo should be from [i, j] = [1, 1] in 3W boundary.
            north_halo_east_value = get_face(field, 1)[32, 33, 1]
            @test face_digit(north_halo_east_value) == 3
            @test i_digits(north_halo_east_value) == 1
            @test j_digits(north_halo_east_value) == 1

            north_halo_values = north_halo(get_face(field, 1), include_corners=false)[:] |> Array
            @test all(face_digit.(north_halo_values) .== 3)
            @test all(i_digits.(north_halo_values) .== 1)
            @test all(j_digits.(north_halo_values) .== reverse(1:32))
        end

        @testset "2W halo <- 1E boundary halo exchange" begin
            # Grid point [i, j] = [0, 1] in 2W halo should be from [i, j] = [32, 1] in 1E boundary.
            west_halo_south_value = get_face(field, 2)[0, 1, 1]
            @test face_digit(west_halo_south_value) == 1
            @test i_digits(west_halo_south_value) == 32
            @test j_digits(west_halo_south_value) == 1

            # Grid point [i, j] = [0, 32] in 2W halo should be from [i, j] = [32, 32] in 1E boundary.
            west_halo_north_value = get_face(field, 2)[0, 32, 1]
            @test face_digit(west_halo_north_value) == 1
            @test i_digits(west_halo_north_value) == 32
            @test j_digits(west_halo_north_value) == 32

            west_halo_values = west_halo(get_face(field, 2), include_corners=false)[:] |> Array
            @test all(face_digit.(west_halo_values) .== 1)
            @test all(i_digits.(west_halo_values) .== 32)
            @test all(j_digits.(west_halo_values) .== 1:32)
        end

        @testset "2E halo <- 4S boundary halo exchange" begin
            # Grid point [i, j] = [33, 1] in 2E halo should be from [i, j] = [32, 1] in 4S boundary.
            east_halo_south_value = get_face(field, 2)[33, 1, 1]
            @test face_digit(east_halo_south_value) == 4
            @test i_digits(east_halo_south_value) == 32
            @test j_digits(east_halo_south_value) == 1

            # Grid point [i, j] = [33, 32] in 2E halo should be from [i, j] = [1, 1] in 4S boundary.
            east_halo_north_value = get_face(field, 2)[33, 32, 1]
            @test face_digit(east_halo_north_value) == 4
            @test i_digits(east_halo_north_value) == 1
            @test j_digits(east_halo_north_value) == 1

            east_halo_values = east_halo(get_face(field, 2), include_corners=false)[:] |> Array
            @test all(face_digit.(east_halo_values) .== 4)
            @test all(i_digits.(east_halo_values) .== reverse(1:32))
            @test all(j_digits.(east_halo_values) .== 1)
        end

        @testset "2S halo <- 6E boundary halo exchange" begin
            # Grid point [i, j] = [1, 0] in 2S halo should be from [i, j] = [32, 32] in 6E boundary.
            south_halo_west_value = get_face(field, 2)[1, 0, 1]
            @test face_digit(south_halo_west_value) == 6
            @test i_digits(south_halo_west_value) == 32
            @test j_digits(south_halo_west_value) == 32

            # Grid point [i, j] = [32, 0] in 2S halo should be from [i, j] = [32, 1] in 6E boundary.
            south_halo_east_value = get_face(field, 2)[32, 0, 1]
            @test face_digit(south_halo_east_value) == 6
            @test i_digits(south_halo_east_value) == 32
            @test j_digits(south_halo_east_value) == 1

            south_halo_values = south_halo(get_face(field, 2), include_corners=false)[:] |> Array
            @test all(face_digit.(south_halo_values) .== 6)
            @test all(i_digits.(south_halo_values) .== 32)
            @test all(j_digits.(south_halo_values) .== reverse(1:32))
        end

        @testset "2N halo <- 3S boundary halo exchange" begin
            # Grid point [i, j] = [1, 33] in 2N halo should be from [i, j] = [1, 1] in 3S boundary.
            north_halo_west_value = get_face(field, 2)[1, 33, 1]
            @test face_digit(north_halo_west_value) == 3
            @test i_digits(north_halo_west_value) == 1
            @test j_digits(north_halo_west_value) == 1

            # Grid point [i, j] = [32, 33] in 2N halo should be from [i, j] = [32, 1] in 3S boundary.
            north_halo_east_value = get_face(field, 2)[32, 33, 1]
            @test face_digit(north_halo_east_value) == 3
            @test i_digits(north_halo_east_value) == 32
            @test j_digits(north_halo_east_value) == 1

            north_halo_values = north_halo(get_face(field, 2), include_corners=false)[:] |> Array
            @test all(face_digit.(north_halo_values) .== 3)
            @test all(i_digits.(north_halo_values) .== 1:32)
            @test all(j_digits.(north_halo_values) .== 1)
        end

        @testset "3W halo <- 1N boundary halo exchange" begin
            # Grid point [i, j] = [0, 1] in 3W halo should be from [i, j] = [32, 32] in 1N boundary.
            west_halo_south_value = get_face(field, 3)[0, 1, 1]
            @test face_digit(west_halo_south_value) == 1
            @test i_digits(west_halo_south_value) == 32
            @test j_digits(west_halo_south_value) == 32

            # Grid point [i, j] = [0, 32] in 3W halo should be from [i, j] = [1, 32] in 1N boundary.
            west_halo_north_value = get_face(field, 3)[0, 32, 1]
            @test face_digit(west_halo_north_value) == 1
            @test i_digits(west_halo_north_value) == 1
            @test j_digits(west_halo_north_value) == 32

            west_halo_values = west_halo(get_face(field, 3), include_corners=false)[:] |> Array
            @test all(face_digit.(west_halo_values) .== 1)
            @test all(i_digits.(west_halo_values) .== reverse(1:32))
            @test all(j_digits.(west_halo_values) .== 32)
        end

        @testset "3E halo <- 4W boundary halo exchange" begin
            # Grid point [i, j] = [33, 1] in 3E halo should be from [i, j] = [1, 1] in 4W boundary.
            east_halo_south_value = get_face(field, 3)[33, 1, 1]
            @test face_digit(east_halo_south_value) == 4
            @test i_digits(east_halo_south_value) == 1
            @test j_digits(east_halo_south_value) == 1

            # Grid point [i, j] = [33, 32] in 3E halo should be from [i, j] = [1, 32] in 4W boundary.
            east_halo_north_value = get_face(field, 3)[33, 32, 1]
            @test face_digit(east_halo_north_value) == 4
            @test i_digits(east_halo_north_value) == 1
            @test j_digits(east_halo_north_value) == 32

            east_halo_values = east_halo(get_face(field, 3), include_corners=false)[:] |> Array
            @test all(face_digit.(east_halo_values) .== 4)
            @test all(i_digits.(east_halo_values) .== 1)
            @test all(j_digits.(east_halo_values) .== 1:32)
        end

        @testset "3S halo <- 2N boundary halo exchange" begin
            # Grid point [i, j] = [1, 0] in 3S halo should be from [i, j] = [1, 32] in 2N boundary.
            south_halo_west_value = get_face(field, 3)[1, 0, 1]
            @test face_digit(south_halo_west_value) == 2
            @test i_digits(south_halo_west_value) == 1
            @test j_digits(south_halo_west_value) == 32

            # Grid point [i, j] = [32, 0] in 3S halo should be from [i, j] = [32, 32] in 2N boundary.
            south_halo_east_value = get_face(field, 3)[32, 0, 1]
            @test face_digit(south_halo_east_value) == 2
            @test i_digits(south_halo_east_value) == 32
            @test j_digits(south_halo_east_value) == 32

            south_halo_values = south_halo(get_face(field, 3), include_corners=false)[:] |> Array
            @test all(face_digit.(south_halo_values) .== 2)
            @test all(i_digits.(south_halo_values) .== 1:32)
            @test all(j_digits.(south_halo_values) .== 32)
        end

        @testset "3N halo <- 5W boundary halo exchange" begin
            # Grid point [i, j] = [1, 33] in 3N halo should be from [i, j] = [1, 32] in 5W boundary.
            north_halo_west_value = get_face(field, 3)[1, 33, 1]
            @test face_digit(north_halo_west_value) == 5
            @test i_digits(north_halo_west_value) == 1
            @test j_digits(north_halo_west_value) == 32

            # Grid point [i, j] = [32, 33] in 3N halo should be from [i, j] = [1, 1] in 5W boundary.
            north_halo_east_value = get_face(field, 3)[32, 33, 1]
            @test face_digit(north_halo_east_value) == 5
            @test i_digits(north_halo_east_value) == 1
            @test j_digits(north_halo_east_value) == 1

            north_halo_values = north_halo(get_face(field, 3), include_corners=false)[:] |> Array
            @test all(face_digit.(north_halo_values) .== 5)
            @test all(i_digits.(north_halo_values) .== 1)
            @test all(j_digits.(north_halo_values) .== reverse(1:32))
        end

        @testset "4W halo <- 3E boundary halo exchange" begin
            # Grid point [i, j] = [0, 1] in 4W halo should be from [i, j] = [32, 1] in 3E boundary.
            west_halo_south_value = get_face(field, 4)[0, 1, 1]
            @test face_digit(west_halo_south_value) == 3
            @test i_digits(west_halo_south_value) == 32
            @test j_digits(west_halo_south_value) == 1

            # Grid point [i, j] = [0, 32] in 4W halo should be from [i, j] = [32, 32] in 1N boundary.
            west_halo_north_value = get_face(field, 4)[0, 32, 1]
            @test face_digit(west_halo_north_value) == 3
            @test i_digits(west_halo_north_value) == 32
            @test j_digits(west_halo_north_value) == 32

            west_halo_values = west_halo(get_face(field, 4), include_corners=false)[:] |> Array
            @test all(face_digit.(west_halo_values) .== 3)
            @test all(i_digits.(west_halo_values) .== 32)
            @test all(j_digits.(west_halo_values) .== 1:32)
        end

        @testset "4E halo <- 6S boundary halo exchange" begin
            # Grid point [i, j] = [33, 1] in 4E halo should be from [i, j] = [32, 1] in 6S boundary.
            east_halo_south_value = get_face(field, 4)[33, 1, 1]
            @test face_digit(east_halo_south_value) == 6
            @test i_digits(east_halo_south_value) == 32
            @test j_digits(east_halo_south_value) == 1

            # Grid point [i, j] = [33, 32] in 4E halo should be from [i, j] = [1, 1] in 6S boundary.
            east_halo_north_value = get_face(field, 4)[33, 32, 1]
            @test face_digit(east_halo_north_value) == 6
            @test i_digits(east_halo_north_value) == 1
            @test j_digits(east_halo_north_value) == 1

            east_halo_values = east_halo(get_face(field, 4), include_corners=false)[:] |> Array
            @test all(face_digit.(east_halo_values) .== 6)
            @test all(i_digits.(east_halo_values) .== reverse(1:32))
            @test all(j_digits.(east_halo_values) .== 1)
        end

        @testset "4S halo <- 2E boundary halo exchange" begin
            # Grid point [i, j] = [1, 0] in 4S halo should be from [i, j] = [32, 32] in 2E boundary.
            south_halo_west_value = get_face(field, 4)[1, 0, 1]
            @test face_digit(south_halo_west_value) == 2
            @test i_digits(south_halo_west_value) == 32
            @test j_digits(south_halo_west_value) == 32

            # Grid point [i, j] = [32, 0] in 4S halo should be from [i, j] = [32, 1] in 2E boundary.
            south_halo_east_value = get_face(field, 4)[32, 0, 1]
            @test face_digit(south_halo_east_value) == 2
            @test i_digits(south_halo_east_value) == 32
            @test j_digits(south_halo_east_value) == 1

            south_halo_values = south_halo(get_face(field, 4), include_corners=false)[:] |> Array
            @test all(face_digit.(south_halo_values) .== 2)
            @test all(i_digits.(south_halo_values) .== 32)
            @test all(j_digits.(south_halo_values) .== reverse(1:32))
        end

        @testset "4N halo <- 5S boundary halo exchange" begin
            # Grid point [i, j] = [1, 33] in 4N halo should be from [i, j] = [1, 1] in 5S boundary.
            north_halo_west_value = get_face(field, 4)[1, 33, 1]
            @test face_digit(north_halo_west_value) == 5
            @test i_digits(north_halo_west_value) == 1
            @test j_digits(north_halo_west_value) == 1

            # Grid point [i, j] = [32, 33] in 4N halo should be from [i, j] = [32, 1] in 5S boundary.
            north_halo_east_value = get_face(field, 4)[32, 33, 1]
            @test face_digit(north_halo_east_value) == 5
            @test i_digits(north_halo_east_value) == 32
            @test j_digits(north_halo_east_value) == 1

            north_halo_values = north_halo(get_face(field, 4), include_corners=false)[:] |> Array
            @test all(face_digit.(north_halo_values) .== 5)
            @test all(i_digits.(north_halo_values) .== 1:32)
            @test all(j_digits.(north_halo_values) .== 1)
        end

        ## TODO: Test faces 5 and 6 (or generalize the test to all faces)!

        CUDA.allowscalar(false)

    end
end

for arch in archs
    @testset "Cubed sphere velocity halo exchange [$(typeof(arch))]" begin
        @info "  Testing cubed sphere velocity halo exchange [$(typeof(arch))]..."

        grid = ConformalCubedSphereGrid(cs32_filepath, arch, Nz=1, z=(-1, 0))

        u_field = XFaceField(arch, grid)
        v_field = YFaceField(arch, grid)

        ## We will fill each grid point with a 6-digit integer "ufiijj" where
        ## the u digit is 1 for u and 2 for v, the f digit is the face number,
        ## the ii digits are the i index, and the jj digits are the j index.
        ## We then check that the halo exchange happened correctly (including
        ## any sign changes).

        U_DIGIT = 1
        V_DIGIT = 2

        uv_digit(n) = digits(abs(Int(n)))[6]
        face_digit(n) = digits(abs(Int(n)))[5]
        i_digits(n) = digits(abs(Int(n)))[3:4] |> undigits
        j_digits(n) = digits(abs(Int(n)))[1:2] |> undigits

        for (face_number, u_field_face) in enumerate(faces(u_field))
            for i in 1:u_field_face.grid.Nx+1, j in 1:u_field_face.grid.Ny
                CUDA.@allowscalar u_field_face[i, j, 1] = parse(Int, @sprintf("%d%d%02d%02d", U_DIGIT, face_number, i, j))
            end
        end

        for (face_number, v_field_face) in enumerate(faces(v_field))
            for i in 1:v_field_face.grid.Nx, j in 1:v_field_face.grid.Ny+1
                CUDA.@allowscalar v_field_face[i, j, 1] = parse(Int, @sprintf("%d%d%02d%02d", V_DIGIT, face_number, i, j))
            end
        end

        fill_horizontal_velocity_halos!(u_field, v_field, arch)

        CUDA.allowscalar(true)

        # @testset "Source and destination faces are correct" begin
        #     for field in (u_field, v_field), (face_number, field_face) in enumerate(faces(field))
        #         west_halo_vals =  west_halo(field_face, include_corners=false)
        #         east_halo_vals =  east_halo(field_face, include_corners=false)
        #         south_halo_vals = south_halo(field_face, include_corners=false)
        #         north_halo_vals = north_halo(field_face, include_corners=false)

        #         @test all(face_digit.(west_halo_vals)  .== grid.face_connectivity[face_number].west.face)
        #         @test all(face_digit.(east_halo_vals)  .== grid.face_connectivity[face_number].east.face)
        #         @test all(face_digit.(south_halo_vals) .== grid.face_connectivity[face_number].south.face)
        #         @test all(face_digit.(north_halo_vals) .== grid.face_connectivity[face_number].north.face)
        #     end
        # end

        @testset "(1W halo <- 5N boundary) horizontal velocity halo exchange" begin
            # Grid point u[i, j] = u[0, 1] in 1W halo should be from +v[i, j] = +v[32, 32] in 5N boundary.
            u_west_halo_south_value = get_face(u_field, 1)[0, 1, 1]
            @test uv_digit(u_west_halo_south_value) == V_DIGIT
            @test sign(u_west_halo_south_value) == +1
            @test face_digit(u_west_halo_south_value) == 5
            @test i_digits(u_west_halo_south_value) == 32
            @test j_digits(u_west_halo_south_value) == 32

            # Grid point u[i, j] = u[0, 32] in 1W halo should be from +v[i, j] = +v[1, 32] in 5N boundary.
            u_west_halo_north_value = get_face(u_field, 1)[0, 32, 1]
            @test uv_digit(u_west_halo_north_value) == V_DIGIT
            @test sign(u_west_halo_north_value) == +1
            @test face_digit(u_west_halo_north_value) == 5
            @test i_digits(u_west_halo_north_value) == 1
            @test j_digits(u_west_halo_north_value) == 32

            u_west_halo_values = west_halo(get_face(u_field, 1), include_corners=false)[:] |> Array
            @test all(uv_digit.(u_west_halo_values) .== V_DIGIT)
            @test all(sign.(u_west_halo_values) .== +1)
            @test all(face_digit.(u_west_halo_values) .== 5)
            @test all(i_digits.(u_west_halo_values) .== reverse(1:32))
            @test all(j_digits.(u_west_halo_values) .== 32)

            # Note: The 1W v halo <- 5N -u boundary exchange involves a sign change and an index shift.

            # Grid point v[i, j] = v[0, 1] in 1W halo should be from -u[i, j] = -u[1, 32] in 6NE corner!
            v_west_halo_south_value = get_face(v_field, 1)[0, 1, 1]
            @test uv_digit(v_west_halo_south_value) == U_DIGIT
            @test sign(v_west_halo_south_value) == -1
            @test face_digit(v_west_halo_south_value) == 6
            @test i_digits(v_west_halo_south_value) == 1
            @test j_digits(v_west_halo_south_value) == 32

            # Grid point v[i, j] = u[0, 33] in 1W halo should be from -u[i, j] = -u[2, 32] in 5N boundary.
            v_west_halo_north_value = get_face(v_field, 1)[0, 32, 1]
            @test uv_digit(v_west_halo_north_value) == U_DIGIT
            @test sign(v_west_halo_north_value) == -1
            @test face_digit(v_west_halo_north_value) == 5
            @test i_digits(v_west_halo_north_value) == 2
            @test j_digits(v_west_halo_north_value) == 32

            v_west_halo_values = west_halo(get_face(v_field, 1), include_corners=false)[:] |> Array
            @test all(uv_digit.(v_west_halo_values) .== U_DIGIT)
            @test all(sign.(v_west_halo_values) .== -1)
            # @test all(face_digit.(v_west_halo_values) .== 5)
            # @test all(i_digits.(v_west_halo_values) .== reverse(2:33))
            @test all(j_digits.(v_west_halo_values) .== 32)
        end

        @testset "(1E halo <- 2W boundary) horizontal velocity halo exchange" begin
            # Grid point u[i, j] = u[33, 1] in 1E halo should be from +u[i, j] = +u[1, 1] in 2W boundary.
            u_east_halo_south_value = get_face(u_field, 1)[33, 1, 1]
            @test uv_digit(u_east_halo_south_value) == U_DIGIT
            @test sign(u_east_halo_south_value) == +1
            @test face_digit(u_east_halo_south_value) == 2
            @test i_digits(u_east_halo_south_value) == 1
            @test j_digits(u_east_halo_south_value) == 1

            # Grid point u[i, j] = u[33, 32] in 1E halo should be from +u[i, j] = +u[1, 32] in 2W boundary.
            u_east_halo_north_value = get_face(u_field, 1)[33, 32, 1]
            @test uv_digit(u_east_halo_north_value) == U_DIGIT
            @test sign(u_east_halo_north_value) == +1
            @test face_digit(u_east_halo_north_value) == 2
            @test i_digits(u_east_halo_north_value) == 1
            @test j_digits(u_east_halo_north_value) == 32

            u_east_halo_values = east_halo(get_face(u_field, 1), include_corners=false)[:] |> Array
            @test all(uv_digit.(u_east_halo_values) .== U_DIGIT)
            @test all(sign.(u_east_halo_values) .== +1)
            @test all(face_digit.(u_east_halo_values) .== 2)
            @test all(i_digits.(u_east_halo_values) .== 1)
            @test all(j_digits.(u_east_halo_values) .== 1:32)

            # Grid point v[i, j] = v[33, 1] in 1E halo should be from +v[i, j] = +v[1, 1] in 2W boundary.
            v_east_halo_south_value = get_face(v_field, 1)[33, 1, 1]
            @test uv_digit(v_east_halo_south_value) == V_DIGIT
            @test sign(v_east_halo_south_value) == +1
            @test face_digit(v_east_halo_south_value) == 2
            @test i_digits(v_east_halo_south_value) == 1
            @test j_digits(v_east_halo_south_value) == 1

            # # Grid point v[i, j] = v[33, 32] in 1E halo should be from +v[i, j] = +v[1, 32] in 2W boundary.
            v_east_halo_north_value = get_face(v_field, 1)[33, 32, 1]
            @test uv_digit(v_east_halo_north_value) == V_DIGIT
            @test sign(v_east_halo_north_value) == +1
            @test face_digit(v_east_halo_north_value) == 2
            @test i_digits(v_east_halo_north_value) == 1
            @test j_digits(v_east_halo_north_value) == 32

            v_east_halo_values = east_halo(get_face(v_field, 1), include_corners=false)[:] |> Array
            @test all(uv_digit.(v_east_halo_values) .== V_DIGIT)
            @test all(sign.(v_east_halo_values) .== +1)
            @test all(face_digit.(v_east_halo_values) .== 2)
            @test all(i_digits.(v_east_halo_values) .== 1)
            @test all(j_digits.(v_east_halo_values) .== 1:32)
        end

        @testset "(1S halo <- 6N boundary) horizontal velocity halo exchange" begin
            # Grid point u[i, j] = u[1, 0] in 1S halo should be from +u[i, j] = +u[1, 32] in 6N boundary.
            u_south_halo_west_value = get_face(u_field, 1)[1, 0, 1]
            @test uv_digit(u_south_halo_west_value) == U_DIGIT
            @test sign(u_south_halo_west_value) == +1
            @test face_digit(u_south_halo_west_value) == 6
            @test i_digits(u_south_halo_west_value) == 1
            @test j_digits(u_south_halo_west_value) == 32

            # Grid point u[i, j] = u[32, 0] in 1S halo should be from +u[i, j] = +u[32, 32] in 6N boundary.
            u_south_halo_east_value = get_face(u_field, 1)[32, 0, 1]
            @test uv_digit(u_south_halo_east_value) == U_DIGIT
            @test sign(u_south_halo_east_value) == +1
            @test face_digit(u_south_halo_east_value) == 6
            @test i_digits(u_south_halo_east_value) == 32
            @test j_digits(u_south_halo_east_value) == 32

            u_south_halo_values = south_halo(get_face(u_field, 1), include_corners=false)[:] |> Array
            @test all(uv_digit.(u_south_halo_values) .== U_DIGIT)
            @test all(sign.(u_south_halo_values) .== +1)
            @test all(face_digit.(u_south_halo_values) .== 6)
            @test all(i_digits.(u_south_halo_values) .== 1:32)
            @test all(j_digits.(u_south_halo_values) .== 32)

            # Grid point v[i, j] = v[1, 0] in 1S halo should be from +v[i, j] = +v[1, 32] in 6N boundary.
            v_south_halo_west_value = get_face(v_field, 1)[1, 0, 1]
            @test uv_digit(v_south_halo_west_value) == V_DIGIT
            @test sign(v_south_halo_west_value) == +1
            @test face_digit(v_south_halo_west_value) == 6
            @test i_digits(v_south_halo_west_value) == 1
            @test j_digits(v_south_halo_west_value) == 32

            # # Grid point v[i, j] = v[32, 0] in 1S halo should be from +v[i, j] = +v[32, 32] in 6N boundary.
            v_south_halo_east_value = get_face(v_field, 1)[32, 0, 1]
            @test uv_digit(v_south_halo_east_value) == V_DIGIT
            @test sign(v_south_halo_east_value) == +1
            @test face_digit(v_south_halo_east_value) == 6
            @test i_digits(v_south_halo_east_value) == 32
            @test j_digits(v_south_halo_east_value) == 32

            v_south_halo_values = south_halo(get_face(v_field, 1), include_corners=false)[:] |> Array
            @test all(uv_digit.(v_south_halo_values) .== V_DIGIT)
            @test all(sign.(v_south_halo_values) .== +1)
            @test all(face_digit.(v_south_halo_values) .== 6)
            @test all(i_digits.(v_south_halo_values) .== 1:32)
            @test all(j_digits.(v_south_halo_values) .== 32)
        end

        @testset "(1N halo <- 3W boundary) horizontal velocity halo exchange" begin

            # Note: The 1N u halo <- 3W -v boundary exchange involves a sign change and an index shift.

            # Grid point u[i, j] = u[1, 33] in 1N halo should be from -v[i, j] = -u[1, 32] in 5NE corner!
            u_nouth_halo_west_value = get_face(u_field, 1)[1, 33, 1]
            @test uv_digit(u_nouth_halo_west_value) == U_DIGIT
            @test sign(u_nouth_halo_west_value) == -1
            @test face_digit(u_nouth_halo_west_value) == 5
            @test i_digits(u_nouth_halo_west_value) == 1
            @test j_digits(u_nouth_halo_west_value) == 32

            # Grid point u[i, j] = u[32, 33] in 1N halo should be from -v[i, j] = -v[1, 2] in 3W boundary.
            u_north_halo_east_value = get_face(u_field, 1)[32, 33, 1]
            @test uv_digit(u_north_halo_east_value) == V_DIGIT
            @test sign(u_north_halo_east_value) == -1
            @test face_digit(u_north_halo_east_value) == 3
            @test i_digits(u_north_halo_east_value) == 1
            @test j_digits(u_north_halo_east_value) == 2

            u_north_halo_values = north_halo(get_face(u_field, 1), include_corners=false)[:] |> Array
            # @test all(uv_digit.(u_north_halo_values) .== V_DIGIT)
            @test all(sign.(u_north_halo_values) .== -1)
            # @test all(face_digit.(u_north_halo_values) .== 3)
            @test all(i_digits.(u_north_halo_values) .== 1)
            # @test all(j_digits.(u_north_halo_values) .== reverse(2:33))

            # Grid point v[i, j] = v[1, 33] in 1N halo should be from +u[i, j] = +u[1, 32] in 3W boundary.
            v_north_halo_west_value = get_face(v_field, 1)[1, 33, 1]
            @test uv_digit(v_north_halo_west_value) == U_DIGIT
            @test sign(v_north_halo_west_value) == +1
            @test face_digit(v_north_halo_west_value) == 3
            @test i_digits(v_north_halo_west_value) == 1
            @test j_digits(v_north_halo_west_value) == 32

            # # Grid point v[i, j] = v[32, 33] in 1N halo should be from +u[i, j] = +u[1, 1] in 3W boundary.
            v_north_halo_east_value = get_face(v_field, 1)[32, 33, 1]
            @test uv_digit(v_north_halo_east_value) == U_DIGIT
            @test sign(v_north_halo_east_value) == +1
            @test face_digit(v_north_halo_east_value) == 3
            @test i_digits(v_north_halo_east_value) == 1
            @test j_digits(v_north_halo_east_value) == 1

            v_north_halo_values = north_halo(get_face(v_field, 1), include_corners=false)[:] |> Array
            @test all(uv_digit.(v_north_halo_values) .== U_DIGIT)
            @test all(sign.(v_north_halo_values) .== +1)
            @test all(face_digit.(v_north_halo_values) .== 3)
            @test all(i_digits.(v_north_halo_values) .== 1)
            @test all(j_digits.(v_north_halo_values) .== reverse(1:32))
        end

        ## TODO: Test the other faces, especially face 2? I hope we can generalize this...

        @testset "Velocities at the corners (MITgcm regression)" begin
            # Face 1
            @test get_face(u_field, 1)[1,  33, 1] == -150132
            @test get_face(u_field, 1)[33, 0,  1] ==  220101

            @test get_face(v_field, 1)[0,   1, 1] == -160132
            @test get_face(v_field, 1)[0,  33, 1] == -150132
            @test get_face(v_field, 1)[33, 33, 1] ==  230101

            # Face 2
            @test get_face(u_field, 2)[1,  0,  1] == -213201
            @test get_face(u_field, 2)[33, 33, 1] ==  140101

            @test get_face(v_field, 2)[33,  1, 1] == -263201
            @test get_face(v_field, 2)[0,  33, 1] ==  130101
            @test get_face(v_field, 2)[33, 33, 1] == -140101

            # Face 3
            @test get_face(u_field, 3)[1,  33, 1] == -110132

            @test get_face(v_field, 3)[0,  33, 1] == -110132
            @test get_face(v_field, 3)[33, 33, 1] ==  250101

            # Face 4
            @test get_face(u_field, 4)[33, 33, 1] ==  160101

            @test get_face(v_field, 4)[33,  1, 1] == -223201
            @test get_face(v_field, 4)[0,  33, 1] ==  150101
            @test get_face(v_field, 4)[33, 33, 1] == -160101

            # Face 5
            @test get_face(v_field, 5)[0,  33, 1] == -130132
            @test get_face(v_field, 5)[33, 33, 1] ==  210101

            # Face 6
            @test get_face(v_field, 6)[0,  33, 1] ==  110101
            @test get_face(v_field, 6)[33, 33, 1] == -120101
        end

        CUDA.allowscalar(false)

    end
end


text:
using BulletsMaggots, Statistics, ProgressMeter

function run_exp(c::Int, reps::Int = 1000)
    println("Solving for code $(string4(c)) for $reps iterations")
    t = [bm_solver(c, false) for _ = 1:reps]
    println("Average steps to solve:   ", mean(t))
    println("Median steps to solve:    ", median(t))
    println("Standard deviation:       ", std(t))
    println("Range in number of steps: ", minimum(t), "-", maximum(t))
end

"""
`run_all(reps::Int=1000)` runs the solver repeatedly on the five possible
digit patterns for Bullets and Maggots.
"""
function run_all(reps::Int = 1000)
    codes = [0000, 0001, 0011, 0012, 0123]
    for c in codes
        run_exp(c, reps)
        println()
    end
end

function experiment(reps::Int = 1000, first_guess::Function = guess_any)
    println("Playing $reps games")
    t = zeros(Int, reps)
    PM = Progress(reps)
    for i = 1:reps
        t[i] = bm_solver(random_code(), false, first_guess)
        next!(PM)
    end
    println("Average steps to solve:   ", mean(t))
    println("Median steps to solve:    ", median(t))
    println("Standard deviation:       ", std(t))
    println("Range in number of steps: ", minimum(t), "-", maximum(t))
end

function first_move_experiment(reps::Int = 1000)

    sep = "-"^60

    println("Any first guess")
    experiment(reps, guess_any)

    println(sep)

    println("All distinct")
    experiment(reps, guess_distinct)

    println(sep)

    println("One pair")
    experiment(reps, guess_one_pair)

    println(sep)

    println("Two pairs")
    experiment(reps, guess_two_pairs)

    println(sep)
end

ext: jl


text:
# -----------------------------------------------------------------------
# Binary coded genetic algorithm
# -----------------------------------------------------------------------
# Required packages
using BitBasis
using Statistics
include("stepper_move.jl")
include("optim_move.jl")
include("random_init.jl")
include("visual.jl")
include("simdata.jl")
# -----------------------------------------------------------------------
# Important onbservations
# These codes were not actually used in the final experiments
# Instead, the decimal coded algorithm was preferred
# It might not work well.
# The simulator implemented works for one cylinder onlys
# -----------------------------------------------------------------------
# Useful functions
# -----------------------------------------------------------------------

"""
    cylpos2bits(cyl::CylType; size=[5,4,3,3,3,3], coord_zero=[0,-5,0,-4,0,-4])

## Description
Convert a position of any CylType to bit.
## Arguments
- cyl: position of type OneCylType, TwoCylType or ThreeCylType
- size: size of the search 2^size (mm)
- coord_zero: initial coordinates (offset)
## Example
```jldoctest
julia> cylpos2bits(ThreeCylType(10,2,5,0,6,-1))
6-element Array{BitArray{1},1}:
 [0, 1, 0, 1, 0]
 [1, 1, 1, 0]
 [1, 0, 1]
 [0, 0, 1]
 [0, 1, 1]
 [1, 1, 0]
```
"""
function cylpos2bits(cyl::CylType; size=[5,4,3,3,3,3], coord_zero=[0,-5,0,-4,0,-4])
    a = convert(Array, cyl)
    n = length(a)
    id = 1:n

    # Applying offset on the coordinates
    @. a[id] = a[id] - coord_zero[id]

    # convert to bit
    a_bit = bitarray.(a, size[id])

    return a_bit
end

"""
    bits2cylpos(b::AbstractArray; coord_zero=[0,-5,0,-4,0,-4])

## Description
Convert a position of any CylType to bit.
## Arguments
- b: position on bits
- coord_zero: initial coordinates (offset)
## Example
```jldoctest
julia> b = cylpos2bits(ThreeCylType(10,2,5,0,6,-1));

julia> bits2cylpos(b)
ThreeCylType(10.0, 2.0, 5.0, 0.0, 6.0, -1.0)
```
"""
function bits2cylpos(b::AbstractArray; coord_zero=[0,-5,0,-4,0,-4])
    n = length(b)
    id = 1:n

    p = packbits.(b)

    # Applying offset on the coordinates
    @. p[id] = p[id] + coord_zero[id]

    if n == 2
        pos = OneCylType(p[1],p[2])
    elseif n == 4
        pos = TwoCylType(p[1],p[2],p[3],p[4])
    elseif n == 6
        pos = ThreeCylType(p[1],p[2],p[3],p[4],p[5],p[6])
    else
        error("Number of cylinders not defined")
    end

    return pos
end

"""
    sigmoid(gen; offset, α)

## Description
Sigmoidal function to control the mutation.
## Arguments
- gen: number of generation
- offset: generation offset parameter (20)
- α: power factor (0.4)
"""
sigmoid(gen; offset, α) = 1 - 1/(1+exp(-α*(gen-offset)))

"""
    dict_init(type=ThreeCylType)

## Description
Initiate a dictionary for the read data.
## Arguments
- type: key type, use OneCylType, TwoCylType or ThreeCylType
"""
dict_init(type=ThreeCylType) = Dict{type, Array{Float64,1}}()

"""
    repeatedpos(dict, pos::CylType)

## Description
Find if the position was not measured. If measured, add one count on the dictionary.
## Arguments
- dict: dictionary with the positions and the data measured
- pos: position of CylType
"""
function repeatedpos(dict::Dict, pos::CylType)
    if haskey(dict, pos)
        dict[pos][2] += 1
        return false
    else
        return true
    end
end

"""
    newposition(dict::Dict, posArray::ArrayCylType)

## Description
Discover which positions of the generation must be measured.
## Arguments
- dict: dictionary with the positions and the data measured
- posArray: array with the positions of CylType
"""
function newposition(dict::Dict, posArray::ArrayCylType)
    id = [repeatedpos(dict,i) for i in posArray]
    new_posArray = posArray[id]
    return new_posArray
end

# -----------------------------------------------------------------------
# Convergence
# -----------------------------------------------------------------------

"""
    criteriarule(new_coord::AbstractVector, old_coord::AbstractVector; μ, σ)

## Description
General criteria rule.
## Arguments
- old_coord: array with floats of the old coordinate
- new_coord: array with floats of the new coordinate
- μ: mean parameter
- σ: standard deviation parameter
"""
function criteriarule(new_coord::AbstractVector, old_coord::AbstractVector; μ, σ)

    μ = float(μ)
    σ = float(σ)

    Δcoord = abs(mean(new_coord)-mean(old_coord))

    if Δcoord < μ && std(new_coord) < σ
        return true
    else
        return false
    end
end

"""
    criteria(new_gen::Array{T,1}, old_gen::Array{T,1}; μ, σ, Δμ, Δσ) where T<:CylType

## Description
Criteria for one, two or three control cylinders.
## Arguments
- old_gen: last generation measured
- new_gen: new generation measured
- μ, Δμ: mean parameters
- σ, Δσ: standard deviation parameters
"""
function criteria(new_gen::Array{OneCylType,1}, old_gen::Array{OneCylType,1}; μ, σ, Δμ=0.0, Δσ=0.0)

    n = length(new_gen)

    x_old = [old_gen[i].x for i in 1:n]
    x_new = [new_gen[i].x for i in 1:n]

    if criteriarule(x_new, x_old , μ=μ, σ=σ)
        y_old = [old_gen[i].y for i in 1:n]
        y_new = [new_gen[i].y for i in 1:n]
        return criteriarule(y_new, y_old, μ=μ, σ=σ)
    else
        return false
    end
end

function criteria(new_gen::Array{TwoCylType,1}, old_gen::Array{TwoCylType,1}; μ, σ, Δμ, Δσ)

    n = length(new_gen)

    x_old = [old_gen[i].x for i in 1:n]
    x_new = [new_gen[i].x for i in 1:n]

    if criteriarule(x_new, x_old , μ=μ, σ=σ)
        y_old = [old_gen[i].y for i in 1:n]
        y_new = [new_gen[i].y for i in 1:n]

        if criteriarule(y_new, y_old, μ=μ, σ=σ)
            dx_old = [old_gen[i].dx for i in 1:n]
            dx_new = [new_gen[i].dx for i in 1:n]

            if criteriarule(dx_new, dx_old, μ=Δμ, σ=Δσ)
                dy_old = [old_gen[i].dy for i in 1:n]
                dy_new = [new_gen[i].dy for i in 1:n]
                return criteriarule(dy_new, dy_old, μ=Δμ, σ=Δσ)
            else
                return false
            end
        else
            return false
        end
    else
        return false
    end
end

function criteria(new_gen::Array{ThreeCylType,1}, old_gen::Array{ThreeCylType,1}; μ, σ, Δμ, Δσ)

    n = length(new_gen)

    x_old = [old_gen[i].x for i in 1:n]
    x_new = [new_gen[i].x for i in 1:n]

    if criteriarule(x_new, x_old , μ=μ, σ=σ)
        y_old = [old_gen[i].y for i in 1:n]
        y_new = [new_gen[i].y for i in 1:n]

        if criteriarule(y_new, y_old, μ=μ, σ=σ)
            dx1_old = [old_gen[i].dx1 for i in 1:n]
            dx1_new = [new_gen[i].dx1 for i in 1:n]

            if criteriarule(dx1_new, dx1_old, μ=Δμ, σ=Δσ)
                dy1_old = [old_gen[i].dy1 for i in 1:n]
                dy1_new = [new_gen[i].dy1 for i in 1:n]

                if criteriarule(dy1_new, dy1_old, μ=Δμ, σ=Δσ)
                    dx2_old = [old_gen[i].dx2 for i in 1:n]
                    dx2_new = [new_gen[i].dx2 for i in 1:n]

                    if criteriarule(dx2_new, dx2_old, μ=Δμ, σ=Δσ)
                        dy2_old = [old_gen[i].dy2 for i in 1:n]
                        dy2_new = [new_gen[i].dy2 for i in 1:n]
                        return criteriarule(dy2_new, dy2_old, μ=Δμ, σ=Δσ)
                    else
                        return false
                    end
                else
                    return false
                end
            else
                return false
            end
        else
            return false
        end
    else
        return false
    end
end

# -----------------------------------------------------------------------
# Genetic Algorithm functions
# -----------------------------------------------------------------------

"""
    remove_individual(rm, pop)

## Description
Remove individual from the population.
## Arguments
- rm: individual removed
- pop: population
"""
function remove_individual(rm, pop)
    id = [i != rm for i in pop]
    return pop[id]
end

"""
    roulette_selection(val, pop)

## Description
Selection of the roulette type.
## Arguments
- val: evaluation of each individual, according to a feature
- pop: population
"""
function roulette_selection(val, pop)
    fitness = round.(cumsum(val/sum(val)), digits=5)
    n = length(val) - 1 # one individual already chosen (eletive selection)
    id = []

    for i in 1:n
        idtemp = sum(rand() .> fitness) + 1
        id = push!(id, idtemp)
    end

    return pop[id]
end

"""
    crossover(parent1::Array{BitArray{1},1}, parent2::Array{BitArray{1},1}; size)

## Description
Crossover.
## Arguments
- parent1: first individual to crossover (bit)
- parent2: second individual to crossover (bit)
- size: bit sizes
"""
function crossover(parent1::Array{BitArray{1},1}, parent2::Array{BitArray{1},1}; size)

    ncoord = length(parent1)
    cut = vcat(rand([1,2,3], 2), rand([1,2],ncoord-2)) # can be changed
    child1 = Array{BitArray{1},1}()
    child2 = Array{BitArray{1},1}()

    for i in 1:ncoord
        k = cut[i] + 1
        push!(child1, vcat(parent1[i][1:cut[i]], parent2[i][k:size[i]]))
        push!(child2, vcat(parent2[i][1:cut[i]], parent1[i][k:size[i]]))
    end

    return [child1, child2]

end

"""
    mutation(ind::Array{BitArray{1},1}; Tm)

## Description
Mutation.
## Arguments
- ind: single individual of the population
- Tm: mutation rate
"""
function mutation(ind::Array{BitArray{1},1}; Tm)

    for i in 1:length(ind)
        idm = rand(length(ind[i])) .< Tm
        ind[i][idm] = .!ind[i][idm]
    end

    return ind
end

"""
    ga_onestep(pop, dict; Tr=0.25, Tm=0.005, size=[5,4,3,3,3,3], coord_zero=[0,-5,0,-4,0,-4])

## Description
Implementation of the genetic algorithm (one iteration).
## Arguments
- pop: population
- Tr: crossover rate
- Tm: mutation rate
- size: size of the search 2^size (mm)
- coord_zero: initial coordinates (offset)
"""
function ga_onestep(pop, dict; Tr=0.25, Tm=0.005, size=[5,4,3,3,3,3], coord_zero=[0,-5,0,-4,0,-4])

    # Fitness of each individual
    val = [dict[i][1] for i in pop]
    #println("fitness ok")

    # Eletive selection
    idmax = argmax(val)
    elite = pop[idmax]
    #println("eletive ok")

    # Roulette selection
    pop_sel = roulette_selection(val, pop)
    #println("roulette ok")

    # Convert to bit
    pop_bit = cylpos2bits.(pop_sel, size=size, coord_zero=coord_zero)
    np_bit = length(pop_bit)
    #println("convert to bit ok")

    # Crossover
    idr = rand(np_bit) .< Tr
    idr_inv = .!(idr)

    if sum(idr) >= 2
        parents = pop_bit[idr]
        pop_bit_r = pop_bit[idr_inv]

        for i in 1:2:sum(idr)
            k = i+1
            if i == sum(idr)
                push!(pop_bit_r, parents[i])
                break
            end
            child1, child2 = crossover(parents[i], parents[k], size=size)
            push!(pop_bit_r, child1, child2)
        end

    else
        pop_bit_r = pop_bit
    end
    #println("crossover ok")

    # Mutation
    pop_bit_m = mutation.(pop_bit_r, Tm=Tm)
    #println("mutation ok")

    # New generation
    new_gen = bits2cylpos.(pop_bit_m, coord_zero=coord_zero)
    push!(new_gen, elite) # adding eletive selection

    #return [elite, pop_sel, pop_bit, pop_bit_r, pop_bit_m, new_gen]
    return new_gen

end

"""
    ga_simulator(;n=30, nc=3, Tr=0.25, Tmmax=0.01, offset=25, α=0.4, μ=1.0, σ=0.5, Δμ=1.0, Δσ=1.0, ngen_max=100)

## Description
GA simulator.
## Arguments
- n: number of individuals in population
- nc: number of cylinders
- Tr: crossover rate
- Tmmax: maximum mutation rate
- offset and α: sigmoid parameters (for mutation decay)
- μ, σ, Δμ, Δσ: convergence criteria parameters
- ngen_max: max number of generations allowed
"""
function ga_simulator(;n=30, nc=3, Tr=0.25, Tmmax=0.01, offset=25, α=0.4, μ=1.0, σ=0.5, Δμ=1.0, Δσ=1.0, ngen_max=100)

    size = [5,4,3,3,3,3]
    coord_zero=[0,-5,0,-4,0,-4]
    pop = population_init(n, nc=nc, size=size[1:2*nc], coord_zero=coord_zero[1:2*nc]) # initializing population
    dict = dict_init(typeof(pop[1])) # initializing dictionary

    data = simDataOne("data_test01") # get simulation data

    ngen = 1 # initializing generation number
    criteriaConv = false

    new_pop = pop
    NI = 0 # total number of individuals (visited at least once)
    NS = 0 # total number of steps
    TT = 0 # total time
    pos = randindividual(nc, size=size) # initial position of the robot

    while ngen < ngen_max && !criteriaConv

        println("Geração: ", ngen)

        toMeasure = newposition(dict, new_pop) # find which individuals must be evaluated
        NI += length(toMeasure)

        if length(toMeasure) > 0
            if length(toMeasure) == 1
                toMeasure_opt = toMeasure
            else
                toMeasure_opt = toMeasure[greedypath(pos, toMeasure)]
            end

            #toMeasure_opt = toMeasure # not applying path optimization

            for i in 1:length(toMeasure_opt)
                NS += sum(abs.(coord2steps(pos,toMeasure_opt[i])[1].step)) # find total steps for each individual
                pos = toMeasure_opt[i]
            end

            plotCyl_Gen(toMeasure_opt; B=(-20,20), L=(-20,70), r=1.2) # plot position measured
            [dict[i] = [data[i], 1] for i in toMeasure_opt] # collecting simulation data
        end

        TT = 2*1.29*nc*NI + 0.002*NS + 15*NI # find total time (seconds)
        println("    new individuals: ", NI)
        println("    total steps (stepper motor): ", NS)
        println("    simulation time: ", round(TT/60,digits=2), " min")

        old_pop = new_pop
        Tm = Tmmax*sigmoid(ngen, offset=offset, α=α) # finding mutation rate
        new_pop = ga_onestep(old_pop, dict, Tr=Tr, Tm=Tm, size=size[1:2*nc], coord_zero=coord_zero[1:2*nc]) # run genetic algorithm

        criteriaConv = criteria(new_pop, old_pop, μ=μ, σ=σ, Δμ=Δμ, Δσ=Δσ) # evaluate criteria
        ngen += 1 # move to next generation
        println("    criteria: ", criteriaConv)

    end

    return [new_pop, dict, TT, NI, NS, criteriaConv]
end

# Determining all the parameters for one cylinder case study
# ga_simulator(;n=30, nc=3, Tr=0.25, Tmmax=0.01, offset=25, α=0.4,
#    μ=1.0, σ=0.5, Δμ=1.0, Δσ=1.0, ngen_max=100)

function params_onecylinder(ncases=100; maxval= -62.74781119, nrep = 200, ngen_max=1000)

    Tmean = []; Tsd = []; V = []; C = []; params = []

    for i in 1:ncases
        pop_size = rand(10:30)
        tr_range = rand(0.1:0.05:0.6)
        tm_range = rand(0.002:0.002:0.03)
        sigmoid_offset = rand(10:50)
        sigmoid_alpha = rand(0.1:0.05:0.9)
        convergence_mu = rand(0.2:0.1:2.0)
        convergence_sigma = rand(0.1:0.05:0.9)

        push!(params, [pop_size,tr_range,tm_range,sigmoid_offset,sigmoid_alpha,convergence_mu,convergence_sigma])

        println("Iteration ", i)

        Ttemp = []; Vtemp = []
        nconv = 0
        for _ in 1:nrep
            _,d,t,conv = ga_simulator(n=pop_size, nc=1, Tr=tr_range, Tmmax=tm_range, offset=sigmoid_alpha,
                        α=sigmoid_alpha, μ=convergence_mu, σ=convergence_sigma, ngen_max=ngen_max)

            if conv
                push!(Ttemp, t/60)
                v = d[argmax(d)][1]
                diff_v = abs(v-maxval)
                push!(Vtemp, diff_v)
                nconv += 1
            end

        end

        println("   End of iteration")
        println("   Summary -> total: ", nrep, ", converged: ", nconv)

        push!(Tmean, mean(Ttemp))
        push!(Tsd, std(Ttemp))
        push!(V, mean(Vtemp))
        push!(C, nconv)
    end

    return [Tmean,Tsd,V,C,params]
end


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ 8c102f1e-e978-4ed3-8620-830c46bb566f
using Pkg; Pkg.develop(path=".."); Pkg.activate("..")

# ╔═╡ 0ae615d5-739f-456f-be76-d7d3307577bb
begin
	using LinearAlgebra
	using StatsBase
	using CSV
	using DataFrames
	using RDatasets
	using Plots
	using StatsPlots
	plotly(size=(680, 500))
end

# ╔═╡ 9bbf3d20-0b30-11ec-1917-ebba9c2ee5f0
md"
# Principal Component Analysis

JingYu Ning

**PCA** is a promissing **decorrelating dimension reduction.**
"

# ╔═╡ d67717f9-b5d6-4f84-ad80-31fe0c93b014
md"
## Let's start from a simple example

dataset: [Height and shoe size](https://osf.io/ja9dw/)

It is said that the length of ones foot is corelated with ones height...
"

# ╔═╡ a4d300fb-f712-4045-a8b8-7bf3b4b55ecc
md"
### Show data
"

# ╔═╡ 820f64a0-287d-40fa-b5f5-2cb35dc91658
SIG = "⚡";

# ╔═╡ 98d1834d-f61e-4b28-a6bc-9dec42a1bff8
begin
	SIG
	raw = CSV.read(joinpath(@__DIR__, "../data/HeightAndShoeSize.csv"), DataFrame)
	raw = raw[!, 2:end]
	data = copy(raw)
end

# ╔═╡ 53400761-545e-4f2a-ae75-80694a43dd92
@df data scatter(
	:height, :shoe_size, group=:sex,
	xlabel="height", ylabel="shoe size",
	legend=:topleft
)

# ╔═╡ 033440c0-efc2-42c4-918a-e0d966e27d3f
md"
### 1. Move the data to the origin
"

# ╔═╡ b3647ffa-d904-4bee-9327-6dec393702e5
begin
	for i in 2:3
		data[!, i] .-= mean(data[!, i])
		raw[!, i] .-= mean(raw[!, i])
	end
	@df data scatter(
		:height, :shoe_size, group=:sex,
		xlabel="height", ylabel="shoe size",
		legend=:topleft,
		xlim=(-40, 40), ylim=(-40, 40)
	)
end

# ╔═╡ 4f771fd7-84be-43f9-a9b2-78d40db4f91c
md"
### 2. SVD

$X = U \Sigma V^\dagger$
"

# ╔═╡ 6b6365a7-486a-4533-bd1e-d6c5726f3b35
u, s, v = svd(Matrix(data[!, 2:3])')

# ╔═╡ ddfc6838-cf88-4aa2-a95c-5d4c95a487ae
md"
### 3. Principal components

Choose the orthogonal basis `u` as the **principal components**
"

# ╔═╡ 66806d5b-55b7-4a4f-b592-2fe6304a7795
begin
	scale = 10
	p = @df data scatter(
		:height, :shoe_size, group=:sex,
		xlabel="height", ylabel="shoe size",
		legend=:topleft,
		xlim=(-40, 40), ylim=(-40, 40)
	)
	for i in 1:size(u, 1)
		pci = scale * u[:, i]
		plot!(p, [0, pci[1]], [0, pci[2]], label="pc$i", lw=5)
	end
	p
end

# ╔═╡ de3026a5-109e-458c-9e38-8e6ab43a81c6
md"
### 4. Apply the basis

$U^T X$
"

# ╔═╡ 5530ddd4-df5e-4308-9dfd-492f68db9f94
begin
	data[!, 2:3] .= (u' * Matrix(data[!, 2:3])')'
	@df data scatter(
		:height, :shoe_size, group=:sex,
		xlabel="pc1", ylabel="pc2",
		legend=:topleft,
		xlim=(-40, 40), ylim=(-40, 40)
	)
end

# ╔═╡ c7490656-5f5f-41a5-8a56-e3341880d855
md"
### 5. Rescale

$\Sigma^{-1} U^T X$
"

# ╔═╡ 89767b3c-63eb-4bdb-9f96-2289954908fa
begin
	data[!, 2:3] .= (inv(diagm(s)) * Matrix(data[!, 2:3])')'
	@df data scatter(
		:height, :shoe_size, group=:sex,
		xlim=(-0.5, 0.5), ylim=(-0.5, 0.5)
	)
end

# ╔═╡ 4ae5b880-1c74-4d37-b1eb-6f6f847b39a0
md"
$I = \Sigma^{-1} U^T X V$
"

# ╔═╡ 6a3f7794-3776-48cf-b188-ee892078c148
Matrix(data[!, 2:3])' * v

# ╔═╡ edcc9e1c-9cc0-4ead-ae6f-4c68bc509470
md"
### 6. Dimension reduction

We want to map those data into lower dimentions with less intormance loss.

Here, $U$ is a basis in the original high dimentional space. We will have least informance loss if project data into one of the axis of $U$, while have most loss if into another.

$U = [v_1, v_2]$
"

# ╔═╡ 18c0c325-c107-42e8-84a3-587580cc5772
md"
By choosing a new basis, we can therefore drop one dimention with little information lossed.
"

# ╔═╡ 88526e35-67c0-430a-a6ef-6454cdab7d2d
begin
	plot(size=(680, 200), yticks=[])
	scatter!((u'*Matrix(raw[!, 2:3])')[1, :], zeros(nrow(raw)), label="pc1")
	scatter!((u'*Matrix(raw[!, 2:3])')[2, :], ones(nrow(raw)), label="pc2")
end

# ╔═╡ 4d5c2e27-8a2c-486b-97bd-65ce4d108322
md"
And the singular value is the information of varience.
"

# ╔═╡ 5312759e-6c75-411b-b28f-9de6f6bfa65b
begin
	plot(size=(680, 200), yticks=[])
	scatter!(
		(inv(diagm(s))*u'*Matrix(raw[!, 2:3])')[1, :], zeros(nrow(raw)),
		label="pc1"
	)
	scatter!(
		(inv(diagm(s))*u'*Matrix(raw[!, 2:3])')[2, :], ones(nrow(raw)),
		label="pc2"
	)
end

# ╔═╡ d5fe909b-784a-443e-9e2e-c00c9e8cd611
md"
We constract a mapping that maps the data into `Z` space

$W = \Sigma^{-1}U^T$
"

# ╔═╡ e3e19d9e-767a-486b-8626-108c5f70a63b
w = inv(diagm(s)) * u'

# ╔═╡ d315c4c8-9945-409a-90b5-f08369665b25
md"The mapping is roughly the overall picture of the data."

# ╔═╡ 9658209d-57b5-4d68-a031-e15711d8cc92
md"
## Iris dataset
"

# ╔═╡ 7475abbd-6fdb-423f-8981-be4461df6648
iris = dataset("datasets", "iris")

# ╔═╡ 663406ed-49d2-439b-8272-69c0b29546da
begin
    cols = names(iris)

    function gen_pots()
        plots = Matrix{Plots.Plot}(undef, 4, 4)
        for i in 1:4, j in 1:4
            if i == j
                plots[i, j] = @df(iris, histogram(cols(i)))
            elseif i < j
                plots[i, j] = @df(iris, histogram2d(
                    cols(i), cols(j),
                    nbins=20, color=:coolwarm, colorbar=false
                ))
            else
                plots[i, j] =  @df(iris, scatter(cols(i), cols(j)))
            end

            (i==1) && (plots[i, j] = plot!(ylabel=cols[j])) # add x label
            (j==4) && (plots[i, j] = plot!(xlabel=cols[i])) # add y label
        end

        return plots
    end
end

# ╔═╡ f79ca361-a2ca-4e0e-b920-36936d01fe8b
plot(gen_pots()..., layout=(4, 4), legend=false)

# ╔═╡ 7a284458-83bb-48fb-863e-ff385a3ba27d
md"
### 1. Move the data to the origin
"

# ╔═╡ da5d583f-b9d8-4c96-892c-c3a018d83890
begin
	for i in 1:(ncol(iris)-1)
		iris[!, i] .-= mean(iris[!, i])
	end
	plot(gen_pots()..., layout=(4, 4), legend=false)
end

# ╔═╡ 221d3ebf-e734-4a7d-8811-d8874311aa21
md"
### 2. SVD
"

# ╔═╡ 58db4367-514b-4ce8-b994-a13d73af93e9
U, S, V = svd(Matrix(iris[!, 1:(ncol(iris)-1)])')

# ╔═╡ 4a0d265d-9c91-4095-899b-20e4c6765fbf
md"
### 3. Construct the mapping

$W = \Sigma^{-1}U^T$
"

# ╔═╡ 1227c477-2fa5-4929-b64b-aeaf5745c52d
W = inv(diagm(S)) * U'

# ╔═╡ c1b885d9-bb2b-486f-ad22-219e4ede6025
begin
	iris[!, 1:(ncol(iris)-1)] .= (W * Matrix(iris[!, 1:(ncol(iris)-1)])')'
	plot(gen_pots()..., layout=(4, 4), legend=false)
end

# ╔═╡ Cell order:
# ╟─9bbf3d20-0b30-11ec-1917-ebba9c2ee5f0
# ╟─8c102f1e-e978-4ed3-8620-830c46bb566f
# ╟─0ae615d5-739f-456f-be76-d7d3307577bb
# ╟─d67717f9-b5d6-4f84-ad80-31fe0c93b014
# ╟─a4d300fb-f712-4045-a8b8-7bf3b4b55ecc
# ╟─98d1834d-f61e-4b28-a6bc-9dec42a1bff8
# ╟─53400761-545e-4f2a-ae75-80694a43dd92
# ╠═820f64a0-287d-40fa-b5f5-2cb35dc91658
# ╟─033440c0-efc2-42c4-918a-e0d966e27d3f
# ╟─b3647ffa-d904-4bee-9327-6dec393702e5
# ╟─4f771fd7-84be-43f9-a9b2-78d40db4f91c
# ╟─6b6365a7-486a-4533-bd1e-d6c5726f3b35
# ╟─ddfc6838-cf88-4aa2-a95c-5d4c95a487ae
# ╟─66806d5b-55b7-4a4f-b592-2fe6304a7795
# ╟─de3026a5-109e-458c-9e38-8e6ab43a81c6
# ╟─5530ddd4-df5e-4308-9dfd-492f68db9f94
# ╟─c7490656-5f5f-41a5-8a56-e3341880d855
# ╟─89767b3c-63eb-4bdb-9f96-2289954908fa
# ╟─4ae5b880-1c74-4d37-b1eb-6f6f847b39a0
# ╟─6a3f7794-3776-48cf-b188-ee892078c148
# ╟─edcc9e1c-9cc0-4ead-ae6f-4c68bc509470
# ╟─18c0c325-c107-42e8-84a3-587580cc5772
# ╟─88526e35-67c0-430a-a6ef-6454cdab7d2d
# ╟─4d5c2e27-8a2c-486b-97bd-65ce4d108322
# ╟─5312759e-6c75-411b-b28f-9de6f6bfa65b
# ╟─d5fe909b-784a-443e-9e2e-c00c9e8cd611
# ╟─e3e19d9e-767a-486b-8626-108c5f70a63b
# ╟─d315c4c8-9945-409a-90b5-f08369665b25
# ╟─9658209d-57b5-4d68-a031-e15711d8cc92
# ╟─663406ed-49d2-439b-8272-69c0b29546da
# ╠═7475abbd-6fdb-423f-8981-be4461df6648
# ╟─f79ca361-a2ca-4e0e-b920-36936d01fe8b
# ╟─7a284458-83bb-48fb-863e-ff385a3ba27d
# ╟─da5d583f-b9d8-4c96-892c-c3a018d83890
# ╟─221d3ebf-e734-4a7d-8811-d8874311aa21
# ╟─58db4367-514b-4ce8-b994-a13d73af93e9
# ╟─4a0d265d-9c91-4095-899b-20e4c6765fbf
# ╟─1227c477-2fa5-4929-b64b-aeaf5745c52d
# ╟─c1b885d9-bb2b-486f-ad22-219e4ede6025


text:
using Test
using ARMA
using LinearAlgebra
using Pope

# Tests:
#
# Verify projectors*basis == eye(...)
# Verify both methods give same projectors if noisecovariance is that of the ARMA model.
# Verify projectors and pcovariance scale properly if you scale the basis.
# Verify projectors are the same if you scale the noise model, but pcovariance changes.
# Error if noisecovariance argument is too short.
# Error if basis columns aren't independent.
# Error if noise is zero.
# Error if basis has more columns than rows.

@testset "projections" begin

# Generate a test basis of fixed size, a trivial noise model, and two versions
# of a non-trivial noise model.
(N, n) = 100, 6
basis = randn(N, n)
white = zeros(Float64, N); white[1] = 1.0

model = ARMA.ARMAModel([2.41644, -4.25252, 1.85175], [1, -1.90467, .905899])
covar = ARMA.model_covariance(model, N)

# Compute projections for all three noise models
(wproj,pwcovar) = Pope.computeprojectors(basis, white)
(mproj,pmcovar) = Pope.computeprojectors(basis, model)
(cproj,pccovar) = Pope.computeprojectors(basis, covar)

@testset "core" begin
    # Verify that projecting the basis (for any noise model) yields identity.
    @test wproj*basis ≈ Matrix(1.0I, n, n)
    @test mproj*basis ≈ Matrix(1.0I, n, n)
    @test cproj*basis ≈ Matrix(1.0I, n, n)

    # Verify projecting combinations of the basis works, too.
    # Yes, I realized that success here is implied by the eye(n) tests, above.
    @test cproj*2basis ≈ 2Matrix(1.0I, n, n)
    truecoef = randn(n,n)
    b = basis*truecoef
    @test truecoef ≈ cproj*b

    # Verify both methods give same projectors if noisecovariance is that of the ARMA model.
    @test mproj ≈ cproj

end

@testset "scaling" begin
    # Verify projectors and pcovariance scale properly if you scale the basis.
    (cproj2,pccovar2) = Pope.computeprojectors(2basis, covar)
    @test cproj ≈ 2cproj2
    @test pccovar ≈ 4pccovar2

    # Verify projectors are the same if you scale the noise model, but pcovariance changes.
    (cproj3,pccovar3) = Pope.computeprojectors(basis, 3covar)
    @test cproj ≈ cproj3
    @test 3pccovar ≈ pccovar3
end

@testset "errors" begin
    # Noise input is too short
    @test_throws ArgumentError Pope.computeprojectors(basis, white[1:N-1])

    # Basis vectors are degenerate
    degeneratebasis = copy(basis)
    degeneratebasis[:,end] = sum(degeneratebasis[:,1:3], dims=2)
    @test_throws ArgumentError Pope.computeprojectors(degeneratebasis, white)
    degeneratebasis[:,end] .= 0.0
    @test_throws ArgumentError Pope.computeprojectors(degeneratebasis, white)

    # Noise is zero
    @test_throws SingularException Pope.computeprojectors(basis, 0*white)

    # Basis has more columns than rows
    fatbasis = randn(6, 10)
    @test_throws ArgumentError Pope.computeprojectors(fatbasis, white)

    # @test_throws ErrorException Pope.computeprojectors(fatbasis, white)
end

end


text:
"""
    centrality_katz(N::Unipartite; a::Float64=0.1, k::Int64=5)

This measure can work on different path length (`k`), and give a different
weight to every subsequent connection (`a`). `k` must be at least 1 (only
immediate neighbors are considered). `a` (being a weight), must be positive.

> Katz, L., 1953. A new status index derived from sociometric analysis.
> Psychometrika 18, 39–43. doi:10.1007/bf02289026

"""
function centrality_katz(N::Union{UnipartiteNetwork, UnipartiteProbabilisticNetwork}; a::Float64=0.1, k::Int64=5)
	@assert 0.0 <= a <= 1.0
  @assert k >= 1

	centr = sum(hcat(map((x) -> vec(sum((a^x).*(N.edges^x); dims=1)), [1:k;])...); dims=2)
	return Dict(zip(species(N), centr ./ sum(centr)))
end

"""
    centrality_closeness(N::UnipartiteNetwork; nmax::Int64=20)

The function calls `shortest_path` internally -- the `nmax` argument is the
maximal path length that will be tried.

#### References

- Bavelas, A., 1950. Communication Patterns in Task‐Oriented Groups. The Journal
  of the Acoustical Society of America 22, 725–730. doi:10.1121/1.1906679
"""
function centrality_closeness(N::UnipartiteNetwork; nmax::Int64=20)
  d = shortest_path(N, nmax=nmax)
  n = richness(N)-1
  d[diagind(d)] .= 0
  interm = sum(d; dims=2)
  interm = vec(n ./ interm)
  for i in eachindex(interm)
    interm[i] = interm[i] == Inf ? 0.0 : interm[i]
  end
  return Dict(zip(species(N), interm))
end

"""
    centrality_harmonic(N::UnipartiteNetwork; nmax::Int64=20)

The function calls `shortest_path` internally -- the `nmax` argument is the
maximal path length that will be tried.
"""
function centrality_harmonic(N::UnipartiteNetwork; nmax::Int64=20)
  d = shortest_path(N, nmax=nmax)
  n = richness(N)-1
  d[diagind(d)] .= 0
  interm = 1.0 ./ d
  for i in eachindex(interm)
    interm[i] = isinf(interm[i]) ? 0.0 : interm[i]
  end
  interm = sum(d; dims=2)./n
  return Dict(zip(species(N), interm))
end

"""
    centrality_degree(N::UnipartiteNetwork)

Degree centrality, corrected by the maximum degree (the most central species has
a degree of 1).
"""
function centrality_degree(N::UnipartiteNetwork)
  d = degree(N)
  dm = maximum(values(d))
  return Dict([p.first=>p.second/dm for p in d])
end

"""
    centrality_eigenvector(N::AbstractUnipartiteNetwork)

Eigen centrality, corrected by the maximum degree (the most central species has
a degree of 1).
"""
function centrality_eigenvector(N::AbstractUnipartiteNetwork)
  evals, evecs = eigen(Array(N.edges))
  emax, epos = findmax(real.(evals))
  cvals = vec(real.(evecs[:,epos]))
  cvals ./= maximum(cvals)
  return Dict(zip(species(N), cvals))
end


text:
"""
Model is an abstract type whose subtypes should provide the following:

* `train(model, data, loss)`
* `test(model, data, loss)`
* `predict(model, data)`
* `setopt!(model; param...)`
"""
abstract Model

setopt!(m::Model; o...)=(for p in params(m); setopt!(p; o...); end)
update!(m::Model; o...)=(for p in params(m); update!(p; o...); end)             # t:19
vnorm(x)=(x==nothing ? 0 : vecnorm(x))
wnorm(m::Model,w=0)=(for p in params(m); w += vnorm(p.out); end; w)           # t:317
gnorm(m::Model,g=0)=(for p in params(m); g += vnorm(p.dif); end; g)           # t:332

# Helper function for train/test
item2xy(item)=(isa(item, Tuple) ? (item[1:end-1],item[end]) : item==nothing ? (nothing,nothing) : ((),item))

# So gradient checking does not mess up random seed:
const gradcheck_rng = MersenneTwister()

# train(m, d; seq=false, a...)=(!seq ? train1(m,d;a...) : train2(m,d;a...))

# function train1(m::Model, data; loss=quadloss, gclip=0, gcheck=0, getnorm=true, getloss=true, a...)
#     numloss = sumloss = maxwnorm = maxgnorm = w = g = 0
#     for item in data
#         (x,ygold) = item2xy(item)
#         gcheck > 0 && (gradcheck(m,ygold,x...; gcheck=gcheck, loss=loss, a...); gcheck=0)
#         ypred = forw(m, x...; trn=true, a...)
#         getloss && (sumloss += loss(ypred, ygold); numloss += 1)
#         back(m, ygold; loss=loss, a...)
#         (getnorm || gclip>0) && (g = gnorm(m); g > maxgnorm && (maxgnorm = g))
#         update!(m; gclip=(g > gclip > 0 ? gclip/g : 0), a...)
#         getnorm && (w = wnorm(m); w > maxwnorm && (maxwnorm = w))
#     end
#     return (sumloss/numloss, maxwnorm, maxgnorm)
# end

# This will not work for MLP!  extra parameterless ops do not effect equality.
# Base.isequal(a::Model,b::Model)=(typeof(a)==typeof(b) && isequal(params(a),params(b)))

# function inittrain(m::Model, d)
#     isempty(params(m)[1]) || return
#     (x,n) = next(d,start(d))
#     init(m, x[1]; trn=true)
# end

# NO: make the model interface more functional:
# back and loss rely on hidden state info.
# forw has to allocate.
# purely functional models are impossible.
# forw needs to compute intermediate values.
# but from user's perspective forw is functional.
# loss/back is not: relying on history.
# we could give them x/y but they would still need internal state.
# if they are going to use internal state they may as well use the one set by forw.

# TODO: this does not work, cannot write back on data
# function predict(m::Model, d)
#     for (x,y) in d
#         forw(m, x; yout=y, trn=false)
#     end
# end

# function backprop(m::Model, x, y; getloss=true, a...)
#     loss1 = forw(m, x; trn=true, ygold=(getloss ? y : nothing), a...)
#     back(m, y; a...)
#     return loss1
# end

# Use test with percloss instead:
#
# function accuracy(m::Model, d) # TODO: this only works if y is a single item
#     numcorr = numinst = 0
#     z = nothing
#     for (x,y) in d
#         z == nothing && (z = Array(eltype(y), size(y)))
#         forw(m, x, z; mode=:test)
#         numinst += ccount(y)
#         numcorr += sum(findmax(convert(Array,y),1)[2] .== findmax(convert(Array,z),1)[2])
#     end
#     return numcorr/numinst
# end

# function train1(m::Model, data; loss=quadloss, getloss=true, getnorm=true, gclip=0, gcheck=0, seq=false, a...)
#     numloss = sumloss = maxwnorm = maxgnorm = w = g = 0
#     for item in data
#         (x,ygold) = item2xy(item)
#         ypred = forw(m, x...; trn=true, seq=false, a...)
#         l = getloss ? loss(ypred, ygold) : 0
#         sumloss += l; numloss += 1
#         back(m, ygold; loss=loss, seq=false, a...)
#         getnorm && (w = wnorm(m); w > maxwnorm && (maxwnorm = w))
#         (getnorm || gclip>0) && (g = gnorm(m); g > maxgnorm && (maxgnorm = g))
#         update!(m; gclip=(g > gclip > 0 ? gclip/g : 0))
#     end
#     return (sumloss/numloss, maxwnorm, maxgnorm)
# end

# function test(m::Model, d; o...)
#     sumloss = numloss = 0
#     for (x,y) in d
#         sumloss += forw(m, x; mode=:test, ygold=y, o...)
#         numloss += 1
#     end
#     return sumloss/numloss
# end


text:
export totalordermapcomponent, totalordermap

"""
$(TYPEDSIGNATURES)

Creates an `HermitemapComponent` of order `order` and dimension `Nx` for the basis `b`.
The features of the created maps are all the tensorial products of the basis elements up to the order `order`.
"""
function totalordermapcomponent(Nx::Int64, order::Int64; withconstant::Bool = false, b::String = "CstProHermiteBasis")
    @assert order >= 0 "Order should be positive"
    if b ∈ ["CstProHermiteBasis"; "CstPhyHermiteBasis"]
        B = MultiBasis(eval(Symbol(b))(order+2), Nx)
    elseif b ∈ ["CstLinProHermiteBasis"; "CstLinPhyHermiteBasis"]
        B = MultiBasis(eval(Symbol(b))(order+3), Nx)
    else
        error("Undefined basis")
    end

    idx = totalorder(order*ones(Int64, Nx))

    if withconstant == false
        idx = idx[2:end,:]
    end

    Nψ = size(idx, 1)

    f = ExpandedFunction(B, idx, zeros(Nψ))
    return HermiteMapComponent(IntegratedFunction(f))
end

"""
$(TYPEDSIGNATURES)

Creates an `Hermitemap` of order `order` for the basis `b`. The size of the map is extracted from the ensemble matrix `X`.
The features of the created maps are all the tensorial products of the basis elements up to the order `order`.
"""
function totalordermap(X::Array{Float64,2}, order::Int64; diag::Bool=true, factor::Float64 = 1.0, withconstant::Bool = false, b::String = "CstProHermite")

    L = LinearTransform(X; diag = diag, factor = factor)

    Nx = size(X,1)
    C = HermiteMapComponent[]
    @inbounds for i=1:Nx
        push!(C, totalordermapcomponent(i, order; withconstant = withconstant, b = b))
    end

    if b ∈ ["CstProHermiteBasis"; "CstPhyHermiteBasis"]
        m = order+2
    elseif b ∈ ["CstLinProHermiteBasis"; "CstLinPhyHermiteBasis"]
        m = order+3
    else
        error("Undefined basis")
    end

    return HermiteMap(m, Nx, L, C)
end


text:
module NonlinearSpectralModels

# package code goes here
using Distributions

abstract NonlinearSpectralModel
abstract TruneOptions
abstract SpectralDataField

export NonlinearSpectralModel, TruneOptions, SpectralDataField
include("helpers/advections.jl")
type TruneSS <: TruneOptions
  n_off :: Int64
end

type TruneAllS <: TruneOptions
  n_off :: Int64
end

type TruneSV <: TruneOptions
  n_off :: Int64
end

immutable TruneNone <: TruneOptions end

export TruneSS, TruneNone,  TruneAllS, TruneSV, TruneOptions

function timegradient(  m:: NonlinearSpectralModel, x :: SpectralDataField ,TR :: TruneOptions)
  error("No function for abstract types")
end

import Base.+
import Base.-
import Base.*
import Base./

function RK4(m:: NonlinearSpectralModel, x :: SpectralDataField ,TR :: TruneOptions, dt :: Float64)
  y1 = timegradient(m,x,TR)
  y2 = timegradient(m,x+0.5*y1*dt,TR)
  y3 = timegradient(m,x+0.5*y2*dt,TR)
  y4 = timegradient(m,x+0.5*y3*dt,TR)
  return x+dt*(y1+2.0*y2+2.0*y3+y4)*(1/6.0)
end

include("Burger/model.jl")
export timegradient, RK4
end # module


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()
##

using DataFrames
using CSV
using Plots
using StatsPlots
using Plots.PlotMeasures
using Statistics
using StatsFuns
##

fpairs = CSV.read("../data/fpairs.txt", DataFrame, header=false)[:,1]

##

theme(:solarized_light)
#
upscale = 2 #8x upscaling in resolution
fntsm = Plots.font("sans-serif", pointsize=round(8.0*upscale))
fntlg = Plots.font("sans-serif", pointsize=round(10.0*upscale))
default(titlefont=fntlg, guidefont=fntlg, tickfont=fntsm, legendfont=fntsm)
default(size=(400*upscale,600*upscale)) #Plot canvas size

##

bs_ = []
for i in 1:28
    push!(bs_, CSV.read("bsResults/$(lpad(i,2,"0")).csv", DataFrame, header=false))
end

bs = vcat(bs_...)
insertcols!(bs,1,:fpair=>fpairs)
rename!(bs, [:fpair, :universal, :lineage])

bs[:,:bf] = bs.universal .- bs.lineage

sort!(bs, :bf)

##

scatter(bs.bf, 1:28,
    legend = false,
    yticks= (1:28,bs.fpair),
    xlab="(log) Bayes factor",
    left_margin = upscale*5mm,
    markersize=upscale*3,
    )

savefig("../data/img/bayesfactor.pdf")

##

sort!(bs, :bf, rev=true)

p0 = logistic.(bs.bf)

bs[:,:cum_pv] = cumsum(1 .- p0)

output = DataFrame(
    "feature pair" => bs.fpair,
    "(log) Bayes Factor" => round.(bs.bf, digits=1),
    "cumulative posterior probability" => rpad.(round.(bs.cum_pv, sigdigits=3), 5, "0"),
)

CSV.write("../data/tables/bf.tex", output, delim="&", newline="\\\\\n")

##

max_stars_repo_name: erathorn/phylogeneticTypology

max_issues_repo_name: erathorn/phylogeneticTypology

max_forks_repo_name: erathorn/phylogeneticTypology


text:
#!/usr/bin/julia

# Copyright (C) 2021, 2022 Matheus Fernandes Bigolin <mfrdrbigolin@disroot.org>
# SPDX-License-Identifier: MIT

"Day One, Sonar Sweep."

include("Utils.jl")

using .Utils

"""Given the `lst` and the `step_size`, calculate the sums which are larger
than previous sums."""
function solve(lst :: Vector{Int}, step_size = 1)
    let Σ = 0
        for (a, b) in zip(lst[1 : end - step_size], lst[step_size + 1 : end])
            Σ += a < b
        end

        return Σ
    end
end

function main()
    usage_and_exit(length(ARGS) != 1)

    input_data = arrange(readlines(ARGS[1]), Int)

    println(solve(input_data))
    println(solve(input_data, 3))
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
@testset "proximal maps" begin
  #
  # proxTV
  p = SnPoint([1.,0.,0.])
  q = SnPoint([0.,1.,0.])
  M = Sphere(2)
  N = Power(M,(2,))
  @test_throws ErrorException proxDistance(M, 1., p, q, 3)
  @test distance(M, proxDistance(M, distance(M,p,q)/2, p,q,1), midPoint(M,p,q)) ≈ 0
  (r,s) = proxTV(M,π/4,(p,q))
  @test norm( getValue(r) - getValue(s) ) < eps(Float64)
  # i.e. they are moved together
  @test distance(M,r,s) < eps(Float64)
  (t,u) = proxTV(M,π/8,(p,q));
  @test_throws ErrorException proxTV(M, π, (p,q), 3)
  # they cross correlate
  @test ( abs(t.value[1]-u.value[2])< eps(Float64) && abs(t.value[2]-u.value[1]) < eps(Float64) && abs(t.value[3]-u.value[3])< eps(Float64) )
  @test distance(M,t,u) == π/4 # and have moved half their distance
  #
  (v,w) = proxTV(M,1.,(p,q),2)
  vC, wC = geodesic(M, p, q, [1/3, 2/3])
  @test distance(M, v, vC) ≈ 0
  @test distance(M, w, wC) ≈ 0
  # proxTV on Power
  T = proxTV(N, π/8, PowPoint([p,q]))
  @test distance(N, T, PowPoint([t, u])) ≈ 0
  # parallelproxTV
  N2 = Power(M,(3,))
  r = midPoint(M,p,q)
  s, t = proxTV(M, π/16, (r, q) )
  u, v = proxTV(M, π/16, (p, r) )
  y = proxParallelTV(N2, π/16, [PowPoint([p,r,q]), PowPoint([p,r,q])])
  @test distance(N2, y[1], PowPoint([p,s,t])) ≈ 0 # even indices in first comp
  @test distance(N2, y[2], PowPoint([u,v,q])) ≈ 0 # odd in second
  # dimensions of x have to fit, here they don't
  @test_throws ErrorException proxParallelTV(N2, π/16, [PowPoint([p,r,q])])
  # proxTV2
  p2, r2, q2 = proxTV2(M,1., (p,r,q) )
  sum(distance.(Ref(M), [p,r,q], [p2, r2, q2] )) ≈ 0
  @test_throws ErrorException proxTV2(M, 1., (p,r,q), 2) # since proxTV is only defined for p=1
  distance(Power(M,3), PowPoint([p2,r2,q2]), proxTV2(Power(M,3),1., PowPoint([p,r,q]))) ≈ 0
  # Circle
  M2 = Circle()
  pS, rS, qS = S1Point.([-0.5, 0.1, 0.5])
  d = dot( getValue.([pS, rS, qS]), [1., -2., 1.] )
  m = min(0.3, abs( symRem(d)/6) )
  s = sign(symRem(d))
  pSc, rSc, qSc = S1Point.( symRem.( getValue.([pS, rS, qS]) .- m .* s .* [1., -2., 1.] ) )
  pSr, rSr, qSr = proxTV2(M2, 0.3, (pS, rS, qS) )
  @test sum( distance.(Ref(M2), [pSc, rSc, qSc], [pSr, rSr, qSr]) ) ≈ 0
  # p=2
  t = 0.3*symRem(d)/(1+0.3*6.)
  @test sum(
    distance.( Ref(M2),
      [proxTV2(M2,0.3,(pS, rS, qS),2)...],
      S1Point.( getValue.([pS, rS, qS]) .- t.*[1., -2., 1.])
  )) ≈ 0
  # others fail
  @test_throws ErrorException proxTV2(M2, 0.3, (pS, rS, qS),3)
  # Rn
  M3 = Euclidean(1)
  pR, rR, qR = RnPoint.( getValue.([pS, rS, qS]) )
  m = min.( Ref(0.3), abs.( getValue.([pR, rR, qR]).*[1., -2., 1.] )/6 )
  s = sign(d)  # we can reuse d
  pRc, rRc, qRc = RnPoint.( getValue.([pR, rR, qR])  .- m .* s .* [1., -2., 1.] )
  pRr, rRr, qRr = proxTV2(M3, 0.3, (pR, rR, qR) )
  @test sum( distance.(Ref(M3), [pRc, rRc, qRc], [pRr, rRr, qRr]) ) ≈ 0
  # p=2
  t = 0.3*d/(1+0.3*6.)
  @test sum(
    distance.( Ref(M3),
      [proxTV2(M3, 0.3, (pR, rR, qR), 2)...],
      RnPoint.( getValue.([pR, rR, qR]) .- t.*[1., -2., 1.])
  )) ≈ 0
  # others fail
  @test_throws ErrorException proxTV2(M3, 0.3, (pR, rR, qR),3)
  #
  # collaborative integer tests
  #
  ξR, ηR, νR = RnTVector.( getValue.([pS, rS, qS]) )
  N3 = Power(M3,3)
  P = PowPoint([pR, rR, qR])
  Ξ = PowTVector([ξR, ηR, νR])
  @test proxCollaborativeTV(N3,0.,P, Ξ,1, 1) == Ξ
  @test proxCollaborativeTV(N3,0.,P, Ξ,1., 1) == Ξ
  @test proxCollaborativeTV(N3,0.,P, Ξ,1, 1.) == Ξ
  @test proxCollaborativeTV(N3,0.,P, Ξ,1., 1.) == Ξ

  @test proxCollaborativeTV(N3,0.,P, Ξ, 2, 1) == Ξ
  @test norm(N3,P, proxCollaborativeTV(N3,0.,P, Ξ, 2, Inf) - zeroTVector(N3,P)) ≈ 0
  @test norm(N3,P, proxCollaborativeTV(N3,0.,P, Ξ, 1, Inf) - zeroTVector(N3,P)) ≈ 0
  @test norm(N3,P, proxCollaborativeTV(N3,0.,P, Ξ, Inf, Inf) - zeroTVector(N3,P)) ≈ 0
  @test_throws ErrorException proxCollaborativeTV(N3,0.,P, Ξ, 3,3)
  @test_throws ErrorException proxCollaborativeTV(N3,0.,P, Ξ, 3,1)
  @test_throws ErrorException proxCollaborativeTV(N3,0.,P, Ξ, 3,Inf)
end


text:
# CPU cross-correlation using files split into Npersplit windows per file
# run splitseis.jl with Npersplit set to 25 before running this
# this assumes a 48-core machine

using Distributed
addprocs(48,topology=:master_worker)
@everywhere begin
using SeisIO, SeisNoise, Dates, Glob
fs = 1.
cc_len, cc_step = 2^15, 26768
freqmin = 0.003
freqmax = 0.005
maxlag = 12000.
Npersplit = 25
splits = glob("*",joinpath(expanduser("~/SeisNoise_paper/DATA/SPLITSEIS"),string(Npersplit)))
Nsplits = length(splits)
SPLITCORR = joinpath(expanduser("~/SeisNoise_paper/CORR/SPLITSEIS"),string(Npersplit))
XML = expanduser("~/SeisNoise_paper/DATA/XML")
N = length(splits)

if !isdir(SPLITCORR)
    mkpath(SPLITCORR)
end

# this function reads a file, applies pre-processing, converts to RawData
#  and then returns an FFTData
function seis2fft(file::String,XML::String,cc_len::Int,cc_step::Int,
                  freqmin::Float64,freqmax::Float64)
    S = SeisData()
    read_data!(S,"sac",file)
    detrend!(S)
    taper!(S)
    filtfilt!(S,fl=0.001,fh=0.4,np=3)

    # get xmlfile
    net,sta,loc,chan = split(S[1].id,'.')
    xmlfile = joinpath(XML,net*'.'*sta*".xml")
    Resp = read_meta("sxml",xmlfile)
    translate_resp!(S,Resp)
    remove_resp!(S)
    filtfilt!(S,fl=freqmin,fh=freqmax,np=3)
    R = RawData(S[1],cc_len,cc_step)
    demean!(R)
    taper!(R)
    return rfft(R)
end

# this function cross-correlates data for a "split"
function splitcorr(SPLITDIR::String,XML::String,cc_len::Int,cc_step::Int,
                  freqmin::Float64,freqmax::Float64,maxlag::Float64,OUTDIR::String)
    files = glob("*",SPLITDIR)
    Nfiles = length(files)
    println("Cross-correlating chunk $(basename(SPLITDIR)) ",now())
    SPLITOUT = joinpath(OUTDIR,basename(SPLITDIR))
    if !isdir(SPLITOUT)
        mkpath(SPLITOUT)
    end

    # I/O + fft
    Fs = Array{FFTData}(undef,Nfiles)
    for ii = 1:Nfiles
        Fs[ii] = seis2fft(files[ii],XML,cc_len,cc_step,freqmin,freqmax)
    end

    # cross-correlate
    for ii = 1:Nfiles -1
        for jj = ii+1:Nfiles
            C = compute_cc(Fs[ii],Fs[jj],maxlag)
            stack!(C,allstack=true)
            save_corr(C,SPLITOUT)
        end
    end
    return nothing

end
end

# precompile
pmap(splitcorr,splits[1:1],[XML],[cc_len],[cc_step],
                  [freqmin],[freqmax],[maxlag],[SPLITCORR])
toremove = glob("*/*",SPLITCORR)
rm.(toremove)

println("Beginning XCORR with $(nprocs()-1) processors...")
tALL = @elapsed begin
pmap(splitcorr,splits,fill(XML,Nsplits),fill(cc_len,Nsplits),fill(cc_step,Nsplits),
     fill(freqmin,Nsplits),fill(freqmax,Nsplits),fill(maxlag,Nsplits),fill(SPLITCORR,Nsplits))
end

println("XCORR took $tALL seconds")


text:
"""
    vtk_write_array(filename, arrays, labels)
    vtk_write_array(filename, array; label = "array")

Write Julia arrays to a VTK image data file (.vti).

Useful for general visualisation of arrays.
The input can be a 2D or 3D array.

Multiple arrays can be given as a tuple.
For instance,

    vtk_write_array(filename, (u, v), ("u", "v"))

In that case, the arrays must have the same dimensions.

"""
function vtk_write_array(filename::AbstractString,
                         arrays::NTuple{M, A},
                         labels::NTuple{M, S}) where
        {T <: Real, M, N, A <: AbstractArray{T,N}, S <: AbstractString}
    if !(2 <= N <= 3)
        throw(DimensionMismatch("input should be a 2D or 3D array."))
    end
    vtk_grid(filename, size(arrays[1])...) do vtkfile
        for (array, label) in Iterators.zip(arrays, labels)
            vtkfile[label, VTKPointData()] = array
        end
    end
end

vtk_write_array(filename::AbstractString,
                array::AbstractArray{T},
                label::AbstractString = "array") where T <: Real =
vtk_write_array(filename, (array, ), (label, ))


text:
module DynamicsPlots

using Plots, RecipesBase
using DataFrames

include("hiscatter.jl")
include("timeevolution.jl")
include("pointcloud.jl")
include("geometry.jl")

export hiscatter, timeevolution, pointcloud, geometry, circle!

end


text:
module SentryClients
using gRPCClient
using ProtoBuf:google

include("dapr.jl")
using .dapr.proto.sentry.v1

import Base: show

# begin service: dapr.proto.sentry.v1.CA

export CABlockingClient, CAClient

struct CABlockingClient
    controller::gRPCController
    channel::gRPCChannel
    stub::CABlockingStub

    function CABlockingClient(api_base_url::String; kwargs...)
        controller = gRPCController(; kwargs...)
        channel = gRPCChannel(api_base_url)
        stub = CABlockingStub(channel)
        new(controller, channel, stub)
    end
end

struct CAClient
    controller::gRPCController
    channel::gRPCChannel
    stub::CAStub

    function CAClient(api_base_url::String; kwargs...)
        controller = gRPCController(; kwargs...)
        channel = gRPCChannel(api_base_url)
        stub = CAStub(channel)
        new(controller, channel, stub)
    end
end

show(io::IO, client::CABlockingClient) = print(io, "CABlockingClient(", client.channel.baseurl, ")")
show(io::IO, client::CAClient) = print(io, "CAClient(", client.channel.baseurl, ")")

import .dapr.proto.sentry.v1: SignCertificate
"""
    SignCertificate

- input: dapr.proto.sentry.v1.SignCertificateRequest
- output: dapr.proto.sentry.v1.SignCertificateResponse
"""
SignCertificate(client::CABlockingClient, inp::dapr.proto.sentry.v1.SignCertificateRequest) = SignCertificate(client.stub, client.controller, inp)
SignCertificate(client::CAClient, inp::dapr.proto.sentry.v1.SignCertificateRequest, done::Function) = SignCertificate(client.stub, client.controller, inp, done)

# end service: dapr.proto.sentry.v1.CA

end # module DaprClients


text:
include("header.jl")

#=
# Dimensionless quantities
In this notebook, we will learn about making quantities dimensionless,
and the physical significance of this.
=#

# ### Set up the module
using MAE103
#-
using Plots

#=
In the previous notebook, we made a plot of the velocity field in a gap
between two walls. We'll regenerate the plot here for reference:
=#
u(y,Uc,H) = 4*Uc/H^2*y*(H-y)
H = 1u"cm"
Uc = 1u"m/s"
y = range(0u"cm",H,length=101)
v = u.(y,Uc,H)
plot(v,y,xlim=(0,Inf),ylim=(0,Inf),xlabel=L"u(y)",ylabel=L"y")

#=
In this plot, the data have units---they are *dimensional* data.
This is useful, since the plot not only conveys the velocities, but
also the gap size and the flow speed.

Let us also note that the plot is quite specific, since it only shows
the plot for *one* choice of gap size (1 cm) and *one* choice of flow speed (1 m/s).
Of course, it is not difficult for us to get plots at different gap sizes
and flow speeds by simply changing the parameters.

But suppose that these data were measured in an experiment. For the experiment,
let's say we constructed one gap $H$, we applied one flow speed $U_c$, and
we measured the resulting profile of fluid velocity across the gap. Then this plot,
at these specific conditions, would not be directly useful to someone else who would like
to know the velocity profile in a different size gap at a different flow speed.

If, instead, we divide the velocity data `v` by flow speed `Uc` and divide the position
data `y` by gap size `H`, then we "scale out" the specific choices for these parameters:
=#
plot(v/Uc,y/H,xlim=(0,Inf),ylim=(0,Inf),xlabel=L"u(y)/U_c",ylabel=L"y/H",yticks=0:0.25:1)

#=
This is a *dimensionless plot*. The values on each axis have no units: they are
called *dimensionless* or *non-dimensional* quantities.

This plot is more directly useful to another person, because they can read off
the values for a particular position *as a fraction of gap height*,
and find the velocity *as a fraction of the centerline velocity*. Then, they can
"scale" the value to their own situation.

For example, if a person whose gap size is $H = 3$ cm and centerline velocity is
$0.5$ m/s wishes to know the velocity at $y = 0.6$ cm, they calculate $y/H = 0.6/3 = 0.2$
and look up on the plot that $y/H = 0.2$ corresponds to $u/U_c = 0.64$. Thus,
for this person's case, $u = 0.64\times0.5$ m/s  $= 0.32$ m/s.
=#

#-
#=
### Unit as a reference quantity.
The principle of "scaling out" the parameters is very important. But before we
go deeper into this, let us make an observation on units.

A unit is a reference value, and this reference is standardized so we can communicate
results with one another on a common foundation. What does it mean to be a reference?

When we measure something, it is important to remember that the thing we're measuring is not
intrinsically tied to any unit. It is more abstract than that. For example, the
distance between two points is a physical fact. Our measurement of that distance with
a particular set of units is simply one way (of many) quantifying that fact.

With that measurement, we are simply comparing that something with another
thing that has been calibrated on a particular unit system. A ruler is the easiest
example of this: the ruler was constructed by marking off lines on a stick,
with each line determined from some other calibrated tool that the ruler manufacturer used.
Ultimately, all of these measurement tools refer back to some ultimate standard
for the unit (e.g. for a meter, it is length of the path travelled by light in a
vacuum in 1/(299,792,458) of a second.) As long as we all agree on this standard,
we can communicate with each other.
=#
#-
#=
However, we don't quite all agree. Most use SI units, but many still use imperial.
So we need conversion factors, and there is opportunity to make mistakes.
=#
#-
#=
### A more physically-meaningful reference quantity.
As we discussed earlier, when we measure something, we are making a comparison
with a standardized unit. We can think of this comparison as a ratio of the
physical quantity with the standard unit. For example, if we were to say the gap $H$ was
0.1 cm high, we could also write this statement as

$$\dfrac{H}{1\,\mathrm{cm}} = 0.1$$

However, this unit is a reference that we (humans) invented, so it is completely irrelevant
to physics.

On the other hand, if we instead choose as a reference some parameter
in the problem itself, then we are much more likely to have a *physically
meaningful* reference. In our example problem, the gap height is clearly physically
meaningful, so it serves as a natural *length scale*. The centerline velocity
is a natural *velocity scale*. These are our *reference scales*.

Actually, the statement that the dimensionless velocity plot has "no units" isn't
quite true. There *are* units on this plot. The units of position are "gap heights"
and the units of velocity are "centerline velocities". We can say, for example,
the: "the velocity at 0.2 gap heights from the lower wall is 0.64 times the centerline
velocity"). The huge advantage of this "local" system of units is that it requires no
standarization or international governing body. It is unambiguous to *anyone*,
including to an alien species.
=#

#-
#=
**In summary**, there are several advantages to the dimensionless view:
- There is no ambiguity about which system we used. It doesn't matter if
  the original parameters and values were in SI or imperial.
- Since the laws of physics hold regardless of our system of units, this sort of plot
    reveals true physical relationships between velocity and location.
- The dimensionless relationship between $u/U_c$ and $y/H$ is *universal*: it holds
  for *every* gap size and
  *every* flow speed. We know that this is true here because the exact velocity
  field was actually known to us:

  $$u(y) = \frac{4U_c}{H^2} y (H - y)$$

  But let's divide this by $U_c$ and group $y$ with $H$:

  $$\dfrac{u}{U_c} = 4 \dfrac{y}{H}\left(1 - \dfrac{y}{H}\right)$$

  This shows that $u/U_c$ depends on $y/H$ *and nothing
else*. There are no "free" parameters in this dimensionless relationship. We
will return to this point in the next notebook.
- We are lucky here because we already knew the mathematical equation. In most
cases, we aren't this lucky. But as we will discuss in the next notebook, once we
  ascertain the number of dimensionless quantities that fully describe a relationship,
  we can develop an *empirical correlation of data*. This correlation should *only* be
  pursued on dimensionless data. A correlation with dimensional coefficients is
  truly useless.
=#
#-
#=
### Dimensional homogeneity
One last point before move on. When we non-dimensionalize a quantity by a reference
quantity, they need to be in the same units. If they're not, then we need to
make sure the conversion factors are accounted for. For example, let's suppose
that the velocites in $u$ are measured in m/s, but the centerline velocity $U_c$ was
measured in ft/s. We obviously have to deal with the conversion between ft and
m to take the ratio $u/U_c$.

A conversion factor between two systems of units
is itself a dimensionless quantity. We can obtain the numerical factor with the
help of the `DimensionlessParameter` function. For example, let's consider the
conversion factor from feet to meters.
=#
DimensionlessParameter(1u"m"/1u"ft")

#=
We would multiply $U_c$ in ft/s by this factor to convert it to m/s.

But we can actually make use of a neat feature to skip this step, because
these software tools will carry out the conversion for us! For example, if
we divide $u = 0.5$ m/s by $U_c = 3$ ft/s:
=#
DimensionlessParameter(0.5u"m/s"/3u"ft/s")

#=
This name `DimensionlessParameter` isn't very descriptive! It might
be nice to use a more specific name for this quantity. Let's create
a new function that we can use to set this dimensionless velocity:
=#
@nondimvar DimensionlessVelocity
#=
Now we can use this to set the ratio of velocities:
=#
u_nd = DimensionlessVelocity(0.5u"m/s"/3u"ft/s")


text:
"""
A Fast and Accurate Dependency Parser using Neural Networks, Chen and Manning, EMNLP 2014
"""

type UnlabeledFeedForward
    word_f
    tag_f
    nonlinear
    W
end

type UExample
    wordids::Vector{Int}
    tagids::Union{Matrix{Float32},Vector{Int}}
    target::Int
    valid::Vector{Float32}
end

function UExample(s::State)
    w, t = unlabeledsparsefeatures(s)
    valid = [Float32(isvalid(s, act)) for act = 1:3]
    UExample(w, t, -1, valid)
end

targetsize(m::UnlabeledFeedForward) = size(m.W[end].w.data)[1]

function initmodel!(parser::DepParser, model::Type{UnlabeledFeedForward}, embed,
    topktags, nonlinear, sparsesizes, embedsizes, hiddensizes)
    T = Float32
    if embed == ""
        info("USING EMBEDDINGS WITH UNIFORM DISTRIBUTION [-0.01, 0.01]")
        word_f = Embedding(T, length(parser.words), embedsizes[1])
    else
        info("USING EMBEDDINGS LOADED FROM $(embed)")
        word_f = Embedding(embed, Float32)
        scale!(word_f)
    end
    if topktags
        tag_f = Linear(T, length(parser.tags), embedsizes[2])
    else
        tag_f = Embedding(T, length(parser.tags), embedsizes[2])
    end
    indim = sum(sparsesizes .* embedsizes)
    outdim = 3
    W = [myLinear(T, indim, hiddensizes[1]),
         myLinear(T, hiddensizes[1], outdim)]
    parser.model = UnlabeledFeedForward(word_f, tag_f, nonlinear, W)
    info("INPUT: [S^word,S^tag] = ", sparsesizes)
    info("EMBED DIMS: [word,tag] = ", embedsizes)
    info("HIDDEN LAYER: ", hiddensizes)
    info("OUTPUT DIM: ", outdim)
    info("NONLINEAR: ", nonlinear)
end

# TODO: make State have id field
# to tell where the State is in a batch
# called from expand(::State ::Int)
@compat function (m::UnlabeledFeedForward){T}(s::State{T}, act::Int)
    Var([0f0])
end

@compat function (m::UnlabeledFeedForward)(batch::AbstractVector{UExample}, istrain=true)
    wordvec, tagvec = [], []
    for s in batch
        push!(wordvec, s.wordids)
        push!(tagvec, s.tagids)
    end
    wordmat = m.word_f(Var(hcat(wordvec...)))
    tagmat = tagvectors(m.tag_f, tagvec)
    x = concat(1, wordmat, tagmat)
    x = m.nonlinear(m.W[1](x))
    x = dropout(x, 0.5, istrain)
    x = m.W[end](x)
    valids = hcat(map(x -> x.valid, batch)...)
    x.data .*= valids
    x
end

@compat function (m::UnlabeledFeedForward){T}(batch::AbstractVector{State{T}}, istrain=true)
    m(map(UExample, batch), istrain)
end

function unlabeledsparsefeatures(s::State)
    # word, tag
    b0 = tokenat(s, s.right)
    b1 = tokenat(s, s.right + 1)
    b2 = tokenat(s, s.right + 2)
    b3 = tokenat(s, s.right + 3)
    s0 = tokenat(s, s.top)
    s0l = tokenat(s, s.lchild)
    s0l2 = tokenat(s, s.lsibl.lchild)
    s0r = tokenat(s, s.rchild)
    s0r2 = tokenat(s, s.rsibl.rchild)
    s02l = tokenat(s, s.lchild.lchild)
    s12r = tokenat(s, s.rchild.rchild)
    s1 = tokenat(s, s.left)
    s1l = tokenat(s, s.left.lchild)
    s1l2 = tokenat(s, s.left.lsibl.lchild)
    s1r = tokenat(s, s.left.rchild)
    s1r2 = tokenat(s, s.left.rsibl.rchild)
    s12l = tokenat(s, s.left.lchild.lchild)
    s12r = tokenat(s, s.left.rchild.rchild)
    s2 = tokenat(s, s.left.left)
    s3 = tokenat(s, s.left.left.left)

    words = [b0.word, b1.word, b2.word, b3.word, s0.word, s0l.word, s0l2.word,
    s0r.word, s0r2.word, s02l.word, s12r.word, s1.word, s1l.word, s1l2.word,
    s1r.word, s1r2.word, s12l.word, s12r.word, s2.word, s3.word]

    tags = [b0.tag, b1.tag, b2.tag, b3.tag, s0.tag, s0l.tag, s0l2.tag,
    s0r.tag, s0r2.tag, s02l.tag, s12r.tag, s1.tag, s1l.tag, s1l2.tag,
    s1r.tag, s1r2.tag, s12l.tag, s12r.tag, s2.tag, s3.tag]

    return words, tags
end

typealias Doc Vector{Vector{Token}}

function train!{T}(::Type{UnlabeledFeedForward}, parser::DepParser{T}, trainsents::Doc,
    testsents::Doc=Vector{Token}[]; embed="", batchsize=32, iter=20, progbar=true,
    nonlinear=tanh, sparsesizes=[20,20] ,embedsizes=[50,50], hiddensizes=[1024],
    topktags=false, opt=MyAdaGrad(0.01), evaliter=100, outfile="parser.dat")
    info("WILL RUN $iter ITERATIONS")

    saver = ModelSaver(outfile)
    initmodel!(parser, UnlabeledFeedForward, embed, topktags,
               nonlinear, sparsesizes, embedsizes, hiddensizes)

    projectivize!(trainsents, false)

    trainsamples = UExample[]
    for s in trainsents
        s = State(s, parser)
        while !isfinal(s)
            ex = UExample(s)
            s = expandgold(s)[1]
            ex.target = s.prevact
            push!(trainsamples, ex)
        end
    end
    samplesize = length(trainsamples)
    labelsize = targetsize(parser.model)

    info("OPTIMIZER: ", typeof(opt))
    info("BATCH SIZE: $batchsize LABEL SIZE: $labelsize")
    info("#BATCHES: $(div(samplesize, batchsize))")
    info("#SAMPLES: $(samplesize)")

    for i = 1:iter
        batch = sub(trainsamples, rand(1:samplesize, batchsize))
        preds = parser.model(batch)
        golds = map(s -> s.target, batch)
        correct = reduce(0, zip(argmax(preds.data, 1), golds)) do v, tup
            v + Int(tup[1] == tup[2])
        end
        accuracy = correct / batchsize
        loss = update!(opt, golds, preds)
        println(now(), " Epoch: $i\tLoss: $loss\tAcc: $accuracy")

        if i % evaliter == 0 && !isempty(testsents)
            println()
            info("**ITER $i TESTING**")
            @time res = decode(UnlabeledFeedForward, parser, testsents)
            uas, las = evaluate(parser, res)
            saver(parser, uas)
            println()
        end
    end
end

function decode{T}(::Type{UnlabeledFeedForward}, parser::DepParser{T}, sents::Doc;
    batchsize=32, progbar=true)
    batches = 1:batchsize:length(sents)
    res = State{T}[]
    for k in batches
        batch = k:min(k+batchsize-1, endof(sents))
        ss = map(s -> State(s, parser), sub(sents, batch))
        parsegreedy!(parser, ss)
        append!(res, ss)
    end
    res
end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP