text:
get_codigo_elemento(::Posicao)::String = "84"

function get_valor_elemento(p::Posicao)
    if p.codigo == :onshore
        return "1"
    elseif p.codigo == :offshore
        return "2"
    else
        error("unicorn")
    end
end

CONST_POSICOES = [:onshore, :offshore]


text:
@testset "Split Sets" begin
    @testset "split_cols" begin
        cols = ["A", "B", "C", "D"]
        target = "A"
        exclude_target = exclude_elem(cols, target)
        @test cols == ["A", "B", "C", "D"]
        @test exclude_target == ["B", "C", "D"]

        target = "C"
        exclude_target2 = exclude_elem(cols, target)
        @test exclude_target2 == ["A", "B", "D"]
    end
end

@testset "Date Validation" begin
    @testset "Invalid date range" begin
        # one month
        _year = 2018
        _month = 1
        from_date = ZonedDateTime(_year, _month, 1, 0, tz"Asia/Seoul")
        to_date = ZonedDateTime(_year, _month, Dates.daysinmonth(DateTime(_year, _month)), 23, tz"Asia/Seoul")
        dates = collect(from_date:Dates.Hour(1):to_date)
        df = DataFrame(date = dates)

        @test_throws ArgumentError validate_dates(to_date, from_date, 1, df)
    end

    @testset "Invalid window size" begin
        # one month
        _year = 2018
        _month = 1
        from_date = ZonedDateTime(_year, _month, 1, 0, tz"Asia/Seoul")
        to_date = ZonedDateTime(_year, _month, 1, 23, tz"Asia/Seoul")
        dates = collect(from_date:Dates.Hour(1):to_date)
        df = DataFrame(date = dates)

        @test_throws BoundsError validate_dates(from_date, to_date, 48, df)
    end
end

@testset "Window" begin
    @testset "Basic window" begin
        df = DataFrame(
            A = 1:12,
            B = 13:24,
            C = 25:36)
        sample_size = 4
        output_size = 2
        dfs = window_df(df, sample_size, output_size)
        @test length(dfs) == 7
        @test dfs[1] == DataFrame(
            A = 1:6,
            B = 13:18,
            C = 25:30)
        @test dfs[7] == DataFrame(
            A = 7:12,
            B = 19:24,
            C = 31:36)
    end

    @testset "Window + Date Range" begin
        # one month
        _year = 2018
        _month = 1
        from_date = ZonedDateTime(_year, _month, 1, 0, tz"Asia/Seoul")
        to_date = ZonedDateTime(_year, _month, Dates.daysinmonth(DateTime(_year, _month)), 23, tz"Asia/Seoul")
        dates = collect(from_date:Dates.Hour(1):to_date)
        n = length(dates)

        df = DataFrame(
            A = 1:n,
            B = n+1:2n,
            C = 2n+1:3n,
            date = dates
        )
        sample_size = 12
        output_size = 2
        dfs = window_df(df, sample_size, output_size,
            ZonedDateTime(_year, _month, 10, 0, 1, tz"Asia/Seoul"),
            ZonedDateTime(_year, _month, 14, 23, 59, tz"Asia/Seoul"))

        @test length(dfs) == 106
        @test dfs[1] == DataFrame(
            A = 218:231,
            B = 962:975,
            C = 1706:1719,
            date = collect(ZonedDateTime(_year, _month, 10, 1, tz"Asia/Seoul"):Hour(1):ZonedDateTime(_year, _month, 10, 14, tz"Asia/Seoul"))
        )
        @test dfs[end] == DataFrame(
            A = 323:336,
            B = 1067:1080,
            C = 1811:1824,
            date = collect(ZonedDateTime(_year, _month, 14, 10, tz"Asia/Seoul"):Hour(1):ZonedDateTime(_year, _month, 14, 23, tz"Asia/Seoul"))
        )
    end

    @testset "Window + Date Range (Last day)" begin
        # one month
        _year = 2018
        _month = 1
        from_date = ZonedDateTime(_year, _month, 1, 0, tz"Asia/Seoul")
        to_date = ZonedDateTime(_year, _month, Dates.daysinmonth(DateTime(_year, _month)), 23, tz"Asia/Seoul")
        dates = collect(from_date:Dates.Hour(1):to_date)
        n = length(dates)

        df = DataFrame(
            A = 1:n,
            B = n+1:2n,
            C = 2n+1:3n,
            date = dates
        )

        sample_size = 12
        output_size = 2
        window_size = sample_size + output_size
        dfs = window_df(df, sample_size, output_size,
            ZonedDateTime(_year, _month, 30, 0, tz"Asia/Seoul"),
            ZonedDateTime(_year, _month, 31, 23, tz"Asia/Seoul"))

        @test length(dfs) == 35
        @test dfs[35] == DataFrame(
            A = 731:744,
            B = 1475:1488,
            C = 2219:2232,
            date = collect(ZonedDateTime(_year, _month, 31, 10, tz"Asia/Seoul"):Hour(1):ZonedDateTime(_year, _month, 31, 23, tz"Asia/Seoul"))
        )
    end

    @testset "Window + Date Range + StationCode" begin
        # one month
        _year = 2018
        _month = 1
        from_date = ZonedDateTime(_year, _month, 1, 0, tz"Asia/Seoul")
        to_date = ZonedDateTime(_year, _month, Dates.daysinmonth(DateTime(_year, _month)), 23, tz"Asia/Seoul")
        dates = collect(from_date:Dates.Hour(1):to_date)
        n = length(dates)
        baseCode = 1000
        tot_stn = 10

        listA = 1:n
        listB = n+1:2n
        listC = 2n+1:3n

        # stationCode 1 1 1 1 2 2 2 2 3 3 3 3 ...
        # A, B, C, date 1 2 3 4 1 2 3 4 1 2 3 4 ...
        df = DataFrame(
            stationCode = repeat(range(baseCode, step=1, length=tot_stn), inner=n),
            A = repeat(listA, outer=tot_stn),
            B = repeat(listB, outer=tot_stn),
            C = repeat(listC, outer=tot_stn),
            date = repeat(dates, outer=tot_stn))

        sample_size = 12
        output_size = 2
        window_size = sample_size + output_size

        target_stncode = baseCode + abs(rand(Int)) % tot_stn
        dfs = window_df(df, sample_size, output_size,
            ZonedDateTime(_year, _month, 30, 0, tz"Asia/Seoul"),
            ZonedDateTime(_year, _month, 31, 23, tz"Asia/Seoul"),
            target_stncode)

        df = dfs[1]
        @test DataFrames.nrow(df) == window_size
        @test df[!, :stationCode] == ones(window_size) * target_stncode
    end
end

@testset "Splitting Sizes" begin
    @testset "split_sizes3_mod0" begin
        total_size = 100
        batch_size = 10
        train_size, valid_size, test_size = split_sizes3(total_size, batch_size)
        @test train_size == 64
        @test valid_size == 16
        @test test_size == 20
    end

    @testset "split_sizes3_prime" begin
        total_size = 97
        batch_size = 10
        train_size, valid_size, test_size = split_sizes3(total_size, batch_size)
        @test train_size == 62
        @test valid_size == 16
        @test test_size == 19
    end

    @testset "split_sizes2_mod0" begin
        total_size = 100
        batch_size = 10
        train_size, valid_size = split_sizes2(total_size, batch_size)
        @test train_size == 80
        @test valid_size == 20
    end

    @testset "split_sizes2_prime" begin
        total_size = 97
        batch_size = 10
        train_size, valid_size = split_sizes2(total_size, batch_size)
        @test train_size == 78
        @test valid_size == 19
    end
end

@testset "Input" begin
    @testset "getX & getY" begin
        sample_size = 4
        output_size = 3
        window_size = sample_size + output_size

        start = 10
        lists = collect(Iterators.partition((start+1):(start + 3*window_size), window_size))

        df = DataFrame(
            A = lists[1],
            B = lists[2],
            C = lists[3])
        features = [:A, :B]

        X = getX(df, features, sample_size)
        @test X == Matrix(hcat([11, 12, 13, 14], [18, 19, 20, 21]))

        Y = getY(df, :C, sample_size)
        @test Y == Vector([29, 30, 31])
    end
end

@testset "DNN Input" begin
    @testset "DNN Pair construction" begin
        sample_size = 24
        output_size = 12

        _year = 2015
        _month = 1
        from_date = ZonedDateTime(_year, _month, 1, 0, tz"Asia/Seoul")
        to_date = ZonedDateTime(_year, _month, Dates.daysinmonth(DateTime(_year, _month)), 23, tz"Asia/Seoul")
        dates = from_date:Dates.Hour(1):to_date
        len_dates = length(dates)
        lists = collect(Iterators.partition(1:3*len_dates, len_dates))
        eltype = Integer

        df = DataFrame(
            date = dates,
            A = lists[1],
            B = lists[2],
            C = lists[3])

        features = [:A, :B]
        ycol = :C

        wd = window_df(df, sample_size, output_size, from_date, to_date)
        # pair = ([1,...,24,len_df,...,len_df+24],[2*len_df,...,2*len_df+12])
        pair = make_pair_DNN(wd[1], ycol, features, sample_size, output_size, eltype)

        # X length
        @test length(pair[1]) == sample_size * length(features)
        # Y length
        @test length(pair[2]) == output_size

        @test pair[1] == reduce(vcat, [collect(1:sample_size), collect(len_dates + 1:len_dates + sample_size)])
        @test pair[2] == collect((sample_size + 2*len_dates + 1):(sample_size + 2*len_dates + output_size))
    end

    @testset "DNN Batch(size 2) construction" begin
        sample_size = 24
        output_size = 12
        batch_size = 2

        _year = 2015
        _month = 1
        from_date = ZonedDateTime(_year, _month, 1, 0, tz"Asia/Seoul")
        to_date = ZonedDateTime(_year, _month, Dates.daysinmonth(DateTime(_year, _month)), 23, tz"Asia/Seoul")
        dates = from_date:Dates.Hour(1):to_date
        len_dates = length(dates)
        lists = collect(Iterators.partition(1:3*len_dates, len_dates))
        eltype = Integer

        df = DataFrame(
            date = dates,
            A = lists[1],
            B = lists[2],
            C = lists[3])

        features = [:A, :B]
        ycol = :C

        # type contraint to Floating point
        for col in [:A, :B, :C]
            df[!, col] = Float32.(df[!, col])
        end

        wd = window_df(df, sample_size, output_size, from_date, to_date)

        batch = make_batch_DNN(wd[1:batch_size], :C, [:A, :B], sample_size, output_size, batch_size,
            0.5, eltype)

        # indices for comparison
        in_chnks = [1:sample_size,
                    2:(sample_size+1)]
        out_chnks = [(sample_size+1):(sample_size+output_size),
                     (sample_size+2):(sample_size+output_size+1)]

        # Test X
        @test batch[1] == hcat(
            vcat(df[in_chnks[1], :A], df[in_chnks[1], :B]),
            vcat(df[in_chnks[2], :A], df[in_chnks[2], :B]))
        # Test Y
        @test batch[2] == hcat(
            vcat(df[out_chnks[1], :C]),
            vcat(df[out_chnks[2], :C]))

        @test size(batch[1]) == (sample_size * 2, batch_size)
        @test size(batch[2]) == (output_size, batch_size)
    end

    @testset "DNN Batch(size 3) construction" begin
        sample_size = 24
        output_size = 12
        batch_size = 3

        _year = 2015
        _month = 1
        from_date = ZonedDateTime(_year, _month, 1, 0, tz"Asia/Seoul")
        to_date = ZonedDateTime(_year, _month, Dates.daysinmonth(DateTime(_year, _month)), 23, tz"Asia/Seoul")
        dates = from_date:Dates.Hour(1):to_date
        len_dates = length(dates)
        lists = collect(Iterators.partition(1:3*len_dates, len_dates))
        eltype = Integer

        df = DataFrame(
            date = dates,
            A = lists[1],
            B = lists[2],
            C = lists[3])

        features = [:A, :B]
        ycol = :C

        # type contraint to Floating point
        for col in [:A, :B, :C]
            df[!, col] = Float32.(df[!, col])
        end

        wd = window_df(df, sample_size, output_size, from_date, to_date)

        batch = make_batch_DNN(wd[1:batch_size], :C, [:A, :B], sample_size, output_size, batch_size,
            0.5, eltype)

        # indices for comparison
        in_chnks = [1:sample_size,
                    2:(sample_size+1),
                    3:(sample_size+2)]
        out_chnks = [(sample_size+1):(sample_size+output_size),
                     (sample_size+2):(sample_size+output_size+1),
                     (sample_size+3):(sample_size+output_size+2)]
        # Test X
        @test batch[1] == hcat(
            vcat(df[in_chnks[1], :A], df[in_chnks[1], :B]),
            vcat(df[in_chnks[2], :A], df[in_chnks[2], :B]),
            vcat(df[in_chnks[3], :A], df[in_chnks[3], :B]))
        @test batch[2] == hcat(
            vcat(df[out_chnks[1], :C]),
            vcat(df[out_chnks[2], :C]),
            vcat(df[out_chnks[3], :C]))
        @test size(batch[1]) == (sample_size * 2, 3)
        @test size(batch[2]) == (output_size, 3)
    end
end

@testset "Validate Pairs" begin
    @testset "validate_pairs_normal" begin
        X = Matrix([ 1  2  3  4;
                     7  8  9 10;
                    13 14 15 16;
                    19 20 21 22;
                    25 26 27 28])
        Y = Matrix([5  6  7  8  9 10])

        remove_sparse_input!(X, Y, 0.5)

        @test X == Matrix([ 1  2  3  4;
                     7  8  9 10;
                    13 14 15 16;
                    19 20 21 22;
                    25 26 27 28])
        @test Y == Matrix([5  6  7  8  9 10])
    end

    @testset "validate_pairs_missings" begin
        X = Matrix([ 1  2  3  4;
                     7  8  9 10;
                    13 14 15 16;
                    19 20 21 22;
                    25 26 27 28])
        Y = Matrix([23 24 25 missing missing missing])

        remove_sparse_input!(X, Y, 0.5)

        @test X == Matrix([ 0.0 0.0 0.0 0.0;
                            0.0 0.0 0.0 0.0;
                            0.0 0.0 0.0 0.0;
                            0.0 0.0 0.0 0.0;
                            0.0 0.0 0.0 0.0;])
        @test Y == Matrix([0.0 0.0 0.0 0.0 0.0 0.0])
    end

    @testset "validate_pairs_missings_0.3" begin
        X = Matrix([ 1  2  3  4;
                     7  8  9 10;
                    13 14 15 16;
                    19 20 21 22;
                    25 26 27 28])
        Y = Matrix([23 24 25 26 missing missing])

        remove_sparse_input!(X, Y, 0.3)

        @test X == Matrix([ 0.0 0.0 0.0 0.0;
                            0.0 0.0 0.0 0.0;
                            0.0 0.0 0.0 0.0;
                            0.0 0.0 0.0 0.0;
                            0.0 0.0 0.0 0.0;])
        @test Y == Matrix([0.0 0.0 0.0 0.0 0.0 0.0])
    end

    @testset "validate_pairs_mx" begin
        X = Matrix([ 1  2  3  4;
                     7  8  9 10;
                    13 14 15 16;
                    19 20 21 22;
                    25 26 27 28])
        Y = Matrix([23 24 25 0.0 missing missing])

        remove_sparse_input!(X, Y, 0.5)

        @test X == Matrix([ 0.0 0.0 0.0 0.0;
                            0.0 0.0 0.0 0.0;
                            0.0 0.0 0.0 0.0;
                            0.0 0.0 0.0 0.0;
                            0.0 0.0 0.0 0.0;])
        @test Y == Matrix([0.0 0.0 0.0 0.0 0.0 0.0])
    end
end

@testset "CNN to RNN Batch Input Conversion" begin
    @testset "Batch Serialization" begin
        m = 3
        n = 5
        a = reshape(1:m*n, m, n)
        @test size(a) == (m, n)

        # 2D -> Arrays of 1D
        b = serializeBatch(a)
        @test size(b) == (n,)
        @test size(collect(b[1])) == (m,)
    end

    @testset "Batch Serialization" begin
        m1 = 3
        m2 = 4
        n = 5
        X = reshape(1:m1*n, m1, n)
        Y = reshape(101:(100+m2*n), m2, n)
        Z = serializeBatch(X, Y)

        @test size(Z) == (n,)
        @test size(collect(Z[1])) == (2,)
        # row size should be maintained
        @test size(Z[1][1]) == (m1,)
        @test size(Z[1][2]) == (m2,)
    end

    @testset "Unpack Sequence 3" begin
        x = reshape(vec(1:24), (2, 3, 4))

        @test unpack_seq(x) == reshape([
            [1 7 13 19; 3 9 15 21; 5 11 17 23],
            [2 8 14 20; 4 10 16 22; 6 12 18 24]], 1, 1, 2)
    end

    @testset "Unpack Sequence 4" begin
        x = reshape(vec(1:24), (2, 1, 3, 4))

        @test unpack_seq(x) == reshape([
            [1 7 13 19; 3 9 15 21; 5 11 17 23],
            [2 8 14 20; 4 10 16 22; 6 12 18 24]], 1, 1, 2, 1)
    end
end

ext: jl


text:
## Field access through these functions to reserve dot-getting for keys
@inline getaxes(x::VarAxes) = getaxes(typeof(x))
@inline getaxes(Ax::Type{<:Axes}) where {Axes<:VarAxes} = map(x->x(), (Ax.types...,))

getaxes(x) = ()

# Get AbstractAxis index
@inline Base.getindex(::AbstractAxis, idx::FlatIdx) = ComponentIndex(idx)
@inline Base.getindex(ax::AbstractAxis, ::Colon) = ComponentIndex(:, ax)
@inline Base.getindex(::AbstractAxis{IdxMap}, s::Symbol) where IdxMap =
    ComponentIndex(getproperty(IdxMap, s))

# Get ComponentArray index
Base.@propagate_inbounds Base.getindex(x::ComponentArray, idx::CartesianIndex) = getdata(x)[idx]
Base.@propagate_inbounds Base.getindex(x::ComponentArray, idx::ComponentArrays.FlatIdx...) = getdata(x)[idx...]
Base.@propagate_inbounds Base.getindex(x::ComponentArray, ::Colon) = getdata(x)[:]
@inline Base.getindex(x::ComponentArray, ::Colon...) = x
Base.@propagate_inbounds Base.getindex(x::ComponentArray, idx...) = getindex(x, toval.(idx)...)
@inline Base.getindex(x::ComponentArray, idx::Val...) = _getindex(x, idx...)

# Set ComponentArray index
Base.@propagate_inbounds Base.setindex!(x::ComponentArray, v, idx::FlatIdx...) = setindex!(getdata(x), v, idx...)
Base.@propagate_inbounds Base.setindex!(x::ComponentArray, v, ::Colon) = setindex!(getdata(x), v, :)
Base.@propagate_inbounds Base.setindex!(x::ComponentArray, v, idx...) = setindex!(x, v, toval.(idx)...)
@inline Base.setindex!(x::ComponentArray, v, idx::Val...) = _setindex!(x, v, idx...)

# Property access for CVectors goes through _get/_setindex
@inline Base.getproperty(x::ComponentVector, s::Symbol) = _getindex(x, Val(s))
@inline Base.getproperty(x::ComponentVector, s::Val) = _getindex(x, s)

@inline Base.setproperty!(x::ComponentVector, s::Symbol, v) = _setindex!(x, v, Val(s))
@inline Base.setproperty!(x::ComponentVector, s::Val, v) = _setindex!(x, v, s)

# Generated get and set index methods to do all of the heavy lifting in the type domain
@generated function _getindex(x::ComponentArray, idx...)
    ci = getindex.(getaxes(x), getval.(idx))
    inds = map(i -> i.idx, ci)
    axs = map(i -> i.ax, ci)
    axs = remove_nulls(axs...)
    # the index must be valid after computing `ci`
    :(Base.@_inline_meta; @inbounds ComponentArray(Base.maybeview(getdata(x), $inds...), $axs...))
end

@generated function _setindex!(x::ComponentArray, v, idx...)
    ci = getindex.(getaxes(x), getval.(idx))
    inds = map(i -> i.idx, ci)
    # the index must be valid after computing `ci`
    return :(Base.@_inline_meta; @inbounds setindex!(getdata(x), v, $inds...))
end

Base.@propagate_inbounds Base.view(x::ComponentArray, idx...) = getindex(x, idx...)


text:
@testset "light" begin
    M.set_interface_mode(M.LightInterface())
    @test M.get_interface_mode() isa M.LightInterface

    # matrix object (:other)
    X   = zeros(3, 4)
    mX  = matrix(X)
    mtX = matrix(X; transpose=true)

    @test mX === X
    @test mtX == permutedims(X)

    # :other but not matrix
    X = (1, 2, 3, 4)
    @test_throws ArgumentError matrix(X)

    # :table
    X = (x=[1,2,3], y=[1,2,3])
    @test M.vtrait(X) isa Val{:table}
    @test_throws M.InterfaceError matrix(X)
end

@testset "full" begin
    M.set_interface_mode(M.FullInterface())
    @test M.get_interface_mode() isa M.FullInterface

    M.matrix(::M.FullInterface, ::Val{:table}, X; kw...) =
        Tables.matrix(X; kw...)

    X = (x=[1,2,3], y=[1,2,3])
    mX = matrix(X)
    @test mX isa Matrix
    @test mX == hcat(X.x, X.y)
end


text:
################################################################################
#
# The following file implements three different methods of performing linear
# transformations on CKKS encrypted vectors.
#
# The different methods that are implemented are:
#   - The naive method.
#   - The square matrix diagonal method by Halevi and Shoup.
#   - The rectangular matrix hybird method by Juvekar, Vaikuntanathan and
#   Chandrakasan.
#
# The documentation for each function should be read before use.
#
################################################################################

module LinearTransformations

    using SEAL

    export naive_matrix_multiply, square_matrix_multiply, hybrid_matrix_multiply

    """
        naive_matrix_multiply(matrix::AbstractMatrix,
                              vector::AbstractArray,
                              evaluator::Evaluator,
                              encoder::CKKSEncoder,
                              relinearization_keys::RelinKeys)

    Compute the matrix multiplication between the provided plaintext matrix and
    array of ciphertext values.

    In the naive approach each entry of the ciphertext vector is encrypted
    independently and element-wise multiplication between the matrix and vector
    entries are performed as usual. When adding the element-wise
    multiplications of a matrix row the addition is performed in a balanced
    (log₂(n)) fashion in order to account for the CKKS scaling differences that
    occur when performing multiplication.
    """
    function naive_matrix_multiply(matrix::AbstractMatrix,
                                   vector::AbstractArray,
                                   evaluator::Evaluator,
                                   encoder::CKKSEncoder,
                                   relinearization_keys::RelinKeys)
        rows = size(matrix, 1)
        cols = size(matrix, 2)

        @debug "Performing naive matrix multiplication" rows cols

        temp = Array{Ciphertext}(undef, rows)

        for row in 1:rows
            list = []

            @debug "Processing row $(row) / $(rows) of matrix"

            # Force the garbage collector to garbage collect the old list.
            GC.gc()

            for col in 1:cols
                m = matrix[row, col]
                c = vector[col]

                # Encode the matrix value with the same initial scale as the
                # ciphertext values scale and mod switch to the same level.
                m_plaintext = Plaintext()
                encode!(m_plaintext, Float64(m), scale(c), encoder)
                mod_switch_to_inplace!(m_plaintext, parms_id(c), evaluator)

                # Since the plaintext and ciphertext have the same scale and
                # are on the same level, we are now able to multiply them.
                r = Ciphertext()
                multiply_plain!(r, c, m_plaintext, evaluator)
                relinearize_inplace!(r, relinearization_keys, evaluator)
                rescale_to_next_inplace!(r, evaluator)

                push!(list, r)
            end

            temp[row] = sum_ciphertexts(list, evaluator)
        end

        temp
    end

    """
        square_matrix_multiply(matrix::AbstractMatrix,
                               vector::Ciphertext,
                               evaluator::Evaluator,
                               encoder::Encoder,
                               galois_keys::GaloisKeys,
                               relinearization_keys::RelinKeys)

    Compute the matrix multiplication between the provided square plaintext
    matrix and ciphertext vector.

    The diagonal matrix multiplication approach uses the method presented by
    Halevi and Shoup in "Algorithms in HElib" and does only work on square
    matrices where both the row and column count are a divisor of the slot
    count (power of 2) in order to handle rotations properly. As usual, the
    addition operations of the method are performed in a balanced (log₂(n))
    fashion.
    """
    function square_matrix_multiply(matrix::AbstractMatrix,
                                    vector::Ciphertext,
                                    evaluator::Evaluator,
                                    encoder::CKKSEncoder,
                                    galois_keys::GaloisKeys,
                                    relinearization_keys::RelinKeys)
        rows = size(matrix, 1)
        cols = size(matrix, 2)

        if rows != cols || ((rows != 0) && ((rows & (rows - 1)) != 0))
            throw(ArgumentError("Expected power 2 sized, square matrix, got: ($rows, $cols)."))
        end

        ciphertexts = diag_matrix_multiply(matrix,
                                           vector,
                                           evaluator,
                                           encoder,
                                           galois_keys,
                                           relinearization_keys)

        sum_ciphertexts(ciphertexts, evaluator)
    end

    """
        hybrid_matrix_multiply(matrix::AbstractMatrix,
                               vector::Ciphertext,
                               evaluator::Evaluator,
                               encoder::Encoder,
                               galois_keys::GaloisKeys,
                               relinearization_keys::RelinKeys)

    Compute the matrix multiplication between the provided rectangular
    plaintext matrix and ciphertext vector.

    The hybrid matrix multiplication approach uses the method presented by
    Juvekar, Vaikuntanathan and Chandrakasan in "GAZELLE: A Low Latency
    Framework for Secure Neural Network Inference" and does only work on
    rectangular matrices where the number of rows are fewer than the number of
    columns and both the row and column count are a divisor of the slot count
    (power of 2) in order to handle rotations properly. As usual, the addition
    operations of the method are performed in a balanced (log₂(n)) fashion.
    """
    function hybrid_matrix_multiply(matrix::AbstractMatrix,
                                    vector::Ciphertext,
                                    evaluator::Evaluator,
                                    encoder::CKKSEncoder,
                                    galois_keys::GaloisKeys,
                                    relinearization_keys::RelinKeys)
        rows = size(matrix, 1)
        cols = size(matrix, 2)

        if cols < rows || (rows != 0) && ((rows & (rows - 1)) != 0)
            throw(ArgumentError("Expected power 2 sized matrix with #rows <= #cols, got: ($rows, $cols)."))
        end

        ciphertexts = diag_matrix_multiply(matrix,
                                           vector,
                                           evaluator,
                                           encoder,
                                           galois_keys,
                                           relinearization_keys)

        temp = sum_ciphertexts(ciphertexts, evaluator)

        # The number of rotations is determined by the rectangular factor.
        rotations = cols / rows

        ciphertexts = []

        # Extract the vector rotations:
        for idx in 1:rotations
            r = Ciphertext()
            rotate_vector!(r, temp, ((idx - 1) * rows), galois_keys, evaluator)

            push!(ciphertexts, r)
        end

        sum_ciphertexts(ciphertexts, evaluator)
    end

    """
        diag_matrix_multiply(matrix::AbstractMatrix,
                             vector::Ciphertext,
                             evaluator::Evaluator,
                             encoder::CKKSEncoder,
                             galois_keys::GaloisKeys,
                             relinearization_keys::RelinKeys)

    Computes the element-wise multiplication between the (wrapped) diagonals of
    a matrix and a vector.
    """
    function diag_matrix_multiply(matrix::AbstractMatrix,
                                  vector::Ciphertext,
                                  evaluator::Evaluator,
                                  encoder::CKKSEncoder,
                                  galois_keys::GaloisKeys,
                                  relinearization_keys::RelinKeys)
        rows = size(matrix, 1)
        cols = size(matrix, 2)

        diagonals = []

        # Extract the diagonals.
        for offset in 0:(rows - 1)
            push!(diagonals, wrapped_diag(matrix, offset))
        end

        ciphertexts = []

        for (idx, diagonal) in enumerate(diagonals)
            d = Plaintext()
            encode!(d, diagonal, scale(vector), encoder)
            mod_switch_to_inplace!(d, parms_id(vector), evaluator)

            r = Ciphertext()
            rotate_vector!(r, vector, idx - 1, galois_keys, evaluator)

            temp = Ciphertext()
            multiply_plain!(temp, r, d, evaluator)
            relinearize_inplace!(temp, relinearization_keys, evaluator)
            rescale_to_next_inplace!(temp, evaluator)

            push!(ciphertexts, temp)
        end

        ciphertexts
    end

    """
        sum_ciphertexts(list:Array)

    Takes an even list of ciphertexts and adds them together.
    """
    function sum_ciphertexts(list::Array, evaluator::Evaluator)
        len = length(list)

        if len == 1
            return list[1]
        end

        result = Ciphertext()

        if len == 2
            add!(result, list[1], list[2], evaluator)
        else
            split = len ÷ 2

            sum1 = sum_ciphertexts(list[1:split], evaluator)
            sum2 = sum_ciphertexts(list[(split + 1):len], evaluator)

            add!(result, sum1, sum2, evaluator)
        end

        result
    end

    """
        wrapped_diag(matrix::AbstractMatrix, offset::Integer)

    Extracts the (wrapped) diagonal of the provided matrix.

    This function is intended to be used by the diagonal and hybrid matrix
    multiplication algorithms and is designed with that in mind (i.e. comes
    without warranty for any other case).
    """
    function wrapped_diag(matrix::AbstractMatrix, offset::Integer=0)
        rows = size(matrix, 1)
        cols = size(matrix, 2)

        if offset < 0 || rows <= offset
            throw(ArgumentError("Offset $(offset) out of range."))
        end

        row_idx = 0
        col_idx = offset

        indices = []

        for _ in 1:cols
            push!(indices, (((row_idx % rows) + 1), (col_idx % cols) + 1))

            row_idx += 1
            col_idx += 1
        end

        [matrix[row, col] for (row, col) in indices]
    end
end


text:
metadata_packages = [
    "IJulia"]
Pkg.init()
Pkg.update()

for package=metadata_packages
    Pkg.add(package)
end

Pkg.resolve()

max_stars_repo_name: jayme-anchante/learning-julia

max_issues_repo_name: jayme-anchante/learning-julia

max_forks_repo_name: jayme-anchante/learning-julia


text:
#Convert xyplots to radial, only works with plots having values for xAxis.data
function radial!(ec::EChart)
    ec.polar = ECharts.Polar()
    ec.angleAxis = ECharts.AngleAxis()
    ec.radiusAxis = ECharts.RadiusAxis(_type = "category", data = ec.xAxis[1].data, z = 10)
    [x.coordinateSystem = "polar" for x in ec.series]
    ec.yAxis = nothing
    ec.xAxis = nothing
    return ec
end


text:
#
#--------------------------------------#
export derivMat
#--------------------------------------#
"""
 Compute derivative matrix for lagrange
 interpolants on points [x]
"""
function derivMat(x)

n = length(x)

a = ones(1,n)
for i=1:n
    for j=1:(i-1) a[i]=a[i]*(x[i]-x[j]) end
    for j=(i+1):n a[i]=a[i]*(x[i]-x[j]) end
end
a = 1 ./ a # Barycentric weights

# diagonal elements
D = x .- x'
for i=1:n D[i,i] = 1. end
D = 1 ./ D
for i=1:n
    D[i,i] = 0.
    D[i,i] = sum(D[i,:])
end

# off-diagonal elements
for j=1:n for i=1:n
    if(i!=j) D[i,j] = a[j] / (a[i]*(x[i]-x[j])) end
end end

return D
end
#--------------------------------------#


text:
import LinearAlgebra: *, exp!, iszero, isone,
    isdiag, issymmetric, ishermitian,
    isposdef, isposdef!, istril, istriu,
    dot, norm, normalize,
    rank, cond, opnorm, norm, det, tr,
    transpose, adjoint, inv, pinv,
    lu, lu!, qr, qr!, schur, schur!,
    cholesky, cholesky!, hessenberg, hessenberg!, factorize,
    eigvals, eigvals!, svdvals, svdvals!, svd,
    sqrt, exp, log,
    sin, cos, tan, csc, sec, cot,
    asin, acos, atan, acsc, asec, acot,
    sinh, cosh, tanh, csch, sech, coth,
    asinh, acosh, atanh, acsch, asech, acoth

norm(a::Vector{Single32}) = reinterpret(Single32, norm(reinterpret(Float64, a)))

normalize(a::Vector{Single32}) = reinterpret(Single32, normalize(reinterpret(Float64, a)))

dot(a::Vector{Single32}, b::Vector{Single32}) =
    reinterpret(Single32, dot(reinterpret(Float64, a), reinterpret(Float64, b)))

(*)(a::Matrix{Single32}, b::Matrix{Single32}) =
    reinterpret(Single32, reinterpret(Float64,a) * reinterpret(Float64,b))

exp!(x::Array{Complex{Single32},2}) = Complex{Single32}.(LinearAlgebra.exp!(Complex{Float64}.(x)))

for Op in (:iszero, :isone, :isdiag, :issymmetric, :ishermitian,
           :isposdef, :isposdef!, :istril, :istriu)
    @eval $Op(x::Matrix{Single32}) = $Op(reinterpret(Float64,x))
end

for Op in (:rank, :cond, :opnorm)
    @eval $Op(x::Array{Single32,2}) =
        reinterpret(Single32, $Op(reinterpret(Float64,x)))
end

for Op in (:norm, :det, :tr)
    @eval $Op(x::Matrix{Single32}) =
        reinterpret(Single32, $Op(reinterpret(Float64,x)))
end

for Op in (:transpose, :adjoint, :inv, :pinv)
    @eval $Op(x::Matrix{Single32}) =
        reinterpret(Single32, ($Op(reinterpret(Float64,x))))
end

function LinearAlgebra.lu(x::Matrix{Single32})
    result = lu(reinterpret(Float64,x))
    return LinearAlgebra.LU{Single32,Array{Single32,2}}(result)
end

#=
function LinearAlgebra.qr(x::Matrix{Single32})
    result = qr(reinterpret(Float64,x))
    if typeof(result) <: LinearAlgebra.QRCompactWY
        LinearAlgebra.QRCompactWY{Single32,Array{Single32,2}}(result)
    else
        LinearAlgebra.QR{Single32,Array{Single32,2}}(result)
    end
end
=#

#=
for Op in (:lu, :lu!, :qr, :qr!, :schur, :schur!)
    @eval $Op(x::Array{Single32,2}) = Single32.($Op(Float64.(x)))
end

for Op in (:cholesky, :cholesky!, :hessenberg, :hessenberg!, :factorize)
    @eval $Op(x::Array{Single32,2}) = Single32.($Op(Float64.(x)))
end
=#

for Op in (:svdvals, :svdvals!)
    @eval $Op(x::Array{Single32,2}) = reinterpret(Single32,($Op(reinterpret(Float64,x))))
end

for Op in (:eigvals, :eigvals!)
    @eval $Op(x::Array{Single32,2}) = (Complex{Single32}).($Op(reinterpret(Float64,x)))
end

for Op in (:sqrt, :exp, :log,
           :sin, :cos, :tan, :csc, :sec, :cot,
           :asin, :acos, :atan, :acsc, :asec, :acot,
           :sinh, :cosh, :tanh, :csch, :sech, :coth,
           :asinh, :acosh, :atanh, :acsch, :asech, :acoth)
    @eval $Op(x::Array{Single32,2}) = Single32.($Op(Float64.(x)))
end


text:
# Autogenerated wrapper script for blis_jll for x86_64-w64-mingw32
export blis

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("blis")
JLLWrappers.@declare_library_product(blis, "libblis.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        blis,
        "bin\\libblis.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using DataFrames, CSV, Dates, Statistics
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing4.09"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
data = dropmissing(CSV.read(datadir("purchaseData.csv"), copycols=true))

# ╔═╡ d4830cb6-2442-11eb-2475-75e064e51c14
begin
	data[!,:Date] = Date.(data[!,:Date], "d/m/y")
	first(sort(data, :Date), 3)
end

# ╔═╡ d4834f32-2442-11eb-0787-f14cebd337da
first(filter(row -> row[:Price] > 50000, data),3 )

# ╔═╡ d483e94c-2442-11eb-1307-3f405b818ba6
begin
	categorical!(data, :Grade)
	first(data, 3)
end

# ╔═╡ d48e90ca-2442-11eb-0d28-95f2cec59a6a
by(data, :Grade, :Price =>
       x -> ( NumSold=length(x), AvgPrice=mean(x)) )

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing4.09"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═d4830cb6-2442-11eb-2475-75e064e51c14
# ╠═d4834f32-2442-11eb-0787-f14cebd337da
# ╠═d483e94c-2442-11eb-1307-3f405b818ba6
# ╠═d48e90ca-2442-11eb-0d28-95f2cec59a6a
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
include("load_kplogic.jl")

include("test_kpbase.jl")

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
x = range(0, stop=1, length=100)
plot(x, x.^2, x.^3, x.^4)


text:
function read_bin(name, ndims = 2)
    """
        reading from a binary file in 3 dimensions
    """

    file = open(name)
    n = Int(read(file, UInt16, 1)[1])
    p = Int(read(file, UInt16, 1)[1])
    q = 1
    if ndims == 3
        q = Int(read(file, UInt16, 1)[1])
    end
    M = read(file, UInt8, n*p*q)
    M = reshape(M, (n, p, q))
    close(file)
    return M
end


text:
using comorbid
using Test

@testset "comorbid.jl" begin
    @test isapprox(exact_burden_term([.1], [.8], [true]), [.08])
end


text:
module FileIO_gpu

import FileIO # save, load
using JLD2: JLD2, JLDWriteSession, jldopen, isgroup, lookup_offset
using Knet.KnetArrays: KnetPtr, KnetArray, Cptr
using Knet.Ops20: RNN
using AutoGrad: Param
using CUDA: CUDA, functional, CuPtr

include("serialize.jl"); export cpucopy, gpucopy
include("jld.jl"); export save, load, @save, @load

end


text:
# The constants below are strings instead of enums because there is a requirement that users
# should be able to define their own without changing PowerSimulations.

# Constraints
const ACTIVE = "active"
const ACTIVE_RANGE = "activerange"
const ACTIVE_RANGE_LB = "activerange_lb"
const ACTIVE_RANGE_UB = "activerange_ub"
const COMMITMENT = "commitment"
const DURATION = "duration"
const DURATION_DOWN = "duration_dn"
const DURATION_UP = "duration_up"
const ENERGY_CAPACITY = "energy_capacity"
const ENERGY_CAPACITY_UP = "energy_capacity_up"
const ENERGY_CAPACITY_DOWN = "energy_capacity_down"
const ENERGY_LIMIT = "energy_limit"
const ENERGY_TARGET = "energy_target"
const FEEDFORWARD = "FF"
const FEEDFORWARD_UB = "FF_ub"
const FEEDFORWARD_BIN = "FF_bin"
const FEEDFORWARD_INTEGRAL_LIMIT = "FF_integral"
const FLOW_LIMIT = "FlowLimit"
const FLOW_LIMIT_FROM_TO = "FlowLimitFT"
const FLOW_LIMIT_TO_FROM = "FlowLimitTF"
const FLOW_REACTIVE_POWER_FROM_TO = "FqFT"
const FLOW_REACTIVE_POWER_TO_FROM = "FqTF"
const FLOW_ACTIVE_POWER_FROM_TO = "FpFT"
const FLOW_ACTIVE_POWER_TO_FROM = "FpTF"
const FLOW_ACTIVE_POWER = "Fp"
const FLOW_REACTIVE_POWER = "Fq"
const INPUT_POWER_RANGE = "inputpower_range"
const OUTPUT_POWER_RANGE = "outputpower_range"
const RAMP = "ramp"
const RAMP_DOWN = "ramp_dn"
const RAMP_UP = "ramp_up"
const RATE_LIMIT = "RateLimit"
const RATE_LIMIT_FT = "RateLimitFT"
const RATE_LIMIT_TF = "RateLimitTF"
const REACTIVE = "reactive"
const REACTIVE_RANGE = "reactiverange"
const REQUIREMENT = "requirement"
const RESERVE_POWER = "reserve_power"
const RESERVE_ENERGY = "reserve_energy"
const INFLOW_RANGE = "inflowrange"
const ACTIVE_RANGE_IC = "active_range_ic"
const START_TYPE = "start_type"
const STARTUP_TIMELIMIT = "startup_timelimit"
const STARTUP_TIMELIMIT_WARM = "startup_timelimit_warm"
const STARTUP_TIMELIMIT_HOT = "startup_timelimit_warm"
const STARTUP_INITIAL_CONDITION = "startup_initial_condition"
const STARTUP_INITIAL_CONDITION_UB = "startup_initial_condition_ub"
const STARTUP_INITIAL_CONDITION_LB = "startup_initial_condition_lb"
const MUST_RUN = "must_run"
const MUST_RUN_LB = "must_run_lb"
const NODAL_BALANCE_ACTIVE = "nodal_balance_active"
const NODAL_BALANCE_REACTIVE = "nodal_balance_reactive"
const NETWORK_FLOW = "network_flow"

abstract type ConstraintType end

abstract type RangeConstraint <: ConstraintType end
abstract type EnergyBalanceConstraint <: ConstraintType end

function make_constraint_name(
    ::Type{T},
    ::Type{U},
    ::Type{V},
) where {T <: ConstraintType, U <: VariableType, V <: PSY.Device}
    return encode_symbol(T, make_variable_name(U, V))
end

make_constraint_name(cons_type, device_type) = encode_symbol(device_type, cons_type)
make_constraint_name(cons_type) = encode_symbol(cons_type)


text:
module PyQuickPOMDPs

using PyCall
inspect = pyimport("inspect")

import QuickPOMDPs

function QuickPOMDPs.preprocess(v::PyObject)
    if pybuiltin("callable")(v)
        return convert(Function, v)
    end

    return v
end

function QuickPOMDPs.preprocess(name::Union{Val{:reward},Val{:observation}}, v::PyObject)
    if pybuiltin("callable")(v)
        n_required, n_optional = argcounts(v)
        return vararg_closure(name, Val(n_required), Val(n_optional), v)
    end

    return v
end

# here is where the bass drops. Probably easier to file an issue and tag @zsunberg than try to understand this
macro provide_vararg_closure(basename::Symbol, n_required::Int, n_optional::Union{Int,Float64})
    name = Symbol(String(basename)*"_pyfunc_closure")
    if n_optional < 0
        argnames = [Symbol("arg"*string(i)) for i in 1:n_required]
        methods = [:($name($(argnames...), a...) = v($(argnames...), a...))]
    else
        methods = [:($name(a::Vararg{Any,$(n_required)}) = v(a...))]
        for i in 1:n_optional
            push!(methods, :($name(a::Vararg{Any,$(n_required+i)}) = v(a...)))
        end
    end
    fname = esc(:vararg_closure)
    return quote
        function $fname(::Val{$(Meta.quot(basename))}, ::Val{$n_required}, ::Val{$n_optional}, v)
            $(methods...)
            return $name
        end
    end
end

@provide_vararg_closure(observation, 1, 0)
@provide_vararg_closure(observation, 1, 1)
@provide_vararg_closure(observation, 1, 2)
@provide_vararg_closure(observation, 1, -1)
@provide_vararg_closure(observation, 2, 0)
@provide_vararg_closure(observation, 2, 1)
@provide_vararg_closure(observation, 2, -1)
@provide_vararg_closure(observation, 3, 0)
@provide_vararg_closure(observation, 3, -1)

@provide_vararg_closure(reward, 2, 0)
@provide_vararg_closure(reward, 2, 1)
@provide_vararg_closure(reward, 2, 2)
@provide_vararg_closure(reward, 2, -1)
@provide_vararg_closure(reward, 3, 0)
@provide_vararg_closure(reward, 3, 1)
@provide_vararg_closure(reward, 3, -1)
@provide_vararg_closure(reward, 4, 0)
@provide_vararg_closure(reward, 4, -1)

"""
Assuming there are no keyword arguments, return the number of required and optional arguments (-1 for arbitrarily many optional).
"""
function argcounts(v::PyObject)
    params = inspect.signature(v).parameters.values()
    n_required = count(params) do p
        if p.kind==p.POSITIONAL_ONLY || p.kind==p.POSITIONAL_OR_KEYWORD
            return p.default == p.empty
        else
            return false
        end
    end
    if any(p->p.kind==p.VAR_POSITIONAL, params)
        n_optional = -1
    else
        n_optional = count(params) do p
            if p.kind==p.POSITIONAL_ONLY || p.kind==p.POSITIONAL_OR_KEYWORD
                return p.default != p.empty
            else
                return false
            end
        end
    end
    return n_required, n_optional
end

end #module

max_stars_repo_name: JuliaPOMDP/QuickPOMDPs

max_issues_repo_name: JuliaPOMDP/QuickPOMDPs

max_forks_repo_name: JuliaPOMDP/QuickPOMDPs


text:
using ITensors,
      LinearAlgebra,
      QuadGK

βc() = 0.5*log(sqrt(2.0)+1.0)

function ising_free_energy(β::Real,J::Real=1.0)
  k = β*J
  c = cosh(2.0*k)
  s = sinh(2.0*k)
  xmin = 0.0
  xmax = π
  integrand(x) = log(c^2+sqrt(s^4+1-2*s^2*cos(x)))
  integral,err = quadgk(integrand, xmin, xmax)::Tuple{Float64,Float64}
  return -(log(2.0)+integral/π)/(2.0*β)
end

function ising_magnetization(β::Real)
 β>βc && return (1.0-sinh(2.0*β)^(-4))^(1/8)
 return 0.0
end

function ising_mpo(sh::Tuple{Index,Index},sv::Tuple{Index,Index},
                   β::Real,J::Real=1.0;
                   sz::Bool=false)
  d = dim(sh[1])
  T = ITensor(sh[1],sh[2],sv[1],sv[2])
    for i = 1:d
      T[i,i,i,i] = 1.0
    end
    sz && (T[1,1,1,1] = -T[1,1,1,1])
    Q = [exp(β*J) exp(-β*J); exp(-β*J) exp(β*J)]
    D,U = eigen(Symmetric(Q))
    sqrtQ = U*Diagonal(sqrt.(D))*U'
    Xh1 = ITensor(vec(sqrtQ),sh[1],sh[1]')
    Xh2 = ITensor(vec(sqrtQ),sh[2],sh[2]')
    Xv1 = ITensor(vec(sqrtQ),sv[1],sv[1]')
    Xv2 = ITensor(vec(sqrtQ),sv[2],sv[2]')
    T = replacetags(T*Xh1*Xh2*Xv1*Xv2,"1","0")
  return T
end

function ising_partition(sh,sv,β)
  ny,nx = size(sh)
  T = Matrix{ITensor}(undef,ny,nx)
  for iy = 1:ny, ix = 1:nx
    ixp = per(ix+1,nx)
    iyp = per(iy+1,ny)
    T[iy,ix] = ising_mpo((sh[iy,ix],sh[iy,ixp]),(sv[iy,ix],sv[iyp,ix]),β)
  end
  return T
end


text:
# ---
# title: 1557. Minimum Number of Vertices to Reach All Nodes
# id: problem1557
# author: Indigo
# date: 2021-01-22
# difficulty: Medium
# categories: Graph
# link: <https://leetcode.com/problems/minimum-number-of-vertices-to-reach-all-nodes/description/>
# hidden: true
# ---
#
# Given a **  directed acyclic graph**, with `n` vertices numbered from `0` to
# `n-1`, and an array `edges` where `edges[i] = [fromi, toi]` represents a
# directed edge from node `fromi` to node `toi`.
#
# Find _the smallest set of vertices from which all nodes in the graph are
# reachable_. It's guaranteed that a unique solution exists.
#
# Notice that you can return the vertices in any order.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/07/07/untitled22.png)
#
#
#
#     Input: n = 6, edges = [[0,1],[0,2],[2,5],[3,4],[4,2]]
#     Output: [0,3]
#     Explanation: It's not possible to reach all the nodes from a single vertex. From 0 we can reach [0,1,2,5]. From 3 we can reach [3,4,2,5]. So we output [0,3].
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/07/07/untitled.png)
#
#
#
#     Input: n = 5, edges = [[0,1],[2,1],[3,1],[1,4],[2,4]]
#     Output: [0,2,3]
#     Explanation: Notice that vertices 0, 3 and 2 are not reachable from any other node, so we must include them. Also any of these vertices can reach nodes 1 and 4.
#
#
#
#
# **Constraints:**
#
#   * `2 <= n <= 10^5`
#   * `1 <= edges.length <= min(10^5, n * (n - 1) / 2)`
#   * `edges[i].length == 2`
#   * `0 <= fromi, toi < n`
#   * All pairs `(fromi, toi)` are distinct.
#
#
## @lc code=start
using LeetCode

function find_smallest_set_of_vertices(n::Int, edges::Vector{Vector{Int}})::Vector{Int}
    set = Set(edge[2] for edge in edges)
    res = Int[]
    for i in 0:n-1
        (i in set) || push!(res, i)
    end
    res
end
## @lc code=end


text:
# Very simple Julia backend which is just for testing the implementation
# and can be used as a reference implementation
import Base: pointer, similar, size, copy!, convert
using Base: RefValue
struct JLArray{T, N} <: GPUArray{T, N}
    data::Array{T, N}
    size::NTuple{N, Int}
end

function showarray(io::IO, A::JLArray, repr::Bool)
    print(io, "CPU: ")
    showarray(io, Array(A), repr)
end
"""
Thread group local memory
"""
struct LocalMem{N, T}
    x::NTuple{N, Vector{T}}
end

size(x::JLArray) = x.size
pointer(x::JLArray) = pointer(x.data)
to_device(state, x::JLArray) = x.data
to_device(state, x::Tuple) = to_device.(state, x)
to_device(state, x::RefValue{<: JLArray}) = RefValue(to_device(state, x[]))
to_device(state, x) = x
# creates a `local` vector for each thread group
to_device(state, x::LocalMemory{T}) where T = LocalMem(ntuple(i-> Vector{T}(x.size), blockdim_x(state)))

to_blocks(state, x) = x
# unpacks local memory for each block
to_blocks(state, x::LocalMem) = x.x[blockidx_x(state)]

function JLArray{T, N}(size::NTuple{N, Integer}) where {T, N}
    JLArray{T, N}(Array{T, N}(size), size)
end

similar(::Type{<: JLArray}, ::Type{T}, size::Base.Dims{N}) where {T, N} = JLArray{T, N}(size)

function unsafe_reinterpret(::Type{T}, A::JLArray{ET}, size::NTuple{N, Integer}) where {T, ET, N}
    JLArray{T, N}(reinterpret(T, A.data, size), size)
end

function copy!(
        dest::Array{T}, d_offset::Integer,
        source::JLArray{T}, s_offset::Integer, amount::Integer
    ) where T
    copy!(dest, d_offset, source.data, s_offset, amount)
end
function copy!(
        dest::JLArray{T}, d_offset::Integer,
        source::Array{T}, s_offset::Integer, amount::Integer
    ) where T
    copy!(dest.data, d_offset, source, s_offset, amount)
    dest
end
function copy!(
        dest::JLArray{T}, d_offset::Integer,
        source::JLArray{T}, s_offset::Integer, amount::Integer
    ) where T
    copy!(dest.data, d_offset, source.data, s_offset, amount)
    dest
end

mutable struct JLState{N}
    blockdim::NTuple{N, Int}
    griddim::NTuple{N, Int}

    blockidx::NTuple{N, Int}
    threadidx::NTuple{N, Int}
    localmem_counter::Int
    localmems::Vector{Vector{Array}}
end

function JLState(threads::NTuple{N}, blockdim::NTuple{N}) where N
    idx = ntuple(i-> 1, Val{N})
    blockcount = prod(blockdim)
    lmems = [Vector{Array}(0) for i in 1:blockcount]
    JLState{N}(threads, blockdim, idx, idx, 0, lmems)
end

function JLState(state::JLState{N}, threadidx::NTuple{N}) where N
    JLState{N}(
        state.blockdim,
        state.griddim,
        state.blockidx,
        threadidx,
        0,
        state.localmems
    )
end

function LocalMemory(state::JLState, ::Type{T}, ::Val{N}, ::Val{C}) where {T, N, C}
    state.localmem_counter += 1
    lmems = state.localmems[blockidx_x(state)]
    # first invokation in block
    if length(lmems) < state.localmem_counter
        lmem = zeros(T, N...)
        push!(lmems, lmem)
        return lmem
    else
        return lmems[state.localmem_counter]
    end
end

function AbstractDeviceArray(ptr::Array, shape::NTuple{N, Integer}) where N
    reshape(ptr, Int.(shape))
end
function AbstractDeviceArray(ptr::Array, shape::Vararg{Integer, N}) where N
    reshape(ptr, Int.(shape))
end

function _gpu_call(f, A::JLArray, args::Tuple, blocks_threads::Tuple{T, T}) where T <: NTuple{N, Integer} where N
    blocks, threads = blocks_threads
    idx = ntuple(i-> 1, length(blocks))
    blockdim = blocks
    state = JLState(threads, blockdim)
    device_args = to_device.(state, args)
    tasks = Array{Task}(threads...)
    for blockidx in CartesianRange(blockdim)
        state.blockidx = blockidx.I
        block_args = to_blocks.(state, device_args)
        for threadidx in CartesianRange(Int.(threads))
            thread_state = JLState(state, threadidx.I)
            tasks[threadidx] = @async f(thread_state, block_args...)
        end
        for t in tasks
            wait(t)
        end
    end
    return
end

# "intrinsics"
struct JLDevice end
device(x::JLArray) = JLDevice()
threads(dev::JLDevice) = 256
blocks(dev::JLDevice) = (256, 256, 256)

@inline function synchronize_threads(::JLState)
    #=
    All threads are getting started asynchronously,so a yield will
    yield to the next execution of the same function, which should call yield
    at the exact same point in the program, leading to a chain of yields  effectively syncing
    the tasks (threads).
    =#
    yield()
    return
end

for (i, sym) in enumerate((:x, :y, :z))
    for f in (:blockidx, :blockdim, :threadidx, :griddim)
        fname = Symbol(string(f, '_', sym))
        @eval $fname(state::JLState) = UInt32(state.$f[$i])
    end
end

blas_module(::JLArray) = Base.LinAlg.BLAS
blasbuffer(A::JLArray) = A.data

import Base: *, plan_ifft!, plan_fft!, plan_fft, plan_ifft, size, plan_bfft, plan_bfft!
# defining our own plan type is the easiest way to pass around the plans in Base interface without ambiguities
struct FFTPlan{T}
    p::T
end
function plan_fft(A::JLArray; kw_args...)
    FFTPlan(plan_fft(A.data; kw_args...))
end
function plan_fft!(A::JLArray; kw_args...)
    FFTPlan(plan_fft!(A.data; kw_args...))
end
function plan_bfft!(A::JLArray; kw_args...)
    FFTPlan(plan_bfft!(A.data; kw_args...))
end
function plan_bfft(A::JLArray; kw_args...)
    FFTPlan(plan_bfft(A.data; kw_args...))
end
function plan_ifft!(A::JLArray; kw_args...)
    FFTPlan(plan_ifft!(A.data; kw_args...))
end
function plan_ifft(A::JLArray; kw_args...)
    FFTPlan(plan_ifft(A.data; kw_args...))
end

function *(plan::FFTPlan, A::JLArray)
    x = plan.p * A.data
    JLArray(x)
end


text:
include((@__DIR__)*"/../preamble.jl")
include("reflecting_model/model_def.jl")

sims = CSV.read((@__DIR__)*"/reflecting_model/transient_distribution/point_mass/data/sim.csv",DataFrame)

rng = StableRNGs.StableRNG(26092021)

sims = simulate(model,fixed_time(0.1),(X=sims.X,φ=sims.φ),rng)
pwd()
sim_data = DataFrame(t=sims.t,φ=sims.φ,X=sims.X,n=sims.n)
CSV.write((@__DIR__)*"/reflecting_model/transient_distribution/point_mass/datat2/sims.csv",
    sim_data,
)

x_vals=range(0,10;length=n_err_evals)

function fast_cdf(x::AbstractVector,data::DataFrame)
    sort!(data,:X)

    x1 = data.X[data.φ.==1]
    x2 = data.X[data.φ.==2]

    c1 = 1
    c2 = 1

    F1 = zeros(length(x))
    F2 = zeros(length(x))

    for (c,xi) in enumerate(x)
        tmp1 = findfirst(y->y>xi,x1[c1:end])
        if !(tmp1===nothing)
            c1 += tmp1-1
        else
            c1 = length(x1)+1
        end
        F1[c] = c1-1

        tmp2 = findfirst(y->y>xi,x2[c2:end])
        if !(tmp2===nothing)
            c2 += tmp2-1
        else
            c2 = length(x2)+1
        end
        F2[c] = c2-1
    end
    return F1./length(data.X), F2./length(data.X)
end

F1, F2 = fast_cdf(x_vals,sim_data)

CSV.write((@__DIR__)*"/reflecting_model/transient_distribution/point_mass/datat2/sims_cdf_evaluated.csv",DataFrame(x=x_vals,phase_1=F1,phase_2=F2))

function bootstrap(sim_data,x_vals,n_boot)
    rng = StableRNGs.StableRNG(1412)
    n_sims = length(sim_data.X)
    for n in 1:n_boot
        @show n
        ind = rand(rng,1:n_sims,n_sims)
        boot_data = sim_data[ind,:]
        F1, F2 = fast_cdf(x_vals,boot_data)
        CSV.write(
            (@__DIR__)*"/reflecting_model/transient_distribution/point_mass/datat2/bootstrap/sample_"*string(n)*".csv",
            DataFrame(x=x_vals,phase_1=F1,phase_2=F2),
        )
    end
    return nothing
end

bootstrap(sim_data,x_vals,n_boot)


text:
"""
    Revise.track(Base)
    Revise.track(Core.Compiler)
    Revise.track(stdlib)

Track updates to the code in Julia's `base` directory, `base/compiler`, or one of its
standard libraries.
"""
function track(mod::Module)
    if mod == Base
        # Determine when the basesrccache was built
        mtcache = mtime(basesrccache)
        # Initialize expression-tracking for files, and
        # note any modified since Base was built
        files = String[]
        for (submod, filename) in Base._included_files
            push!(fileinfos, filename=>FileInfo(submod, basesrccache))
            push!(files, filename)
            if mtime(filename) > mtcache
                push!(revision_queue, filename)
            end
        end
        # Add the files to the watch list
        init_watching(files)
    elseif mod == Core.Compiler
        compilerdir = joinpath(juliadir, "base", "compiler")
        track_subdir_from_git(Core.Compiler, compilerdir)
    elseif nameof(mod) ∈ stdlib_names
        stdlibdir = joinpath(juliadir, "stdlib")
        libdir = joinpath(stdlibdir, String(nameof(mod)), "src")
        track_subdir_from_git(mod, normpath(libdir))
    else
        error("no Revise.track recipe for module ", mod)
    end
    nothing
end

# For tracking subdirectories of Julia itself (base/compiler, stdlibs)
function track_subdir_from_git(mod::Module, subdir::AbstractString)
    # diff against files at the same commit used to build Julia
    repo, repo_path = git_repo(subdir)
    if repo == nothing
        error("could not find git repository at $subdir")
    end
    prefix = subdir[length(repo_path)+2:end]   # git-relative path of this subdir
    tree = git_tree(repo, Base.GIT_VERSION_INFO.commit)
    files = Iterators.filter(file->startswith(file, prefix) && endswith(file, ".jl"), keys(tree))
    ccall((:giterr_clear, :libgit2), Cvoid, ())  # necessary to avoid errors like "the global/xdg file 'attributes' doesn't exist: No such file or directory"
    wfiles = String[]  # files to watch
    for file in files
        fullpath = joinpath(repo_path, file)
        src = git_source(file, tree)
        if src != read(fullpath, String)
            push!(revision_queue, fullpath)
        end
        fmod = get(juliaf2m, fullpath, Core.Compiler)  # Core.Compiler is not cached
        fi = FileInfo(fmod)
        if parse_source!(fi.fm, src, Symbol(file), 1, fmod) === nothing
            warn("failed to parse Git source text for ", file)
        else
            instantiate_sigs!(fi.fm)
        end
        fileinfos[fullpath] = fi
        push!(wfiles, fullpath)
    end
    init_watching(wfiles)
    return nothing
end

# For tracking Julia's own stdlibs
const stdlib_names = Set([
    :Base64, :CRC32c, :Dates, :DelimitedFiles, :Distributed,
    :FileWatching, :Future, :InteractiveUtils, :Libdl,
    :LibGit2, :LinearAlgebra, :Logging, :Markdown, :Mmap,
    :OldPkg, :Pkg, :Printf, :Profile, :Random, :REPL,
    :Serialization, :SHA, :SharedArrays, :Sockets, :SparseArrays,
    :Statistics, :SuiteSparse, :Test, :Unicode, :UUIDs])

# This replacement is needed because the path written during compilation differs from
# the git source path
const stdlib_rep = joinpath("usr", "share", "julia", "stdlib", "v$(VERSION.major).$(VERSION.minor)") => "stdlib"

const juliaf2m = Dict(normpath(replace(file, stdlib_rep))=>mod
    for (mod,file) in Base._included_files)


text:
function saveFEMasVTK( filename,
                       x,
                       eConn,
                       scalarLabels, scalars,
                       vectorLabels, vectors )
#  Saves a finite element solution using the WriteVTK package
#
#  Author: Jeff Borggaard, Virginia Tech
#          part of FEMfunctions.jl
#
#  Licensing:
#     This code is distributed under the MIT license.
#
#  requires the package WriteVTK

  #  Get mesh dimensions
  nElements, nNodesPerElement = size(eConn)
  nNodes, nDim = size(x)

  #  Place the element connectivity in the correct format
  if (nNodesPerElement==3)
    cellType = VTKCellTypes.VTK_TRIANGLE              # =  5 in the VTK 4.2 format
  elseif ( nNodesPerElement==4 && nDim==2)
    cellType = VTKCellTypes.VTK_QUAD                  # =  9 in the VTK 4.2 format
  elseif ( nNodesPerElement==4 && nDim==3)
    cellType = VTKCellTypes.VTK_TETRA                 # = 10 in the VTK 4.2 format
  elseif ( nNodesPerElement==6 )
    cellType = VTKCellTypes.VTK_QUADRATIC_TRIANGLE    # = 22 in the VTK 4.2 format
  elseif ( nNodesPerElement==8 )
    cellType = VTKCellTypes.VTK_HEXAHEDRON            # = 12 in the VTK 4.2 format
  elseif ( nNodesPerElement==10 )
    cellType = VTKCellTypes.VTK_QUADRATIC_TETRA       # = 24 in the VTK 4.2 format
  end

  cells = MeshCell[]
  for i=1:nElements
    c = MeshCell(cellType,eConn[i,:])
    cells = push!(cells, c)
  end

  #  Open the VTK file and write out the geometry data
  if ( nDim==2 )
    vtk = vtk_grid(filename, x[:,1], x[:,2], cells)
  elseif ( nDim==3 )
    vtk = vtk_grid(filename, x[:,1], x[:,2], x[:,3], cells)
  else
    # warning about mesh
  end

  #  Write out scalar data
  nScalars = length(scalarLabels)
  for i=1:nScalars
    vtk_point_data(vtk, scalars[:,i], scalarLabels[i])
  end

  #  Write out vector data
  nVectors = length(vectorLabels)
  if ( nDim==2 )
    for i=1:nVectors
      i1,i2 = 2*(i-1)+1, 2*(i-1)+2
      vec = (vectors[:,i1],vectors[:,i2])
      vtk_point_data(vtk, vec, vectorLabels[i])
    end
  else
    for i=1:nVectors
      i1,i2,i3 = 3*(i-1)+1, 3*(i-1)+2, 3*(i-1)+3
      vec = (vectors[:,i1],vectors[:,i2],vectors[:,i3])
      vtk_point_data(vtk, vec, vectorLabels[i])
    end
  end

  #  Close the file
  outfiles = String[]
  append!(outfiles, vtk_save(vtk))

end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

struct DensityWithDiff{D<:AbstractDensity,VJP<:DifferentiationAlgorithm} <: AbstractDensity
    vjpalg::VJP
    density::D
end

@inline Base.parent(density::DensityWithDiff) = density.density

vjp_algorithm(density::DensityWithDiff) = density.vjpalg

Base.:(|)(vjpalg::DifferentiationAlgorithm, density::AbstractDensity) = DensityWithDiff(vjpalg, density)
Base.:(|)(vjpalg::DifferentiationAlgorithm, density::DensityWithDiff) = DensityWithDiff(vjpalg, parent(density))
Base.:(|)(vjpalg::DifferentiationAlgorithm, density::PosteriorDensity) = PosteriorDensity(vjpalg | density.likelihood, vjpalg | density.prior)

function Base.show(io::IO, density::DensityWithDiff)
    vjpalg = vjp_algorithm(density)
    Base.show(io, vjpalg)
    print(io, "|")
    Base.show(io, parent(density))
end

var_bounds(density::DensityWithDiff) = var_bounds(density.density)

ValueShapes.varshape(density::DensityWithDiff) = varshape(density.density)

ValueShapes.unshaped(density::DensityWithDiff) = DensityWithDiff(density.vjpalg, unshaped(density.density))

@inline DensityInterface.logdensityof(density::DensityWithDiff, v::Any) = logdensityof(density.density, v)

@inline checked_logdensityof(density::DensityWithDiff, v::Any) = checked_logdensityof(density.density, v)

@inline get_initsrc_from_target(target::DensityWithDiff) = get_initsrc_from_target(target.density)

function bat_transform_impl(target::NoDensityTransform, density::DensityWithDiff, algorithm::DensityIdentityTransform)
    (result = density, trafo = identity)
end

function bat_transform_impl(target::AbstractDensityTransformTarget, density::DensityWithDiff, algorithm::TransformAlgorithm)
    transformed_density, trafo = bat_transform_impl(target, parent(density), algorithm)
    (result = DensityWithDiff(vjp_algorithm(density), transformed_density), trafo = trafo)
end


text:
using ArchGDAL
using DelimitedFiles
using DataFrames
import GDAL
# You need to load rastersweep.jl before running this.

"""
Gets the canonical names and ids of each admin unit.
```julia
admins = unique(df[!, ["admin1_name", "admin1_id"]])
```
"""
function uganda_canonical()
   csv_path = expanduser(joinpath(
      "~", "data", "projects", "uganda2020", "outputs", "uga_canonical",
      "uga_canonical_names.csv"))
   csv, header = readdlm(csv_path, ',', '\n'; header = true)
   df = DataFrame(csv, reshape(header, length(header)))
   for scol in ["admin1_name", "admin2_name", "admin3_name"]
      df[scol] = String.(df[scol])
   end
   for icol in ["admin1_id", "admin2_id", "admin3_id"]
      df[icol] = Int.(df[icol])
   end
   df
end

"""
Given a single feature, this finds all overlapping pixels
in a grid. It then marks those pixels as overlapping this feature.
"""
function build_identify_district(rasterlist)
   function identify_district(geometry, featuredict)::UInt16
      admin_cnt = length(rasterlist)
      admin = [
         UInt16(featuredict["admin1id"]),
         UInt16(featuredict["admin2id"]),
         UInt16(featuredict["admin3id"])
         ]
      envelope = ArchGDAL.envelope(geometry)
      # upper-left on the left column. lower-right in the right column.
      rect = [envelope.MinX envelope.MaxX; envelope.MaxY envelope.MinY]
      raster = rasterlist[1]
      indices = cover_rect(rect, raster.geo, size(raster.data))
      total = 0.0
      # We use a GDAL point to test for inside polygon.
      geom_point = ArchGDAL.createpoint()
      ArchGDAL.addpoint!(geom_point, 0.0, 0.0)
      for idx in indices
         xy = pixel_center(raster.geo, [idx[1], idx[2]])
         ArchGDAL.setpoint!(geom_point, 0, xy[1], xy[2])
         if ArchGDAL.within(geom_point, geometry)
            for aidx in 1:admin_cnt
               rasterlist[aidx].data[idx] = UInt16(admin[aidx])
            end
         end
      end
      admin[3]
   end
end

"""
    layerinfo(layer)

Prints a bunch of stuff about an ArchGDAL layer, if you're curious.
"""
function layerinfo(layer)
   @show ArchGDAL.getname(layer)
   @show ArchGDAL.nfeature(layer)
   featuredefn = ArchGDAL.layerdefn(layer)
   nfield = ArchGDAL.nfield(featuredefn)
   ngeom = ArchGDAL.nfield(featuredefn)
   @show ngeom
   for fieldidx in 1:nfield
      field = ArchGDAL.getfielddefn(featuredefn, fieldidx - 1)
      name = ArchGDAL.getname(field)
      type = ArchGDAL.gettype(field)
      @show (name, type, fieldidx)
   end
   for geomidx in 1:1
      geom = ArchGDAL.getgeomdefn(featuredefn, geomidx - 1)
      name = ArchGDAL.getname(geom)
      type = ArchGDAL.gettype(geom)
      @show (name, type, geomidx)
   end
end

"""
    featureinfo(featuredefn, feature)

Makes a dictionary with all fields for this feature.
If you think of a shapefile as a dataframe, this returns a row
of the dataframe.

```julia
featuredefn = ArchGDAL.layerdefn(layer)
feature = ArchGDAL.getfeature(layer, 7)
featuredict = featurinfo(featuredefn, feature)
```
"""
function featureinfo(featuredefn, feature)
   nfield = ArchGDAL.nfield(featuredefn)
   values = Dict{String,Any}()
   for fieldidx in 1:nfield
      field = ArchGDAL.getfielddefn(featuredefn, fieldidx - 1)
      name = ArchGDAL.getname(field)
      value = ArchGDAL.getfield(feature, fieldidx - 1)
      values[name] = value
   end
   values
end

"""
Loops through a vector file, feature-by-feature.
"""
function overgeometry(toapply, layer; limit = typemax(Int))
   featuredefn = ArchGDAL.layerdefn(layer)
   nfeature = min(ArchGDAL.nfeature(layer), limit)
   result = nothing
   for featureidx in 1:nfeature
      ArchGDAL.getfeature(layer, featureidx - 1) do feature
         infodict = featureinfo(featuredefn, feature)
         geom = ArchGDAL.getgeom(feature, 0)
         envelope = ArchGDAL.envelope(geom)
         status = toapply(geom, infodict)
         if result === nothing
            result = fill(("", status), nfeature)
         end
         admin = infodict["admin3name"]
         @show (featureidx, admin)
         result[featureidx] = (admin, status)
      end
   end
   result
end

asrect(envelope) = [envelope.MinX envelope.MaxX; envelope.MaxY envelope.MinY]
unionrect(a, b) = [
      min(a[1,1], b[1,1]) max(a[1,2], b[1,2]); max(a[2,1], b[2,1]) min(a[2,2], b[2,2])
      ]

"""
Loops through a vector file, iterating over bounding boxes.
You choose the admin area over which to loop.
"""
function overenvelope(toapply, layer; limit = typemax(Int))
   featuredefn = ArchGDAL.layerdefn(layer)
   nfeature = min(ArchGDAL.nfeature(layer), limit)
   # Find the type of the envelope.
   admin1_envelope = Dict{Int,Array{Float64,2}}()
   # Combine feature envelopes to get admin1 envelopes.
   for featureidx in 1:nfeature
      ArchGDAL.getfeature(layer, featureidx - 1) do feature
         # The feature is at admin3, but it's in an admin1.
         admin1 = featureinfo(featuredefn, feature)["admin1id"]
         envelope = asrect(ArchGDAL.envelope(ArchGDAL.getgeom(feature, 0)))
         if haskey(admin1_envelope, admin1)
            admin1_envelope[admin1] = unionrect(admin1_envelope[admin1], envelope)
         else
            admin1_envelope[admin1] = envelope
         end
      end
   end
   @show length(admin1_envelope)

   # Iterate over envelopes.
   result = nothing
   admin_idx = 0
   for (admin1id, envelope) in admin1_envelope
         status = toapply(envelope, admin1id)
         if result === nothing
            result = fill((0, status), length(admin1_envelope))
         end
         @show (admin1id, status)
         admin_idx += 1
         result[admin_idx] = (admin1id, status)
   end
   result
end

# You need to include rastersweep.jl
struct RasterRead
   data
   no_data
   geo::Geo
end

"""
Returns a functor that assigns a weight to each HRSL pixel.
That weight is the fraction of the district population that lives
in this pixel. This searches by asking whether the pixel is inside
the feature, using computational geometry from GDAL.
"""
function build_countup(raster)
   function count_up(geometry, featuredict)::Float64
      envelope = ArchGDAL.envelope(geometry)
      # upper-left on the left column. lower-right in the right column.
      rect = [envelope.MinX envelope.MaxX; envelope.MaxY envelope.MinY]
      indices = cover_rect(rect, raster.geo, size(raster.data))
      total = 0.0
      # We use a GDAL point to test for inside polygon.
      geom_point = ArchGDAL.createpoint()
      ArchGDAL.addpoint!(geom_point, 0.0, 0.0)
      within = fill(CartesianIndex(0, 0), length(indices))
      inside_cnt = 0
      for idx in indices
         if raster.data[idx] != raster.no_data
            xy = pixel_center(raster.geo, [idx[1], idx[2]])
            ArchGDAL.setpoint!(geom_point, 0, xy[1], xy[2])
            if ArchGDAL.within(geom_point, geometry)
               value = raster.data[idx]
               @assert value >= 0
               total += value
               inside_cnt += 1
               within[inside_cnt] = idx
            end
         end
      end
      largest = 0.0
      for inside_idx in 1:inside_cnt
         idx = within[inside_idx]
         # All values are 1000 times the actual fraction for numerical accuracy.
         pop_rate = 1000 * raster.data[idx] / total
         @assert pop_rate >= 0
         raster.data[idx] = pop_rate
         largest = max(largest, pop_rate)
      end
      largest
   end
end

"""
Returns a functor that assigns a weight to each HRSL pixel.
That weight is the fraction of the district population that lives
in this pixel. This knows a pixel is in a district by checking a large
lookup table, called the idraster.
"""
function build_fractional(popraster, idraster)
   function count_fractions(rect, admin_id)::Float64
      # rect is upper-left on the left column. lower-right in the right column.
      indices = cover_rect(rect, popraster.geo, size(popraster.data))
      total = 0.0
      geom_point = ArchGDAL.createpoint()
      ArchGDAL.addpoint!(geom_point, 0.0, 0.0)
      within = fill(CartesianIndex(0, 0), length(indices))
      inside_cnt = 0
      for idx in indices
         if popraster.data[idx] != popraster.no_data
            if idraster.data[idx] == admin_id
               value = popraster.data[idx]
               @assert value >= 0
               total += value
               inside_cnt += 1
               within[inside_cnt] = idx
            end
         end
      end
      largest = 0.0
      for inside_idx in 1:inside_cnt
         idx = within[inside_idx]
         # All values are 1000 times the actual fraction for numerical accuracy.
         pop_rate = 1000 * popraster.data[idx] / total
         @assert pop_rate >= 0
         popraster.data[idx] = pop_rate
         largest = max(largest, pop_rate)
      end
      largest
   end
end

"""
Make a dataset with the same size and transform as the given dataset
but it has two bands, one with floats and one with integers.
Use 0 as the nodata value.
"""
function create_fraction(dataset, path, admin_cnt = 3)
   b = ArchGDAL.getband(dataset, 1)
   w = ArchGDAL.width(b)
   h = ArchGDAL.height(b)
   driver = ArchGDAL.getdriver("GTiff")
   options = [
      "TILED=YES",
      "BLOCKXSIZE=256",
      "BLOCKYSIZE=256",
      "COMPRESS=LZW",
   ]
   ds = ArchGDAL.create(
         path;
         driver = driver,
         width = w,
         height = h,
         nbands = admin_cnt,
         dtype = UInt16,
         options = options,
         )
   ArchGDAL.setgeotransform!(ds, ArchGDAL.getgeotransform(dataset))
   ArchGDAL.setproj!(ds, ArchGDAL.getproj(dataset))
   for bandidx in 1:admin_cnt
      ArchGDAL.getband(ds, bandidx) do band
         ArchGDAL.setnodatavalue!(band, 0)
         ArchGDAL.fillraster!(band, 0)
         ArchGDAL.setcolorinterp!(band, GDAL.GCI_GrayIndex)
         access = ArchGDAL.accessflag(band)
         @assert access == GDAL.GA_Update
      end
   end
   ds
end

"""
Make a TIFF that tells you, for each pixel, which admin units it belongs to.
There is a band each for admin1, admin2, and admin3. The values are the index
of that admin in the canonical list of admins.

You can read them in R this way.
```R
admin2id = raster::raster("~/.julia/dev/simplexpop/hrsl_admin_id.tif", band = 2)
```
"""
function make_index_tiff(; limit = typemax(Int))
   hrsl_path = "hrsl_ls_adjusted.tif"
   admin3_path = expanduser(joinpath("~", "data", "projects", "uganda2020", "outputs",
      "uga_canonical", "uga_canonical_names.shp"))
   copy_path = abspath("hrsl_admin_id.tif")
   @show copy_path
   admin_cnt = 3
   if ispath(copy_path)
      println("$(copy_path) exists, deleting.")
      rm(copy_path)
   end
   id_raster_ds = ArchGDAL.read(hrsl_path) do hrsl_ds
      create_fraction(hrsl_ds, copy_path, admin_cnt)
   end
   geo = band_geo(ArchGDAL.getgeotransform(id_raster_ds))
   hrsl_no_data = 0
   hrsl_raster = [
      RasterRead(ArchGDAL.read(ArchGDAL.getband(id_raster_ds, i)), hrsl_no_data, geo)
      for i in 1:admin_cnt
      ]
   identify_district = build_identify_district(hrsl_raster)

   ArchGDAL.read(admin3_path) do sc_dataset
      @show ArchGDAL.nlayer(sc_dataset)
      admin = ArchGDAL.getlayer(sc_dataset, 0)
      layerinfo(admin)
      totals = overgeometry(identify_district, admin; limit = limit)
      for aidx in 1:admin_cnt
         ArchGDAL.write!(id_raster_ds, hrsl_raster[aidx].data, aidx)
      end
   end
   ArchGDAL.destroy(id_raster_ds)
   ispath(copy_path)
end

function assign_to_admin_unit(; limit = typemax(Int))
   hrsl_path = expanduser("~/data/projects/uganda2020/outputs/hrsl_ls_adjusted.tif")
   admin3_path = expanduser(joinpath("~", "data", "projects", "uganda2020", "outputs",
      "uga_canonical", "uga_canonical_names.shp"))
   admin1_path = expanduser(joinpath("~", "data", "projects", "uganda2020", "inputs",
      "uganda_districts_2019-wgs84", "uganda_districts_2019-wgs84.shp"))
   id_path = expanduser("~/data/projects/uganda2020/outputs/hrsl_admin_id.tif")

   @show ispath(admin3_path)
   ArchGDAL.read(hrsl_path) do hrsl_ds
      geo = band_geo(ArchGDAL.getgeotransform(hrsl_ds))
      hrsl_band = ArchGDAL.getband(hrsl_ds, 1)
      hrsl = ArchGDAL.read(hrsl_band)
      hrsl_no_data = ArchGDAL.getnodatavalue(hrsl_band)
      hrsl_raster = RasterRead(hrsl, hrsl_no_data, geo)
      id_data = ArchGDAL.read(id_path) do id_ds
         ArchGDAL.read(id_ds, 1)
      end
      id_raster = RasterRead(id_data, 0, geo)
      count_up = build_fractional(hrsl_raster, id_raster)

      ArchGDAL.read(admin3_path) do sc_dataset
         @show ArchGDAL.nlayer(sc_dataset)
         admin = ArchGDAL.getlayer(sc_dataset, 0)
         layerinfo(admin)
         totals = overenvelope(count_up, admin; limit = limit)
         writedlm("admin1_largest_fraction.txt", totals)
         largest = maximum(x[2] for x in totals)
         @show largest

         copy_path = "hrsl_proportion.tif"
         if ispath(copy_path)
            rm(copy_path)
         end
         ArchGDAL.copy(hrsl_ds, filename = copy_path) do hrsl_copy_ds
            ArchGDAL.write!(hrsl_copy_ds, hrsl, 1)
         end
      end
   end
   return nothing
end