text:
using Test, InMemoryDatasets, PooledArrays, Random, CategoricalArrays
const ≅ = isequal

# from DataFrames.jl
refds = Dataset(reshape(1.5:15.5, (3, 5)), :auto)

@testset "CartesianIndex" begin
    ds = Dataset(rand(2, 3), :auto)
    for i in axes(ds, 1), j in axes(ds, 2)
        @test ds[i, j] == ds[CartesianIndex(i, j)]
        r = rand()
        ds[CartesianIndex(i, j)] = r
        @test ds[i, j] == r
    end
    @test_throws BoundsError ds[CartesianIndex(0, 1)]
    @test_throws BoundsError ds[CartesianIndex(0, 0)]
    @test_throws BoundsError ds[CartesianIndex(1, 0)]
    @test_throws BoundsError ds[CartesianIndex(5, 1)]
    @test_throws BoundsError ds[CartesianIndex(5, 5)]
    @test_throws BoundsError ds[CartesianIndex(1, 5)]

    @test_throws BoundsError ds[CartesianIndex(0, 1)] = 1
    @test_throws ArgumentError ds[CartesianIndex(0, 0)] = 1
    @test_throws ArgumentError ds[CartesianIndex(1, 0)] = 1
    @test_throws BoundsError ds[CartesianIndex(5, 1)] = 1
    @test_throws ArgumentError ds[CartesianIndex(5, 5)] = 1
    @test_throws ArgumentError ds[CartesianIndex(1, 5)] = 1
end

@testset "broadcasting of AbstractDataset objects" begin
    for ds in (copy(refds), view(copy(refds), :, :))
        @test identity.(ds) == refds
        @test identity.(ds) !== ds
        @test (x->x).(ds) == refds
        @test (x->x).(ds) !== ds
        @test (ds .+ ds) ./ 2 == refds
        @test (ds .+ ds) ./ 2 !== ds
        @test ds .+ Matrix(ds) == 2 .* ds
        @test Matrix(ds) .+ ds == 2 .* ds
        @test (Matrix(ds) .+ ds .== 2 .* ds) == Dataset(trues(size(ds)), names(ds))
        @test ds .+ 1 == ds .+ ones(size(ds))
        @test ds .+ axes(ds, 1) == Dataset(Matrix(ds) .+ axes(ds, 1), names(ds))
        @test ds .+ permutedims(axes(ds, 2)) == Dataset(Matrix(ds) .+ permutedims(axes(ds, 2)), names(ds))
    end

    ds1 = copy(refds)
    ds2 = view(copy(refds), :, :)
    @test (ds1 .+ ds2) ./ 2 == refds
    @test (ds1 .- ds2) == Dataset(zeros(size(refds)), names(refds))
    @test (ds1 .* ds2) == refds .^ 2
    @test (ds1 ./ ds2) == Dataset(ones(size(refds)), names(refds))
end

@testset "broadcasting of AbstractDataset objects thrown exceptions" begin
    ds = copy(refds)
    dsv = view(ds, :, 2:ncol(ds))

    @test_throws DimensionMismatch ds .+ dsv
    @test_throws DimensionMismatch ds .+ ds[2:end, :]

    @test_throws DimensionMismatch ds .+ [1, 2]
    @test_throws DimensionMismatch ds .+ [1 2]
    @test_throws DimensionMismatch ds .+ rand(2, 2)
    @test_throws DimensionMismatch dsv .+ [1, 2]
    @test_throws DimensionMismatch dsv .+ [1 2]
    @test_throws DimensionMismatch dsv .+ rand(2, 2)

    ds2 = copy(ds)
    rename!(ds2, [:x1, :x2, :x3, :x4, :y])
    @test_throws ArgumentError ds .+ ds2
    @test_throws ArgumentError ds .+ 1 .+ ds2
end

@testset "broadcasting data sets" begin
    ds1 = Dataset(x=1, y=2)
    ds2 = Dataset(x=[1, 11], y=[2, 12])
    @test ds1 .+ ds2 == Dataset(x=[2, 12], y=[4, 14])

    ds1 = Dataset(x=1, y=2)
    ds2 = Dataset(x=[1, 11], y=[2, 12])
    x = ds2.x
    y = ds2.y
    t1 = IMD._get_lastmodified(IMD._attributes(ds2))
    sleep(.1)
    ds2 .+= ds1
    @test IMD._get_lastmodified(IMD._attributes(ds2)) != t1
    @test ds2.x === x
    @test ds2.y === y
    @test ds2 == Dataset(x=[2, 12], y=[4, 14])

    ds = Dataset(x=[1, 11], y=[2, 12])
    dsv = view(ds, 1:1, 1:2)
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)
    ds .-= dsv
    @test IMD._get_lastmodified(IMD._attributes(ds)) != t1
    @test ds == Dataset(x=[0, 10], y=[0, 10])

    @test Dataset() .+ Dataset() == Dataset()
    @test_throws ArgumentError Dataset(a=1, b=1) .+ Dataset(b=1, a=1)

    ds = Dataset(a=1, b=2)
    @test_throws ArgumentError ds .= Dataset(b=1, a=2)
    @test_throws ArgumentError ds .= Dataset(a=1, c=2)
    @test_throws ArgumentError ds[!, [:a, :b]] .= Dataset(b=1, a=2)
    @test_throws ArgumentError ds[!, [:a, :b]] .= Dataset(a=1, c=2)
end

@testset "broadcasting of AbstractDataset objects corner cases" begin
    ds = Dataset(c11=categorical(["a", "b"]), c12=categorical([missing, "b"]), c13=categorical(["a", missing]),
                   c21=categorical([1, 2]), c22=categorical([missing, 2]), c23=categorical([1, missing]),
                   p11=PooledArray(["a", "b"]), p12=PooledArray([missing, "b"]), p13=PooledArray(["a", missing]),
                   p21=PooledArray([1, 2]), p22=PooledArray([missing, 2]), p23=PooledArray([1, missing]),
                   b1=[true, false], b2=[missing, false], b3=[true, missing],
                   f1=[1.0, 2.0], f2=[missing, 2.0], f3=[1.0, missing],
                   s1=["a", "b"], s2=[missing, "b"], s3=["a", missing])

    ds2 = Dataset(c11=categorical(["a", "b"]), c12=[nothing, "b"], c13=["a", nothing],
                    c21=categorical([1, 2]), c22=[nothing, 2], c23=[1, nothing],
                    p11=["a", "b"], p12=[nothing, "b"], p13=["a", nothing],
                    p21=[1, 2], p22=[nothing, 2], p23=[1, nothing],
                    b1=[true, false], b2=[nothing, false], b3=[true, nothing],
                    f1=[1.0, 2.0], f2=[nothing, 2.0], f3=[1.0, nothing],
                    s1=["a", "b"], s2=[nothing, "b"], s3=["a", nothing])

    @test ds ≅ identity.(ds)
    @test ds ≅ (x->x).(ds)
    ds3 = coalesce.(ds, nothing)
    @test ds2 == ds3
    @test eltype.(eachcol(ds2)) == eltype.(eachcol(ds3))
    for i in axes(ds, 2)
        @test typeof(ds2[!, i]) == typeof(ds3[!, i])
    end
    ds4 = (x -> ds[1, 1]).(ds)
    @test names(ds4) == names(ds)
    @test all(isa.(eachcol(ds4), Ref(CategoricalArray)))
    @test all(eachcol(ds4) .== Ref(categorical(["a", "a"])))

    ds5 = Dataset(x=Any[1, 2, 3], y=Any[1, 2.0, big(3)])
    @test identity.(ds5) == ds5
    @test (x->x).(ds5) == ds5
    @test ds5 .+ 1 == Dataset(Matrix(ds5) .+ 1, names(ds5))
    @test eltype.(eachcol(identity.(ds5))) == [Union{Missing, Int}, Union{Missing, BigFloat}]
    @test eltype.(eachcol((x->x).(ds5))) == [Union{Missing, Int}, Union{Missing, BigFloat}]
    @test eltype.(eachcol(ds5 .+ 1)) == [Union{Missing, Int}, Union{Missing, BigFloat}]
end

@testset "normal data set and data set row in broadcasted assignment - one column" begin
    ds = copy(refds)
    ds[!, 1] .+= 1
    @test ds.x1 == [2.5, 3.5, 4.5]
    @test ds[:, 2:end] == refds[:, 2:end]
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)

    dsv = @view ds[1:2, 2:end]
    dsv[!, 1] .+= 100
    @test IMD._get_lastmodified(IMD._attributes(ds)) != t1
    @test ds.x2 == [104.5, 105.5, 6.5]
    # reverse the performed operations
    ds.x1 -= [1, 1, 1]
    ds.x2 -= [100, 100, 0]
    @test ds == Dataset(reshape(1.5:15.5, (3, 5)), :auto)

    ds = copy(refds)
    ds[:, 1] .+= 1
    @test ds.x1 == [2.5, 3.5, 4.5]
    @test ds[:, 2:end] == refds[:, 2:end]

    dsv = @view ds[1:2, 2:end]
    dsv.x2 .+= 1
    @test dsv.x2 == [5.5, 6.5]
    @test dsv[:, 2:end] == refds[1:2, 3:end]
    @test Matrix(ds) == [2.5  5.5  7.5  10.5  13.5
                         3.5  6.5  8.5  11.5  14.5
                         4.5  6.5  9.5  12.5  15.5]

    dsr = ds[1, 3:end]
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)

    dsr[end-1:end] .= 10
    @test t1 != IMD._get_lastmodified(IMD._attributes(ds))

    @test Vector(dsr) == [7.5, 10.0, 10.0]
    @test Matrix(ds) == [2.5  5.5  7.5  10.0  10.0
                         3.5  6.5  8.5  11.5  14.5
                         4.5  6.5  9.5  12.5  15.5]

    ds = copy(refds)
    ds[!, 1] .+= [1, 2, 3]
    @test ds.x1 == [2.5, 4.5, 6.5]
    @test ds[:, 2:end] == refds[:, 2:end]

    ds = copy(refds)
    ds[!, :x1] .+= [1, 2, 3]
    @test ds.x1.val == [2.5, 4.5, 6.5]
    @test ds[:, 2:end] == refds[:, 2:end]

    dsv = @view ds[1:2, 2:end]
    dsv[!, :x2] .+= [1, 2]
    @test dsv.x2 == [5.5, 7.5]
    @test dsv[:, 2:end] == refds[1:2, 3:end]
    @test Matrix(ds) == [2.5  5.5  7.5  10.5  13.5
                         4.5  7.5  8.5  11.5  14.5
                         6.5  6.5  9.5  12.5  15.5]

    dsr = ds[1, 3:end]
    dsr[end-1:end] .= [10, 11]
    @test Vector(dsr) == [7.5, 10.0, 11.0]
    @test Matrix(ds) == [2.5  5.5  7.5  10.0  11.0
                         4.5  7.5  8.5  11.5  14.5
                         6.5  6.5  9.5  12.5  15.5]

    ds = copy(refds)
    ds[:, 1] .+= [1, 2, 3]
    @test ds.x1 == [2.5, 4.5, 6.5]
    @test ds[:, 2:end] == refds[:, 2:end]
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)
    dsv = @view ds[1:2, 2:end]
    dsv[:, 1] .+= [1, 2]
    @test t1 != IMD._get_lastmodified(IMD._attributes(ds))
    @test dsv.x2 == [5.5, 7.5]
    @test dsv[:, 2:end] == refds[1:2, 3:end]
    @test Matrix(ds) == [2.5  5.5  7.5  10.5  13.5
                         4.5  7.5  8.5  11.5  14.5
                         6.5  6.5  9.5  12.5  15.5]

    # test a more complex broadcasting pattern
    ds = copy(refds)
    ds[!, 1] .+= [0, 1, 2] .+ 1
    @test ds.x1 == ds[!, 1] == [2.5, 4.5, 6.5]
    @test ds[:, 2:end] == refds[:, 2:end]

    ds = copy(refds)
    ds[!, "x1"] .+= [0, 1, 2] .+ 1
    @test ds."x1" == ds[!, 1] == [2.5, 4.5, 6.5]
    @test ds[:, Not("x1")] == refds[:, 2:end]

    ds = copy(refds)
    dsv = @view ds[1:2, 2:end]
    dsv[!, 1] .+= [0, 1] .+ 1
    @test ds == Dataset([1.5  5.5  7.5  10.5  13.5
                           2.5  7.5  8.5  11.5  14.5
                           3.5  6.5  9.5  12.5  15.5], :auto)

    dsv = @view ds[1:2, 2:end]
    @test_throws ArgumentError dsv[!, "x1"] .+= [0, 1] .+ 1
    @test ds == Dataset([1.5  5.5  7.5  10.5  13.5
                           2.5  7.5  8.5  11.5  14.5
                           3.5  6.5  9.5  12.5  15.5], :auto)

    ds = copy(refds)
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)
    ds[!, :x1] .+= [0, 1, 2] .+ 1
    @test t1 != IMD._get_lastmodified(IMD._attributes(ds))
    @test ds.x1 == [2.5, 4.5, 6.5]
    @test ds[:, 2:end] == refds[:, 2:end]

    dsv = @view ds[1:2, 2:end]
    dsv[!, :x2] .+= [0, 1] .+ 1
    @test dsv.x2 == [5.5, 7.5]
    @test dsv[:, 2:end] == refds[1:2, 3:end]
    @test Matrix(ds) == [2.5  5.5  7.5  10.5  13.5
                         4.5  7.5  8.5  11.5  14.5
                         6.5  6.5  9.5  12.5  15.5]

    ds = copy(refds)
    ds[!, "x1"] .+= [0, 1, 2] .+ 1
    @test ds."x1" == [2.5, 4.5, 6.5]
    @test ds[:, 2:end] == refds[:, 2:end]

    dsv = @view ds[1:2, 2:end]
    dsv[!, "x2"] .+= [0, 1] .+ 1
    @test dsv."x2" == [5.5, 7.5]
    @test dsv[:, 2:end] == refds[1:2, 3:end]
    @test Matrix(ds) == [2.5  5.5  7.5  10.5  13.5
                         4.5  7.5  8.5  11.5  14.5
                         6.5  6.5  9.5  12.5  15.5]

    dsr = ds[1, 3:end]
    dsr[end-1:end] .= [9, 10] .+ 1
    @test Vector(dsr) == [7.5, 10.0, 11.0]
    @test Matrix(ds) == [2.5  5.5  7.5  10.0  11.0
                         4.5  7.5  8.5  11.5  14.5
                         6.5  6.5  9.5  12.5  15.5]

    ds = copy(refds)
    ds[:, 1] .+= [0, 1, 2] .+ 1
    @test ds.x1 == [2.5, 4.5, 6.5]
    @test ds[:, 2:end] == refds[:, 2:end]

    dsv = @view ds[1:2, 2:end]
    dsv[:, 1] .+= [0, 1] .+ 1
    @test dsv.x2 == [5.5, 7.5]
    @test dsv[:, 2:end] == refds[1:2, 3:end]
    @test Matrix(ds) == [2.5  5.5  7.5  10.5  13.5
                         4.5  7.5  8.5  11.5  14.5
                         6.5  6.5  9.5  12.5  15.5]

    ds = copy(refds)
    ds[:, "x1"] .+= [0, 1, 2] .+ 1
    @test ds."x1" == [2.5, 4.5, 6.5]
    @test ds[:, 2:end] == refds[:, 2:end]

    dsv = @view ds[1:2, 2:end]
    dsv[:, "x2"] .+= [0, 1] .+ 1
    @test dsv."x2" == [5.5, 7.5]
    @test dsv[:, 2:end] == refds[1:2, 3:end]
    @test Matrix(ds) == [2.5  5.5  7.5  10.5  13.5
                         4.5  7.5  8.5  11.5  14.5
                         6.5  6.5  9.5  12.5  15.5]

    ds = copy(refds)
    dsv = @view ds[1:2, 2:end]
    dsr = ds[1, 3:end]
    @test_throws DimensionMismatch dsv[!, 1] .= fill(100, 2, 2)
    @test_throws DimensionMismatch dsv[!, 1] .= reshape(fill(200, 4), :, 2)
    @test_throws DimensionMismatch ds[!, 1] .= rand(1, 2)
    @test_throws DimensionMismatch dsr[end-1:end] .= rand(3, 1)
    @test_throws DimensionMismatch ds[:, 1] .= rand(1, 3)
    @test_throws DimensionMismatch dsv[:, 1] .= rand(1, 2)
    @test_throws DimensionMismatch ds[!, 1] .= reshape(rand(3), 1, :)
    @test_throws DimensionMismatch dsr[end-1:end] .= reshape(rand(3), :, 1)
    @test_throws DimensionMismatch ds[:, 1] .= reshape(rand(3), 1, :, 1)
    @test_throws DimensionMismatch dsv[:, 1] .= reshape(rand(2), 1, :, 1)

    ds = copy(refds)
    ds[!, :x1] .+= 1
    @test ds.x1 == [2.5, 3.5, 4.5]
    @test ds[:, 2:end] == refds[:, 2:end]

    dsv = @view ds[1:2, 2:end]
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)
    dsv[!, :x2] .+= 1
    @test t1 != IMD._get_lastmodified(IMD._attributes(ds))
    dsr = ds[1, 3:end]
    dsr[[:x4, :x5]] .= 10
    @test Vector(dsr) == [7.5, 10.0, 10.0]
    @test Matrix(ds) == [2.5  5.5  7.5  10.0  10.0
                         3.5  6.5  8.5  11.5  14.5
                         4.5  6.5  9.5  12.5  15.5]

    ds = copy(refds)
    ds[!, "x1"] .+= 1
    @test ds."x1" == [2.5, 3.5, 4.5]
    @test ds[:, 2:end] == refds[:, 2:end]

    dsv = @view ds[1:2, 2:end]
    dsv[!, "x2"] .+= 1

    dsr = ds[1, 3:end]
    dsr[["x4", "x5"]] .= 10
    @test Vector(dsr) == [7.5, 10.0, 10.0]
    @test Matrix(ds) == [2.5  5.5  7.5  10.0  10.0
                         3.5  6.5  8.5  11.5  14.5
                         4.5  6.5  9.5  12.5  15.5]

    ds = copy(refds)
    ds[:, :x1] .+= 1
    @test ds.x1 == [2.5, 3.5, 4.5]
    @test ds[:, 2:end] == refds[:, 2:end]

    dsv = @view ds[1:2, 2:end]
    dsv[:, :x2] .+= 1
    @test dsv.x2 == [5.5, 6.5]
    @test dsv[:, 2:end] == refds[1:2, 3:end]
    @test Matrix(ds) == [2.5  5.5  7.5  10.5  13.5
                         3.5  6.5  8.5  11.5  14.5
                         4.5  6.5  9.5  12.5  15.5]

    ds = copy(refds)
    ds[:, "x1"] .+= 1
    @test ds."x1" == [2.5, 3.5, 4.5]
    @test ds[:, 2:end] == refds[:, 2:end]

    dsv = @view ds[1:2, 2:end]
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)
    dsv[:, "x2"] .+= 1
    @test t1 != IMD._get_lastmodified(IMD._attributes(ds))
    @test dsv."x2" == [5.5, 6.5]
    @test dsv[:, 2:end] == refds[1:2, 3:end]
    @test Matrix(ds) == [2.5  5.5  7.5  10.5  13.5
                         3.5  6.5  8.5  11.5  14.5
                         4.5  6.5  9.5  12.5  15.5]

    ds = copy(refds)
    ds[!, :x1] .+= [1, 2, 3]
    @test ds.x1 == [2.5, 4.5, 6.5]
    @test ds[:, 2:end] == refds[:, 2:end]

    dsv = @view ds[1:2, 2:end]
    dsv[!, :x2] .+= [1, 2]

    dsr = ds[1, 3:end]
    dsr[[:x4, :x5]] .= [10, 11]
    @test Vector(dsr) == [7.5, 10.0, 11.0]
    @test Matrix(ds) == [2.5  5.5  7.5  10.0  11.0
                         4.5  7.5  8.5  11.5  14.5
                         6.5  6.5  9.5  12.5  15.5]

    ds = copy(refds)
    ds[!, "x1"] .+= [1, 2, 3]
    @test ds."x1" == [2.5, 4.5, 6.5]
    @test ds[:, 2:end] == refds[:, 2:end]

    dsv = @view ds[1:2, 2:end]
    dsv[!, :x2] .+= [1, 2]

    dsr = ds[1, 3:end]
    dsr[["x4", "x5"]] .= [10, 11]
    @test Vector(dsr) == [7.5, 10.0, 11.0]
    @test Matrix(ds) == [2.5  5.5  7.5  10.0  11.0
                         4.5  7.5  8.5  11.5  14.5
                         6.5  6.5  9.5  12.5  15.5]

    ds = copy(refds)
    ds[:, :x1] .+= [1, 2, 3]
    @test ds.x1 == [2.5, 4.5, 6.5]
    @test ds[:, 2:end] == refds[:, 2:end]

    dsv = @view ds[1:2, 2:end]
    dsv[:, :x2] .+= [1, 2]
    @test dsv.x2 == [5.5, 7.5]
    @test dsv[:, 2:end] == refds[1:2, 3:end]
    @test Matrix(ds) == [2.5  5.5  7.5  10.5  13.5
                         4.5  7.5  8.5  11.5  14.5
                         6.5  6.5  9.5  12.5  15.5]

    ds = copy(refds)
    ds[:, "x1"] .+= [1, 2, 3]
    @test ds."x1" == [2.5, 4.5, 6.5]
    @test ds[:, 2:end] == refds[:, 2:end]

    dsv = @view ds[1:2, 2:end]
    dsv[:, "x2"] .+= [1, 2]
    @test dsv."x2" == [5.5, 7.5]
    @test dsv[:, 2:end] == refds[1:2, 3:end]
    @test Matrix(ds) == [2.5  5.5  7.5  10.5  13.5
                         4.5  7.5  8.5  11.5  14.5
                         6.5  6.5  9.5  12.5  15.5]

    ds = copy(refds)
    dsv = @view ds[1:2, 2:end]
    dsr = ds[1, 3:end]

    @test_throws DimensionMismatch dsv[!, :x2] .= fill(100, 2, 2)
    @test_throws DimensionMismatch dsv[!, 1] .= reshape(fill(200, 4), :, 2)
    @test_throws DimensionMismatch dsv[!, "x2"] .= fill(100, 2, 2)
    @test_throws DimensionMismatch ds[!, :x1] .= rand(1, 3)
    @test_throws DimensionMismatch dsr[[:x4, :x5]] .= rand(3, 1)
    @test_throws DimensionMismatch ds[:, :x1] .= rand(1, 3)
    @test_throws DimensionMismatch dsv[:, :x2] .= rand(1, 2)
    @test_throws DimensionMismatch ds[!, 1] .= reshape(rand(3), 1, :)
    @test_throws DimensionMismatch dsr[end-1:end] .= reshape(rand(3), :, 1)
    @test_throws DimensionMismatch ds[:, 1] .= reshape(rand(3), 1, :)
    @test_throws DimensionMismatch dsv[:, 1] .= reshape(rand(2), 1, :)
    @test_throws DimensionMismatch ds[!, "x1"] .= rand(1, 3)
    @test_throws DimensionMismatch dsr[["x4", "x5"]] .= rand(3, 1)
    @test_throws DimensionMismatch ds[:, "x1"] .= rand(1, 3)
    @test_throws DimensionMismatch dsv[:, "x2"] .= rand(1, 2)
end

@testset "normal data set and data set view in broadcasted assignment - two columns" begin
    ds = copy(refds)
    ds[:, [1, 2]] .= Matrix(ds[:, [1, 2]]) .+ 1
    @test ds.x1 == [2.5, 3.5, 4.5]
    @test ds.x2 == [5.5, 6.5, 7.5]
    @test ds[:, 3:end] == refds[:, 3:end]

    dsv = @view ds[1:2, 3:end]
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)
    dsv[:, [1, 2]] .= Matrix(dsv[:, [1, 2]]) .+ 1
    @test t1 != IMD._get_lastmodified(IMD._attributes(ds))

    @test dsv.x3 == [8.5, 9.5]
    @test dsv.x4 == [11.5, 12.5]
    @test dsv[:, 3:end] == refds[1:2, 5:end]
    @test Matrix(ds) == [2.5  5.5  8.5  11.5  13.5
                         3.5  6.5  9.5  12.5  14.5
                         4.5  7.5  9.5  12.5  15.5]

    ds = copy(refds)
    ds[:, [1, 2]] .= Matrix(ds[:, [1, 2]]) .+ 1
    @test ds.x1 == [2.5, 3.5, 4.5]
    @test ds.x2 == [5.5, 6.5, 7.5]
    @test ds[:, 3:end] == refds[:, 3:end]

    dsv = @view ds[1:2, 3:end]
    dsv[:, [1, 2]] .= Matrix(dsv[:, [1, 2]]) .+ 1
    @test dsv.x3 == [8.5, 9.5]
    @test dsv.x4 == [11.5, 12.5]
    @test dsv[:, 3:end] == refds[1:2, 5:end]
    @test Matrix(ds) == [2.5  5.5  8.5  11.5  13.5
                         3.5  6.5  9.5  12.5  14.5
                         4.5  7.5  9.5  12.5  15.5]

    ds = copy(refds)
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)

    ds[:, [1, 2]] .= Matrix(ds[:, [1, 2]]) .+ [1 4
                                               2 5
                                               3 6]
   @test t1 != IMD._get_lastmodified(IMD._attributes(ds))

    @test ds.x1 == [2.5, 4.5, 6.5]
    @test ds.x2 == [8.5, 10.5, 12.5]
    @test ds[:, 3:end] == refds[:, 3:end]

    dsv = @view ds[1:2, 3:end]
    dsv[:, [1, 2]] .= Matrix(dsv[:, [1, 2]]) .+ [1 3
                                               2 4]
    @test dsv.x3 == [8.5, 10.5]
    @test dsv.x4 == [13.5, 15.5]
    @test dsv[:, 3:end] == refds[1:2, 5:end]
    @test Matrix(ds) == [2.5   8.5   8.5  13.5  13.5
                         4.5  10.5  10.5  15.5  14.5
                         6.5  12.5   9.5  12.5  15.5]

    ds = copy(refds)
    ds[:, [1, 2]] .= Matrix(ds[:, [1, 2]]) .+ [1 4
                                               2 5
                                               3 6]
    @test ds.x1 == [2.5, 4.5, 6.5]
    @test ds.x2 == [8.5, 10.5, 12.5]
    @test ds[:, 3:end] == refds[:, 3:end]

    dsv = @view ds[1:2, 3:end]
    dsv[:, [1, 2]] .= Matrix(dsv[:, [1, 2]]) .+ [1 3
                                               2 4]
    @test dsv.x3 == [8.5, 10.5]
    @test dsv.x4 == [13.5, 15.5]
    @test dsv[:, 3:end] == refds[1:2, 5:end]
    @test Matrix(ds) == [2.5   8.5   8.5  13.5  13.5
                         4.5  10.5  10.5  15.5  14.5
                         6.5  12.5   9.5  12.5  15.5]

    ds = copy(refds)
    dsv = @view ds[1:2, 2:end]
    @test_throws DimensionMismatch ds[:, [1, 2]] .= rand(3, 10)
    @test_throws DimensionMismatch dsv[:, [1, 2]] .= rand(2, 10)
    @test_throws DimensionMismatch ds[:, [1, 2]] .= rand(3, 10)
    @test_throws DimensionMismatch dsv[:, [1, 2]] .= rand(2, 10)

    ds = copy(refds)
    ds[:, [:x1, :x2]] .= Matrix(ds[:, [:x1, :x2]]) .+ 1
    @test ds.x1 == [2.5, 3.5, 4.5]
    @test ds.x2 == [5.5, 6.5, 7.5]
    @test ds[:, 3:end] == refds[:, 3:end]

    dsv = @view ds[1:2, 3:end]
    dsv[:, [:x3, :x4]] .= Matrix(dsv[:, [:x3, :x4]]) .+ 1
    @test dsv.x3 == [8.5, 9.5]
    @test dsv.x4 == [11.5, 12.5]
    @test dsv[:, 3:end] == refds[1:2, 5:end]
    @test Matrix(ds) == [2.5  5.5  8.5  11.5  13.5
                         3.5  6.5  9.5  12.5  14.5
                         4.5  7.5  9.5  12.5  15.5]

    ds = copy(refds)
    ds[:, ["x1", "x2"]] .= Matrix(ds[:, [:x1, :x2]]) .+ 1
    @test ds.x1 == [2.5, 3.5, 4.5]
    @test ds.x2 == [5.5, 6.5, 7.5]
    @test ds[:, 3:end] == refds[:, 3:end]

    dsv = @view ds[1:2, 3:end]
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)

    dsv[:, ["x3", "x4"]] .= Matrix(dsv[:, [:x3, :x4]]) .+ 1
    @test t1 != IMD._get_lastmodified(IMD._attributes(ds))

    @test dsv.x3 == [8.5, 9.5]
    @test dsv.x4 == [11.5, 12.5]
    @test dsv[:, 3:end] == refds[1:2, 5:end]
    @test Matrix(ds) == [2.5  5.5  8.5  11.5  13.5
                         3.5  6.5  9.5  12.5  14.5
                         4.5  7.5  9.5  12.5  15.5]

    ds = copy(refds)
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)

    ds[:, [:x1, :x2]] .= Matrix(ds[:, [:x1, :x2]]) .+ 1
    @test t1 != IMD._get_lastmodified(IMD._attributes(ds))

    @test ds.x1 == [2.5, 3.5, 4.5]
    @test ds.x2 == [5.5, 6.5, 7.5]
    @test ds[:, 3:end] == refds[:, 3:end]

    dsv = @view ds[1:2, 3:end]
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)

    dsv[:, [:x3, :x4]] .= Matrix(dsv[:, [:x3, :x4]]) .+ 1
    @test t1 != IMD._get_lastmodified(IMD._attributes(ds))

    @test dsv.x3 == [8.5, 9.5]
    @test dsv.x4 == [11.5, 12.5]
    @test dsv[:, 3:end] == refds[1:2, 5:end]
    @test Matrix(ds) == [2.5  5.5  8.5  11.5  13.5
                         3.5  6.5  9.5  12.5  14.5
                         4.5  7.5  9.5  12.5  15.5]

    ds = copy(refds)
    ds[:, ["x1", "x2"]] .= Matrix(ds[:, [:x1, :x2]]) .+ 1
    @test ds.x1 == [2.5, 3.5, 4.5]
    @test ds.x2 == [5.5, 6.5, 7.5]
    @test ds[:, 3:end] == refds[:, 3:end]

    dsv = @view ds[1:2, 3:end]
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)

    dsv[:, ["x3", "x4"]] .= Matrix(dsv[:, [:x3, :x4]]) .+ 1
    @test t1 != IMD._get_lastmodified(IMD._attributes(ds))

    @test dsv.x3 == [8.5, 9.5]
    @test dsv.x4 == [11.5, 12.5]
    @test dsv[:, 3:end] == refds[1:2, 5:end]
    @test Matrix(ds) == [2.5  5.5  8.5  11.5  13.5
                         3.5  6.5  9.5  12.5  14.5
                         4.5  7.5  9.5  12.5  15.5]

    ds = copy(refds)
    ds[:, [:x1, :x2]] .= Matrix(ds[:, [:x1, :x2]]) .+ [1 4
                                                       2 5
                                                       3 6]
    @test ds.x1 == [2.5, 4.5, 6.5]
    @test ds.x2 == [8.5, 10.5, 12.5]
    @test ds[:, 3:end] == refds[:, 3:end]

    dsv = @view ds[1:2, 3:end]
    dsv[:, [:x3, :x4]] .= Matrix(dsv[:, [:x3, :x4]]) .+ [1 3
                                                         2 4]
    @test dsv.x3 == [8.5, 10.5]
    @test dsv.x4 == [13.5, 15.5]
    @test dsv[:, 3:end] == refds[1:2, 5:end]
    @test Matrix(ds) == [2.5   8.5   8.5  13.5  13.5
                         4.5  10.5  10.5  15.5  14.5
                         6.5  12.5   9.5  12.5  15.5]

    ds = copy(refds)
    ds[:, ["x1", "x2"]] .= Matrix(ds[:, ["x1", "x2"]]) .+ [1 4
                                                           2 5
                                                           3 6]
    @test ds.x1 == [2.5, 4.5, 6.5]
    @test ds.x2 == [8.5, 10.5, 12.5]
    @test ds[:, 3:end] == refds[:, 3:end]

    dsv = @view ds[1:2, 3:end]
    dsv[:, ["x3", "x4"]] .= Matrix(dsv[:, ["x3", "x4"]]) .+ [1 3
                                                             2 4]
    @test dsv.x3 == [8.5, 10.5]
    @test dsv.x4 == [13.5, 15.5]
    @test dsv[:, 3:end] == refds[1:2, 5:end]
    @test Matrix(ds) == [2.5   8.5   8.5  13.5  13.5
                         4.5  10.5  10.5  15.5  14.5
                         6.5  12.5   9.5  12.5  15.5]

    ds = copy(refds)
    ds[:, [:x1, :x2]] .= Matrix(ds[:, [:x1, :x2]]) .+ [1 4
                                                       2 5
                                                       3 6]
    @test ds.x1 == [2.5, 4.5, 6.5]
    @test ds.x2 == [8.5, 10.5, 12.5]
    @test ds[:, 3:end] == refds[:, 3:end]

    dsv = @view ds[1:2, 3:end]
    dsv[:, [:x3, :x4]] .= Matrix(dsv[:, [:x3, :x4]]) .+ [1 3
                                                         2 4]
    @test dsv.x3 == [8.5, 10.5]
    @test dsv.x4 == [13.5, 15.5]
    @test dsv[:, 3:end] == refds[1:2, 5:end]
    @test Matrix(ds) == [2.5   8.5   8.5  13.5  13.5
                         4.5  10.5  10.5  15.5  14.5
                         6.5  12.5   9.5  12.5  15.5]

    ds = copy(refds)
    ds[:, ["x1", "x2"]] .= Matrix(ds[:, ["x1", "x2"]]) .+ [1 4
                                                           2 5
                                                           3 6]
    @test ds.x1 == [2.5, 4.5, 6.5]
    @test ds.x2 == [8.5, 10.5, 12.5]
    @test ds[:, 3:end] == refds[:, 3:end]

    dsv = @view ds[1:2, 3:end]
    dsv[:, ["x3", "x4"]] .= Matrix(dsv[:, ["x3", "x4"]]) .+ [1 3
                                                             2 4]
    @test dsv.x3 == [8.5, 10.5]
    @test dsv.x4 == [13.5, 15.5]
    @test dsv[:, 3:end] == refds[1:2, 5:end]
    @test Matrix(ds) == [2.5   8.5   8.5  13.5  13.5
                         4.5  10.5  10.5  15.5  14.5
                         6.5  12.5   9.5  12.5  15.5]

    ds = copy(refds)
    dsv = @view ds[1:2, 2:end]
    @test_throws DimensionMismatch ds[:, [:x1, :x2]] .= rand(3, 10)
    @test_throws DimensionMismatch dsv[:, [:x3, :x4]] .= rand(2, 10)
    @test_throws DimensionMismatch ds[:, [:x1, :x2]] .= rand(3, 10)
    @test_throws DimensionMismatch dsv[:, [:x3, :x4]] .= rand(2, 10)
    @test_throws DimensionMismatch ds[:, ["x1", "x2"]] .= rand(3, 10)
    @test_throws DimensionMismatch dsv[:, ["x3", "x4"]] .= rand(2, 10)
    @test_throws DimensionMismatch ds[:, ["x1", "x2"]] .= rand(3, 10)
    @test_throws DimensionMismatch dsv[:, ["x3", "x4"]] .= rand(2, 10)

    ds = copy(refds)
    ds[:, [1, 2]] .= [1 2
                      3 4
                      5 6]
    @test Matrix(ds) == [1.0  2.0  7.5  10.5  13.5
                         3.0  4.0  8.5  11.5  14.5
                         5.0  6.0  9.5  12.5  15.5]

    ds = copy(refds)
    ds[:, [1, 2]] .= [1, 3, 5]
    @test Matrix(ds) == [1.0  1.0  7.5  10.5  13.5
                         3.0  3.0  8.5  11.5  14.5
                         5.0  5.0  9.5  12.5  15.5]

    ds = copy(refds)
    ds[:, [1, 2]] .= reshape([1, 3, 5], 3, 1)
    @test Matrix(ds) == [1.0  1.0  7.5  10.5  13.5
                         3.0  3.0  8.5  11.5  14.5
                         5.0  5.0  9.5  12.5  15.5]

    ds = copy(refds)
    ds[:, [1, 2]] .= 1
    @test Matrix(ds) == [1.0  1.0  7.5  10.5  13.5
                         1.0  1.0  8.5  11.5  14.5
                         1.0  1.0  9.5  12.5  15.5]

    ds = copy(refds)
    dsv = view(ds, 2:3, 2:4)
    dsv[:, [1, 2]] .= [1 2
                       3 4]
    @test Matrix(ds) == [1.5  4.5  7.5  10.5  13.5
                         2.5  1.0  2.0  11.5  14.5
                         3.5  3.0  4.0  12.5  15.5]

    ds = copy(refds)
    dsv = view(ds, 2:3, 2:4)
    dsv[:, [1, 2]] .= [1, 3]
    @test Matrix(ds) == [1.5  4.5  7.5  10.5  13.5
                         2.5  1.0  1.0  11.5  14.5
                         3.5  3.0  3.0  12.5  15.5]

    ds = copy(refds)
    dsv = view(ds, 2:3, 2:4)
    dsv[:, [1, 2]] .= reshape([1, 3], 2, 1)
    @test Matrix(ds) == [1.5  4.5  7.5  10.5  13.5
                         2.5  1.0  1.0  11.5  14.5
                         3.5  3.0  3.0  12.5  15.5]

    ds = copy(refds)
    dsv = view(ds, 2:3, 2:4)
    dsv[:, [1, 2]] .= 1
    @test Matrix(ds) == [1.5  4.5  7.5  10.5  13.5
                         2.5  1.0  1.0  11.5  14.5
                         3.5  1.0  1.0  12.5  15.5]
end

@testset "assignment to a whole data set and data set row" begin
    ds = copy(refds)
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)

    ds .= 10
    @test t1 != IMD._get_lastmodified(IMD._attributes(ds))

    @test all(Matrix(ds) .== 10)
    dsv = view(ds, 1:2, 1:4)
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)

    dsv .= 100
    @test t1 != IMD._get_lastmodified(IMD._attributes(ds))

    @test Matrix(ds) == [100.0  100.0  100.0  100.0  10.0
                        100.0  100.0  100.0  100.0  10.0
                         10.0   10.0   10.0   10.0  10.0]
    dsr = ds[1, 1:2]
    dsr .= 1000
    @test Matrix(ds) == [1000.0  1000.0  100.0  100.0  10.0
                         100.0   100.0  100.0  100.0  10.0
                         10.0    10.0   10.0   10.0  10.0]

    ds = copy(refds)
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)

    ds[:, :] .= 10
    @test t1 != IMD._get_lastmodified(IMD._attributes(ds))

    @test all(Matrix(ds) .== 10)
    dsv = view(ds, 1:2, 1:4)
    dsv[:, :] .= 100
    @test Matrix(ds) == [100.0  100.0  100.0  100.0  10.0
                         100.0  100.0  100.0  100.0  10.0
                         10.0   10.0   10.0   10.0  10.0]
    dsr = ds[1, 1:2]
    dsr[:] .= 1000
    @test Matrix(ds) == [1000.0  1000.0  100.0  100.0  10.0
                         100.0   100.0   100.0  100.0  10.0
                         10.0    10.0    10.0   10.0   10.0]

    ds = copy(refds)
    ds[:, :] .= 10
    @test all(Matrix(ds) .== 10)
    dsv = view(ds, 1:2, 1:4)
    dsv[:, :] .= 100
    @test Matrix(ds) == [100.0  100.0  100.0  100.0  10.0
                         100.0  100.0  100.0  100.0  10.0
                         10.0   10.0   10.0   10.0  10.0]
end

@testset "extending data set in broadcasted assignment - one column" begin
    ds = copy(refds)
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)

    ds[!, :a] .= 1
    @test t1 != IMD._get_lastmodified(IMD._attributes(ds))

    @test Matrix(ds) == [1.5  4.5  7.5  10.5  13.5  1.0
                         2.5  5.5  8.5  11.5  14.5  1.0
                         3.5  6.5  9.5  12.5  15.5  1.0]
    @test names(ds)[end] == "a"
    @test ds[:, 1:end-1] == refds
    ds[!, :b] .= [1, 2, 3]
    @test Matrix(ds) == [1.5  4.5  7.5  10.5  13.5  1.0 1.0
                         2.5  5.5  8.5  11.5  14.5  1.0 2.0
                         3.5  6.5  9.5  12.5  15.5  1.0 3.0]
    @test names(ds)[end] == "b"
    @test ds[:, 1:end-2] == refds
    cds = copy(ds)
    @test_throws DimensionMismatch ds[!, :c] .= ones(1, 3)
    @test ds == cds
    @test_throws DimensionMismatch ds[!, :x] .= ones(4)
    @test ds == cds
    @test_throws ArgumentError ds[!, 10] .= ones(3)
    @test ds == cds

    dsv = @view ds[1:2, 2:end]
    @test_throws BoundsError dsv[!, 10] .= ones(3)
    @test_throws ArgumentError dsv[!, :z] .= ones(3)
    @test ds == cds
    dsr = ds[1, 3:end]
    @test_throws BoundsError dsr[10] .= ones(3)
    @test_throws ArgumentError dsr[:z] .= ones(3)
    @test ds == cds

    ds = Dataset()
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)

    @test_throws DimensionMismatch ds[!, :a] .= sin.(1:3)
    ds[!, :b] .= sin.(1)
    ds[!, :c] .= sin(1) .+ 1
    @test t1 != IMD._get_lastmodified(IMD._attributes(ds))

    @test ds == Dataset(b=Float64[], c=Float64[])

    ds = copy(refds)
    ds[!, "a"] .= 1
    @test Matrix(ds) == [1.5  4.5  7.5  10.5  13.5  1.0
                         2.5  5.5  8.5  11.5  14.5  1.0
                         3.5  6.5  9.5  12.5  15.5  1.0]
    @test names(ds)[end] == "a"
    @test ds[:, 1:end-1] == refds
    ds[!, "b"] .= [1, 2, 3]
    @test Matrix(ds) == [1.5  4.5  7.5  10.5  13.5  1.0 1.0
                         2.5  5.5  8.5  11.5  14.5  1.0 2.0
                         3.5  6.5  9.5  12.5  15.5  1.0 3.0]
    @test names(ds)[end] == "b"
    @test ds[:, 1:end-2] == refds
    cds = copy(ds)
    @test_throws DimensionMismatch ds[!, "c"] .= ones(1, 3)
    @test ds == cds
    @test_throws DimensionMismatch ds[!, "x"] .= ones(4)
    @test ds == cds
    @test_throws ArgumentError ds[!, 10] .= ones(3)
    @test ds == cds

    dsv = @view ds[1:2, 2:end]
    @test_throws BoundsError dsv[!, 10] .= ones(3)
    @test_throws ArgumentError dsv[!, "z"] .= ones(3)
    @test ds == cds
    dsr = ds[1, 3:end]
    @test_throws BoundsError dsr[10] .= ones(3)
    @test_throws ArgumentError dsr["z"] .= ones(3)
    @test ds == cds

    ds = Dataset()
    @test_throws DimensionMismatch ds[!, "a"] .= sin.(1:3)
    ds[!, "b"] .= sin.(1)
    ds[!, "c"] .= sin(1) .+ 1
    @test ds == Dataset(b=Float64[], c=Float64[])
end

@testset "empty data frame corner case" begin
    ds = Dataset()
    @test_throws ArgumentError ds[!, 1] .= 1
    @test_throws ArgumentError ds[!, 2] .= 1
    @test_throws ArgumentError ds[!, [:a, :b]] .= [1]
    @test_throws ArgumentError ds[!, [:a, :b]] .= 1
    @test_throws DimensionMismatch ds[!, :a] .= [1 2]
    @test_throws DimensionMismatch ds[!, :a] .= [1, 2]
    @test_throws DimensionMismatch ds[!, :a] .= sin.(1) .+ [1, 2]
    @test_throws ArgumentError ds[!, ["a", "b"]] .= [1]
    @test_throws ArgumentError ds[!, ["a", "b"]] .= 1
    @test_throws DimensionMismatch ds[!, "a"] .= [1 2]
    @test_throws DimensionMismatch ds[!, "a"] .= [1, 2]
    @test_throws DimensionMismatch ds[!, "a"] .= sin.(1) .+ [1, 2]

    for rhs in [1, [1], Int[], "abc", ["abc"]]
        ds = Dataset()
        t1 = IMD._get_lastmodified(IMD._attributes(ds))
        sleep(.1)

        ds[!, :a] .= rhs
        @test size(ds) == (0, 1)
        @test nonmissingtype(eltype(ds[!, 1])) == (rhs isa AbstractVector ? eltype(rhs) : nonmissingtype(typeof(rhs)))
        @test t1 != IMD._get_lastmodified(IMD._attributes(ds))

        ds = Dataset()
        t1 = IMD._get_lastmodified(IMD._attributes(ds))
        sleep(.1)

        ds[!, :a] .= length.(rhs)
        @test size(ds) == (0, 1)
        @test eltype(ds[!, 1]) == Union{Missing, Int}
        @test t1 != IMD._get_lastmodified(IMD._attributes(ds))

        ds = Dataset()
        ds[!, :a] .= length.(rhs) .+ 1
        @test size(ds) == (0, 1)
        @test eltype(ds[!, 1]) == Union{Missing, Int}

        ds = Dataset()
        @. ds[!, :a] = length(rhs) + 1
        @test size(ds) == (0, 1)
        @test eltype(ds[!, 1]) == Union{Missing, Int}

        ds = Dataset(x=Int[])
        ds[!, :a] .= rhs
        @test size(ds) == (0, 2)
        @test nonmissingtype(eltype(ds[!, 2])) == (rhs isa AbstractVector ? eltype(rhs) : nonmissingtype(typeof(rhs)))

        ds = Dataset(x=Int[])
        ds[!, :a] .= length.(rhs)
        @test size(ds) == (0, 2)
        @test eltype(ds[!, 2]) == Union{Missing, Int}

        ds = Dataset(x=Int[])
        ds[!, :a] .= length.(rhs) .+ 1
        @test size(ds) == (0, 2)
        @test eltype(ds[!, 2]) == Union{Missing, Int}

        ds = Dataset(x=Int[])
        @. ds[!, :a] = length(rhs) + 1
        @test size(ds) == (0, 2)
        @test eltype(ds[!, 2]) == Union{Missing, Int}

        ds = Dataset()
        ds[!, "a"] .= rhs
        @test size(ds) == (0, 1)
        @test nonmissingtype(eltype(ds[!, 1])) == (rhs isa AbstractVector ? eltype(rhs) : nonmissingtype(typeof(rhs)))

        ds = Dataset()
        ds[!, "a"] .= length.(rhs)
        @test size(ds) == (0, 1)
        @test eltype(ds[!, 1]) == Union{Missing, Int}

        ds = Dataset()
        ds[!, "a"] .= length.(rhs) .+ 1
        @test size(ds) == (0, 1)
        @test eltype(ds[!, 1]) == Union{Missing, Int}

        ds = Dataset()
        @. ds[!, "a"] = length(rhs) + 1
        @test size(ds) == (0, 1)
        @test eltype(ds[!, 1]) == Union{Missing, Int}

        ds = Dataset(x=Int[])
        ds[!, "a"] .= rhs
        @test size(ds) == (0, 2)
        @test nonmissingtype(eltype(ds[!, 2])) == (rhs isa AbstractVector ? eltype(rhs) : nonmissingtype(typeof(rhs)))

        ds = Dataset(x=Int[])
        ds[!, "a"] .= length.(rhs)
        @test size(ds) == (0, 2)
        @test eltype(ds[!, 2]) == Union{Missing, Int}

        ds = Dataset(x=Int[])
        ds[!, "a"] .= length.(rhs) .+ 1
        @test size(ds) == (0, 2)
        @test eltype(ds[!, 2]) == Union{Missing, Int}

        ds = Dataset(x=Int[])
        @. ds[!, "a"] = length(rhs) + 1
        @test size(ds) == (0, 2)
        @test eltype(ds[!, 2]) == Union{Missing, Int}
    end

    ds = Dataset()
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)

    ds .= 1
    @test t1 != IMD._get_lastmodified(IMD._attributes(ds))

    @test ds == Dataset()
    ds .= [1]
    @test ds == Dataset()
    ds .= ones(1, 1)
    @test ds == Dataset()
    @test_throws DimensionMismatch ds .= ones(1, 2)
    @test_throws DimensionMismatch ds .= ones(1, 2, 1)

    ds = Dataset(a=[])
    ds[!, :b] .= sin.(1)
    @test eltype(ds.b) == Union{Missing, Float64}
    ds[!, :b] .= [1]
    @test eltype(ds.b) == Union{Missing, Int}
    ds[!, :b] .= 'a'
    @test eltype(ds.b) == Union{Missing, Char}
    @test names(ds) == ["a", "b"]

    c = categorical(["a", "b", "c"])
    ds = Dataset()
    @test_throws DimensionMismatch ds[!, :a] .= c

    ds[!, :b] .= c[1]
    @test nrow(ds) == 0
    @test ds.b.val isa CategoricalVector{Union{Missing, String}, UInt32, String, CategoricalValue{String, UInt32}, Missing}

    ds = Dataset(a=[])
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)

    ds[!, "b"] .= sin.(1)
    @test eltype(ds."b") == Union{Missing, Float64}
    @test t1 != IMD._get_lastmodified(IMD._attributes(ds))

    ds[!, "b"] .= [1]
    @test eltype(ds."b") == Union{Int, Missing}
    ds[!, "b"] .= 'a'
    @test eltype(ds."b") == Union{Char, Missing}
    @test names(ds) == ["a", "b"]

    c = categorical(["a", "b", "c"])
    ds = Dataset()
    @test_throws DimensionMismatch ds[!, "a"] .= c

    ds[!, "b"] .= c[1]
    @test nrow(ds) == 0
    @test ds."b".val isa CategoricalVector{Union{Missing, String}}
end

@testset "test categorical values" begin
    for v in Any[categorical([1, 2, 3]), categorical([1, 2, missing]),
              categorical([missing, 1, 2]),
              categorical(["1", "2", "3"]), categorical(["1", "2", missing]),
              categorical([missing, "1", "2"])]
        ds = copy(refds)
        t1 = IMD._get_lastmodified(IMD._attributes(ds))
        sleep(.1)

        ds[!, :c1] .= v
        @test t1 != IMD._get_lastmodified(IMD._attributes(ds))

        @test ds.c1 ≅ v
        @test ds.c1.val !== v
        @test ds.c1.val isa CategoricalVector
        @test levels(ds.c1.val) == levels(v)
        @test levels(ds.c1.val) !== levels(v)
        ds[!, :c2] .= v[2]
        @test ds.c2.val == fill(v[2], 3)
        @test ds.c2.val isa CategoricalVector
        @test levels(ds.c2.val) == levels(v)
        ds[!, :c3] .= (x->x).(v)
        @test ds.c3.val ≅ v
        @test ds.c3.val !== v
        @test ds.c3.val isa CategoricalVector
        @test levels(ds.c3.val) == levels(v)
        @test levels(ds.c3.val) !== levels(v)
        ds[!, :c4] .= identity.(v)
        @test ds.c4.val ≅ v
        @test ds.c4.val !== v
        @test ds.c4.val isa CategoricalVector
        @test levels(ds.c4.val) == levels(v)
        @test levels(ds.c4.val) !== levels(v)
        ds[!, :c5] .= (x->v[2]).(v)
        @test unique(ds.c5.val) == [unwrap(v[2])]
        @test ds.c5.val isa CategoricalVector
        @test levels(ds.c5.val) == levels(v)
    end
end

@testset "scalar broadcasting" begin
    a = Dataset(x=zeros(2))
    a .= 1 ./ (1 + 2)
    @test a.x == [1/3, 1/3]
    a .= 1 ./ (1 .+ 3)
    @test a.x == [1/4, 1/4]
    a .= sqrt.(1 ./ 2)
    @test a.x == [sqrt(1/2), sqrt(1/2)]
end

@testset "tuple broadcasting" begin
    X = Dataset(zeros(2, 3), :auto)
    X .= (1, 2)
    @test X == Dataset([1 1 1; 2 2 2], :auto)

    X = Dataset(zeros(2, 3), :auto)
    X .= (1, 2) .+ 10 .- X
    @test X == Dataset([11 11 11; 12 12 12], :auto)

    X = Dataset(zeros(2, 3), :auto)
    X .+= (1, 2) .+ 10
    @test X == Dataset([11 11 11; 12 12 12], :auto)

    ds = Dataset(rand(2, 3), :auto)
    @test floor.(Int, ds ./ (1,)) == Dataset(zeros(Int, 2, 3), :auto)
    ds .= floor.(Int, ds ./ (1,))
    @test ds == Dataset(zeros(2, 3), :auto)

    ds = Dataset(rand(2, 3), :auto)
    @test_throws InexactError convert.(Int, ds)
    ds2 = convert.(Int, floor.(ds))
    @test ds2 == Dataset(zeros(Int, 2, 3), :auto)
    @test eltype.(eachcol(ds2)) == [Union{Missing, Int}, Union{Missing, Int}, Union{Missing, Int}]
end

@testset "scalar on assignment side" begin
    ds = Dataset(rand(2, 3), :auto)
    @test_throws MethodError ds[1, 1] .= ds[1, 1] .- ds[1, 1]
    ds[1, 1:1] .= ds[1, 1] .- ds[1, 1]
    @test ds[1, 1] == 0
    @test_throws MethodError ds[1, 2] .-= ds[1, 2]
    ds[1:1, 2] .-= ds[1, 2]
    @test ds[1, 2] == 0
end

@testset "nothing test" begin
    X = Dataset(Any[1 2; 3 4], :auto)
    X .= nothing
    @test (X .== nothing) == Dataset(trues(2, 2), :auto)

    X = Dataset([1 2; 3 4], :auto)
    @test_throws MethodError X .= nothing
    @test X == Dataset([1 2; 3 4], :auto)

    X = Dataset([1 2; 3 4], :auto)
    foreach(i -> X[!, i] .= nothing, axes(X, 2))
    @test (X .== nothing) == Dataset(trues(2, 2), :auto)
end

######
@testset "aliasing test" begin
    ds = Dataset(x=[1, 2])
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)

    y = view(ds.x, [2, 1])
    ds .= y
    @test t1 != IMD._get_lastmodified(IMD._attributes(ds))

    @test ds[!, :x].val == [2, 1]

    ds = Dataset(x=[1, 2])
    y = view(ds.x, [2, 1])
    dsv = view(ds, :, :)
    dsv .= y
    @test ds.x == [2, 1]

    ds = Dataset(x=2, y=1, z=1)
    dsr = ds[1, :]
    y = view(ds.x.val, 1)
    dsr .= 2 .* y
    @test Vector(dsr) == [4, 4, 4]

    ds = Dataset(x=[1, 2], y=[11, 12])
    ds2 = Dataset()
    ds2.x = [-1, -2]
    ds2.y = ds.x
    ds3 = copy(ds2)
    ds .= ds2
    @test ds == ds3

    Random.seed!(1234)
    for i in 1:10
        ds1 = Dataset(rand(100, 100), :auto)
        ds2 = copy(ds1)
        for i in 1:100
            ds2[!, rand(1:100)] = ds1[!, i]
        end
        ds3 = copy(ds2)
        ds1 .= ds2
        @test ds1 == ds3
        @test ds2 != ds3
    end

    for i in 1:10
        ds1 = Dataset(rand(100, 100), :auto)
        ds2 = copy(ds1)
        for i in 1:100
            ds2[!, rand(1:100)] = ds1[!, i]
        end
        ds3 = copy(ds2)
        ds1 .= view(ds2, :, :)
        @test ds1 == ds3
        @test ds2 != ds3
    end

    for i in 1:10
        ds1 = Dataset(rand(100, 100), :auto)
        ds2 = copy(ds1)
        for i in 1:100
            ds2[!, rand(1:100)] = ds1[!, i]
        end
        ds3 = copy(ds2)
        view(ds1, :, :) .= ds2
        @test ds1 == ds3
        @test ds2 != ds3
    end

    for i in 1:10
        ds1 = Dataset(rand(100, 100), :auto)
        ds2 = copy(ds1)
        ds3 = copy(ds1)
        for i in 1:100
            ds2[!, rand(1:100)] = ds1[!, i]
            ds3[!, rand(1:100)] = ds1[!, i]
        end
        ds6 = copy(ds2)
        ds7 = copy(ds3)
        ds4 = Dataset(sin.(ds1[1, 1] .+ copy(ds1[!, 1]) .+ Matrix(ds2) ./ Matrix(ds3)), names(ds3))
        ds5 = sin.(view(ds1, 1, 1) .+ ds1[!, 1] .+ ds2 ./ ds3)
        ds1 .= sin.(view(ds1, 1, 1) .+ ds1[!, 1] .+ ds2 ./ ds3)
        @test ds1 == ds4 == ds5
        @test ds2 != ds6
        @test ds3 != ds7
    end

    for i in 1:10
        ds1 = Dataset(rand(100, 100), :auto)
        ds2 = copy(ds1)
        ds3 = copy(ds1)
        for i in 1:100
            ds2[!, rand(1:100)] = ds1[!, i]
            ds3[!, rand(1:100)] = ds1[!, i]
        end
        ds6 = copy(ds2)
        ds7 = copy(ds3)
        ds4 = Dataset(sin.(ds1[1, 1] .+ copy(ds1[!, 1]) .+ Matrix(ds2) ./ Matrix(ds3)), names(ds3))
        ds5 = sin.(view(ds1, 1, 1) .+ ds1[!, 1] .+ view(ds2, :, :) ./ ds3)
        ds1 .= sin.(view(ds1[!, 1], 1) .+ view(ds1[!, 1], :) .+ ds2 ./ view(ds3, :, :))
        @test ds1 == ds4 == ds5
        @test ds2 != ds6
        @test ds3 != ds7
    end

    for i in 1:10
        ds1 = Dataset(rand(100, 100), :auto)
        t1 = IMD._get_lastmodified(IMD._attributes(ds1))
        sleep(.1)

        ds2 = copy(ds1)
        ds3 = copy(ds1)
        for i in 1:100
            ds2[!, rand(1:100)] = ds1[!, i]
            ds3[!, rand(1:100)] = ds1[!, i]
        end
        ds6 = copy(ds2)
        ds7 = copy(ds3)
        ds4 = Dataset(sin.(ds1[1, 1] .+ copy(ds1[!, 1]) .+ Matrix(ds2) ./ Matrix(ds3)), names(ds3))
        ds5 = sin.(view(ds1, 1, 1) .+ ds1[!, 1] .+ view(ds2, :, :) ./ ds3)
        view(ds1, :, :) .= sin.(view(ds1[!, 1], 1) .+ view(ds1[!, 1], :) .+ ds2 ./ view(ds3, :, :))
        @test t1 != IMD._get_lastmodified(IMD._attributes(ds1))

        @test ds1 == ds4 == ds5
        @test ds2 != ds6
        @test ds3 != ds7
    end
end

@testset "@. test" begin
    ds = Dataset(rand(2, 3), :auto)
    sds = view(ds, 1:1, :)
    dsm = Matrix(ds)
    sdsm = Matrix(sds)

    r1 = @. (ds + sds + 5) / sds
    @test r1 isa Dataset

    @. ds = sin(sds / (ds + 1))
    @. dsm = sin(sdsm / (dsm + 1))
    @test ds == Dataset(dsm, names(ds))
end

@testset "test common cases" begin
    m = rand(1000, 10)
    ds = Dataset(m, :auto)
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)

    @test ds .+ 1 == Dataset(m .+ 1, names(ds))
    @test ds .+ transpose(1:10) == Dataset(m .+ transpose(1:10), names(ds))
    @test ds .+ (1:1000) == Dataset(m .+ (1:1000), names(ds))
    @test ds .+ m == Dataset(m .+ m, names(ds))
    @test m .+ ds == Dataset(m .+ m, names(ds))
    @test ds .+ ds == Dataset(m .+ m, names(ds))

    ds .+= 1
    @test t1 != IMD._get_lastmodified(IMD._attributes(ds))

    m .+= 1
    @test ds == Dataset(m, names(ds))
    ds .+= transpose(1:10)
    m .+= transpose(1:10)
    @test ds == Dataset(m, names(ds))
    ds .+= (1:1000)
    m .+= (1:1000)
    @test ds == Dataset(m, names(ds))
    ds .+= ds
    m .+= m
    @test ds == Dataset(m, names(ds))
    ds2 = copy(ds)
    m2 = copy(m)
    ds .+= ds .+ ds2 .+ m2 .+ 1
    m .+= m .+ ds2 .+ m2 .+ 1
    @test ds == Dataset(m, names(ds))
end

@testset "data set only on left hand side broadcasting assignment" begin
    Random.seed!(1234)

    m = rand(3, 4);
    m = allowmissing(m)
    m2 = copy(m);
    m3 = copy(m);
    ds = Dataset(a=view(m, :, 1), b=view(m, :, 1),
                   c=view(m, :, 1), d=view(m, :, 1), copycols=false);
    ds2 = copy(ds)
    mds = Matrix(ds)

    @test m .+ ds == m2 .+ ds
    @test Matrix(m .+ ds) == m .+ mds
    @test sin.(m .+ ds) .+ 1 .+ m2 == sin.(m2 .+ ds) .+ 1 .+ m
    @test Matrix(m .+ ds ./ 2 .* ds2) == m .+ mds ./ 2 .* mds

    m2 .+= ds .+ 1 ./ ds2
    m .+= ds .+ 1 ./ ds2
    @test m2 == m
    for col in eachcol(ds)
        @test col.val == m[:, 1]
    end
    for col in eachcol(ds2)
        @test col.val == m3[:, 1]
    end

    m = rand(3, 4);
    m = allowmissing(m)
    m2 = copy(m);
    m3 = copy(m);
    ds = view(Dataset(a=view(m, :, 1), b=view(m, :, 1),
                        c=view(m, :, 1), d=view(m, :, 1), copycols=false),
              [3, 2, 1], :)
    ds2 = copy(ds)
    mds = Matrix(ds)

    @test m .+ ds == m2 .+ ds
    @test Matrix(m .+ ds) == m .+ mds
    @test sin.(m .+ ds) .+ 1 .+ m2 == sin.(m2 .+ ds) .+ 1 .+ m
    @test Matrix(m .+ ds ./ 2 .* ds2) == m .+ mds ./ 2 .* mds

    m2 .+= ds .+ 1 ./ ds2
    m .+= ds .+ 1 ./ ds2
    @test m2 == m
    for col in eachcol(ds)
        @test IMD.__!(col) == m[3:-1:1, 1]
    end
    for col in eachcol(ds2)
        @test IMD.__!(col) == m3[3:-1:1, 1]
    end
end

@testset "broadcasting with 3-dimensional object" begin
    y = zeros(4, 3, 2)
    ds = Dataset(ones(4, 3), :auto)
    @test_throws DimensionMismatch ds .+ y
    @test_throws DimensionMismatch y .+ ds
    @test_throws DimensionMismatch ds .+= y
    y .+= ds
    @test y == ones(4, 3, 2)
end

@testset "additional checks of post-! broadcasting rules" begin
    ds = copy(refds)
    v1 = ds[!, 1]
    @test_throws MethodError ds[CartesianIndex(1, 1)] .= 1
    @test_throws MethodError ds[CartesianIndex(1, 1)] .= "d"
    @test_throws DimensionMismatch ds[CartesianIndex(1, 1)] .= [1, 2]

    ds = copy(refds)
    v1 = ds[!, 1]
    @test_throws MethodError ds[1, 1] .= 1
    @test_throws MethodError ds[1, 1] .= "d"
    @test_throws DimensionMismatch ds[1, 1] .= [1, 2]

    ds = copy(refds)
    v1 = ds[!, 1]
    @test_throws MethodError ds[1, :x1] .= 1
    @test_throws MethodError ds[1, :x1] .= "d"
    @test_throws DimensionMismatch ds[1, :x1] .= [1, 2]

    ds = copy(refds)
    v1 = ds[!, 1]
    @test_throws MethodError ds[1, "x1"] .= 1
    @test_throws MethodError ds[1, "x1"] .= "d"
    @test_throws DimensionMismatch ds[1, "x1"] .= [1, 2]

    ds = copy(refds)
    v1 = ds[!, 1]
    v2 = ds[!, 2]
    ds[1, 1:2] .= 'd'
    @test v1 == [100.0, 2.5, 3.5]
    @test v2 == [100.0, 5.5, 6.5]
    @test_throws MethodError ds[1, 1:2] .= "d"
    @test v1 == [100.0, 2.5, 3.5]
    @test v2 == [100.0, 5.5, 6.5]
    ds[1, 1:2] .= 'e':'f'
    @test v1 == [101.0, 2.5, 3.5]
    @test v2 == [102.0, 5.5, 6.5]
    @test_throws DimensionMismatch ds[1, 1:2] .= ['d' 'd']
    @test v1 == [101.0, 2.5, 3.5]
    @test v2 == [102.0, 5.5, 6.5]

    ds = copy(refds)
    v1 = ds[!, 1]
    ds[:, 1] .= 'd'
    @test v1 == [100.0, 100.0, 100.0]
    @test_throws MethodError ds[:, 1] .= "d"
    @test v1 == [100.0, 100.0, 100.0]
    @test_throws DimensionMismatch ds[:, 1] .= [1 2 3]
    @test v1 == [100.0, 100.0, 100.0]

    ds = copy(refds)
    v1 = ds[!, 1]
    ds[:, :x1] .= 'd'
    @test v1 == [100.0, 100.0, 100.0]
    @test_throws MethodError ds[:, :x1] .= "d"
    @test v1 == [100.0, 100.0, 100.0]
    @test_throws DimensionMismatch ds[:, :x1] .= [1 2 3]
    @test v1 == [100.0, 100.0, 100.0]

    ds = copy(refds)
    v1 = ds[!, 1]
    ds[:, 1] .= 'd':'f'
    @test v1 == [100.0, 101.0, 102.0]
    @test_throws MethodError ds[:, 1] .= ["d", "e", "f"]
    @test v1 == [100.0, 101.0, 102.0]

    ds = copy(refds)
    t1 = IMD._get_lastmodified(IMD._attributes(ds))
    sleep(.1)

    v1 = ds[!, 1]
    v2 = ds[!, 2]
    ds[:, 1:2] .= 'd'
    @test t1 != IMD._get_lastmodified(IMD._attributes(ds))

    @test v1 == [100.0, 100.0, 100.0]
    @test v2 == [100.0, 100.0, 100.0]
    @test_throws MethodError ds[:, 1:2] .= "d"
    @test v1 == [100.0, 100.0, 100.0]
    @test v2 == [100.0, 100.0, 100.0]
    @test_throws DimensionMismatch ds[:, 1:2] .= [1 2 3]
    @test v1 == [100.0, 100.0, 100.0]
    @test v2 == [100.0, 100.0, 100.0]

    ds = copy(refds)
    v1 = ds[!, 1]
    v2 = ds[!, 2]
    ds[:, 1:2] .= 'd':'f'
    @test v1 == [100.0, 101.0, 102.0]
    @test v2 == [100.0, 101.0, 102.0]
    @test_throws MethodError ds[:, 1:2] .= ["d", "e", "f"]
    @test v1 == [100.0, 101.0, 102.0]
    @test v2 == [100.0, 101.0, 102.0]

    ds = copy(refds)
    v1 = ds[!, 1]
    v2 = ds[!, 2]
    ds[:, 1:2] .= permutedims('d':'e')
    @test v1 == [100.0, 100.0, 100.0]
    @test v2 == [101.0, 101.0, 101.0]

    ds = copy(refds)
    v1 = ds[!, 1]
    v2 = ds[!, 2]
    ds[:, 1:2] .= reshape('d':'i', 3, :)
    @test v1 == [100.0, 101.0, 102.0]
    @test v2 == [103.0, 104.0, 105.0]
    @test_throws DimensionMismatch ds[:, 1:2] .= reshape('d':'i', 1, :, 3)
    @test v1 == [100.0, 101.0, 102.0]
    @test v2 == [103.0, 104.0, 105.0]

    ds = copy(refds)
    v1 = ds[!, 1]
    v1′ = ds[:, 1]
    ds[!, 1] .= 100.0
    @test ds.x1 == [100.0, 100.0, 100.0]
    @test v1 == v1′
    ds[!, 1] .= 'd'
    @test ds.x1 == ['d', 'd', 'd']
    @test v1 == v1′
    @test_throws DimensionMismatch ds[!, 1] .= [1 2 3]
    @test ds.x1 == ['d', 'd', 'd']
    @test v1 == v1′

    ds = copy(refds)
    v1 = ds[!, 1]
    v1′ = ds[:, 1]
    ds[!, :x1] .= 100.0
    @test ds.x1 == [100.0, 100.0, 100.0]
    @test v1 == v1′
    ds[!, :x1] .= 'd'
    @test ds.x1 == ['d', 'd', 'd']
    @test v1 == v1′
    @test_throws DimensionMismatch ds[!, :x1] .= [1 2 3]
    @test ds.x1 == ['d', 'd', 'd']
    @test v1 == v1′

    ds = copy(refds)
    ds[!, :newcol] .= 100.0
    @test ds.newcol == [100.0, 100.0, 100.0]
    @test ds[:, 1:end-1] == refds

    ds = copy(refds)
    ds[!, "newcol"] .= 100.0
    @test ds.newcol == [100.0, 100.0, 100.0]
    @test ds[:, 1:end-1] == refds

    ds = copy(refds)
    ds[!, :newcol] .= 'd'
    @test ds.newcol == ['d', 'd', 'd']
    @test ds[:, 1:end-1] == refds

    ds = copy(refds)
    ds[!, "newcol"] .= 'd'
    @test ds.newcol == ['d', 'd', 'd']
    @test ds[:, 1:end-1] == refds

    ds = copy(refds)
    @test_throws DimensionMismatch ds[!, :newcol] .= [1 2 3]
    @test ds == refds

    ds = copy(refds)
    @test_throws DimensionMismatch ds[!, "newcol"] .= [1 2 3]
    @test ds == refds

    ds = copy(refds)
    @test_throws ArgumentError ds[!, 10] .= 'a'
    @test ds == refds
    @test_throws ArgumentError ds[!, 10] .= [1, 2, 3]
    @test ds == refds
    @test_throws ArgumentError ds[!, 10] .= [1 2 3]
    @test ds == refds

    ds = copy(refds)
    ds[!, 1:2] .= 'a'
    @test Matrix(ds) == ['a'  'a'  7.5  10.5  13.5
                         'a'  'a'  8.5  11.5  14.5
                         'a'  'a'  9.5  12.5  15.5]

    ds = copy(refds)
    v1 = ds[!, 1]
    ds[:, :x1] .= 'd'
    @test v1.val == [100.0, 100.0, 100.0]
    @test_throws MethodError ds[:, 1] .= "d"
    @test v1.val == [100.0, 100.0, 100.0]
    @test_throws DimensionMismatch ds[:, 1] .= [1 2 3]
    @test v1.val == [100.0, 100.0, 100.0]

    ds = copy(refds)
    if isdefined(Base, :dotgetproperty)
        ds.newcol .= 'd'
        @test ds == [refds Dataset(newcol=fill('d', 3))]
    else
        @test_throws ArgumentError ds.newcol .= 'd'
        @test ds == refds
    end

    ds = view(copy(refds), :, :)
    v1 = ds[!, 1]
    @test_throws MethodError ds[CartesianIndex(1, 1)] .= 1
    @test_throws MethodError ds[CartesianIndex(1, 1)] .= "d"
    @test_throws DimensionMismatch ds[CartesianIndex(1, 1)] .= [1, 2]

    ds = view(copy(refds), :, :)
    v1 = ds[!, 1]
    @test_throws MethodError ds[1, 1] .= 1
    @test_throws MethodError ds[1, 1] .= "d"
    @test_throws DimensionMismatch ds[1, 1] .= [1, 2]

    ds = view(copy(refds), :, :)
    v1 = ds[!, 1]
    @test_throws MethodError ds[1, :x1] .= 1
    @test_throws MethodError ds[1, :x1] .= "d"
    @test_throws DimensionMismatch ds[1, :x1] .= [1, 2]

    ds = view(copy(refds), :, :)
    v1 = ds[!, 1]
    v2 = ds[!, 2]
    ds[1, 1:2] .= 'd'
    @test v1 == [100.0, 2.5, 3.5]
    @test v2 == [100.0, 5.5, 6.5]
    @test_throws MethodError ds[1, 1:2] .= "d"
    @test v1 == [100.0, 2.5, 3.5]
    @test v2 == [100.0, 5.5, 6.5]
    ds[1, 1:2] .= 'e':'f'
    @test v1 == [101.0, 2.5, 3.5]
    @test v2 == [102.0, 5.5, 6.5]
    @test_throws DimensionMismatch ds[1, 1:2] .= ['d' 'd']
    @test v1 == [101.0, 2.5, 3.5]
    @test v2 == [102.0, 5.5, 6.5]

    ds = view(copy(refds), :, :)
    v1 = ds[!, 1]
    ds[:, 1] .= 'd'
    @test v1 == [100.0, 100.0, 100.0]
    @test_throws MethodError ds[:, 1] .= "d"
    @test v1 == [100.0, 100.0, 100.0]
    @test_throws DimensionMismatch ds[:, 1] .= [1 2 3]
    @test v1 == [100.0, 100.0, 100.0]

    ds = view(copy(refds), :, :)
    v1 = ds[!, 1]
    ds[:, :x1] .= 'd'
    @test v1 == [100.0, 100.0, 100.0]
    @test_throws MethodError ds[:, :x1] .= "d"
    @test v1 == [100.0, 100.0, 100.0]
    @test_throws DimensionMismatch ds[:, :x1] .= [1 2 3]
    @test v1 == [100.0, 100.0, 100.0]

    ds = view(copy(refds), :, :)
    v1 = ds[!, 1]
    ds[:, 1] .= 'd':'f'
    @test v1 == [100.0, 101.0, 102.0]
    @test_throws MethodError ds[:, 1] .= ["d", "e", "f"]
    @test v1 == [100.0, 101.0, 102.0]

    ds = view(copy(refds), :, :)
    v1 = ds[!, 1]
    v2 = ds[!, 2]
    ds[:, 1:2] .= 'd'
    @test v1 == [100.0, 100.0, 100.0]
    @test v2 == [100.0, 100.0, 100.0]
    @test_throws MethodError ds[:, 1:2] .= "d"
    @test v1 == [100.0, 100.0, 100.0]
    @test v2 == [100.0, 100.0, 100.0]
    @test_throws DimensionMismatch ds[:, 1:2] .= [1 2 3]
    @test v1 == [100.0, 100.0, 100.0]
    @test v2 == [100.0, 100.0, 100.0]

    ds = view(copy(refds), :, :)
    v1 = ds[!, 1]
    v2 = ds[!, 2]
    ds[:, 1:2] .= 'd':'f'
    @test v1 == [100.0, 101.0, 102.0]
    @test v2 == [100.0, 101.0, 102.0]
    @test_throws MethodError ds[:, 1:2] .= ["d", "e", "f"]
    @test v1 == [100.0, 101.0, 102.0]
    @test v2 == [100.0, 101.0, 102.0]

    ds = view(copy(refds), :, :)
    v1 = ds[!, 1]
    v2 = ds[!, 2]
    ds[:, 1:2] .= permutedims('d':'e')
    @test v1 == [100.0, 100.0, 100.0]
    @test v2 == [101.0, 101.0, 101.0]

    ds = view(copy(refds), :, :)
    v1 = ds[!, 1]
    v2 = ds[!, 2]
    ds[:, 1:2] .= reshape('d':'i', 3, :)
    @test v1 == [100.0, 101.0, 102.0]
    @test v2 == [103.0, 104.0, 105.0]
    @test_throws DimensionMismatch ds[:, 1:2] .= reshape('d':'i', 1, :, 3)
    @test v1 == [100.0, 101.0, 102.0]
    @test v2 == [103.0, 104.0, 105.0]

    ds = view(copy(refds), :, :)
    ds[!, 1] .= 100
    @test parent(ds).x1 == [100, 100, 100]
    @test eltype(parent(ds).x1) == Union{Missing, Float64}

    ds = view(copy(refds), :, :)
    ds[!, :x1] .= 100.0
    @test parent(ds).x1 == [100, 100, 100]
    @test eltype(parent(ds).x1) == Union{Missing, Float64}

    # we don't support adding :newcol to views
    # ds = view(copy(refds), :, :)
    # ds[!, :newcol] .= 100.0
    # @test parent(ds).newcol == [100, 100, 100]
    # @test eltype(parent(ds).newcol) == Union{Float64, Missing}

    ds = view(copy(refds), :, :)
    @test_throws BoundsError ds[!, 10] .= 'a'
    @test ds == refds
    @test_throws BoundsError ds[!, 10] .= [1, 2, 3]
    @test ds == refds
    @test_throws BoundsError ds[!, 10] .= [1 2 3]
    @test ds == refds

    ds = view(copy(refds), :, :)
    ds[!, 1:2] .= 'a'
    @test parent(ds).x1.val == parent(ds).x2.val == [97.0,97,97]

    ds = view(copy(refds), :, :)
    v1 = ds[!, 1]
    ds.x1 .= 'd'
    @test v1 == [100.0, 100.0, 100.0]
    @test_throws MethodError ds[:, 1] .= "d"
    @test v1 == [100.0, 100.0, 100.0]
    @test_throws DimensionMismatch ds[:, 1] .= [1 2 3]
    @test v1 == [100.0, 100.0, 100.0]

    # ds = view(copy(refds), :, :)
    # if VERSION >= v"1.7"
    #     ds.newcol .= 'd'
    #     @test ds.newcol == fill('d', 3)
    # else
    #     @test_throws ArgumentError ds.newcol .= 'd'
    #     @test ds == refds
    # end
end

@testset "DatasetRow getproperty broadcasted assignment" begin
    ds = Dataset(a=[[1, 2], [3, 4]], b=[[5, 6], [7, 8]])
    dsr = ds[1, :]
    dsr.a .= 10
    @test ds == Dataset(a=[[10, 10], [3, 4]], b=[[5, 6], [7, 8]])
    @test_throws MethodError dsr.a .= ["a", "b"]

    ds = Dataset(a=[[1, 2], [3, 4]], b=[[5, 6], [7, 8]])
    dsr = ds[1, 1:1]
    dsr.a .= 10
    @test ds == Dataset(a=[[10, 10], [3, 4]], b=[[5, 6], [7, 8]])
    @test_throws MethodError dsr.a .= ["a", "b"]

    ds = Dataset(a=[[1, 2], [3, 4]], b=[[5, 6], [7, 8]])
    dsr = ds[1, :]
    dsr."a" .= 10
    @test ds == Dataset(a=[[10, 10], [3, 4]], b=[[5, 6], [7, 8]])
    @test_throws MethodError dsr."a" .= ["a", "b"]

    ds = Dataset(a=[[1, 2], [3, 4]], b=[[5, 6], [7, 8]])
    dsr = ds[1, 1:1]
    dsr."a" .= 10
    @test ds == Dataset(a=[[10, 10], [3, 4]], b=[[5, 6], [7, 8]])
    @test_throws MethodError dsr."a" .= ["a", "b"]
end

@testset "make sure that : is in place and ! allocates" begin
    ds = Dataset(a=[1, 2, 3])
    a = ds.a
    ds[:, :a] .+= 1
    @test a == [2, 3, 4]
    @test ds.a === a
    ds[!, :a] .+= 1
    @test a == [2, 3, 4]
    @test ds.a == [3, 4, 5]
    @test ds.a !== a

    ds = Dataset(a=[1, 2, 3])
    a = ds.a
    ds[:, "a"] .+= 1
    @test a == [2, 3, 4]
    @test ds.a === a
    ds[!, "a"] .+= 1
    @test a == [2, 3, 4]
    @test ds.a == [3, 4, 5]
    @test ds.a !== a
end

@testset "add new correct rules for ds[row, col] .= v broadcasting" begin
    for v in [:a, "a"]
        ds = Dataset(a=1)
        @test_throws MethodError ds[1, 1] .= 10
        @test_throws MethodError ds[1, v] .= 10
        @test_throws MethodError ds[CartesianIndex(1, 1)] .= 10
        ds = Dataset(a=[[1, 2, 3]])
        ds[1, 1] .= 10
        @test ds == Dataset(a=[[10, 10, 10]])
        ds[1, v] .= 100
        @test ds == Dataset(a=[[100, 100, 100]])
        ds[CartesianIndex(1, 1)] .= 1000
        @test ds == Dataset(a=[[1000, 1000, 1000]])
    end
end

@testset "broadcasting into ds[!, cols]" begin
    for selector in [1:2, Between(:x1, :x2), Not(r"x3"), [:x1, :x2],
                     ["x1", "x2"], Between("x1", "x2")]
        ds = Dataset(x1=1:3, x2=4:6)
        ds[!, selector] .= "a"
        @test ds == Dataset(fill("a", 3, 2), :auto)
        @test ds.x1 !== ds.x2

        ds = Dataset(x1=1:3, x2=4:6)
        ds[!, selector] .= Ref((a=1, b=2))
        @test ds == Dataset(fill((a=1, b=2), 3, 2), :auto)
        @test ds.x1 !== ds.x2

        ds = Dataset(x1=1:3, x2=4:6)
        ds[!, selector] .= ["a" "b"]
        @test ds == Dataset(["a" "b"
                               "a" "b"
                               "a" "b"], :auto)
        @test ds.x1 !== ds.x2

        ds = Dataset(x1=1:3, x2=4:6)
        ds[!, selector] .= ["a", "b", "c"]
        @test ds == Dataset(["a" "a"
                               "b" "b"
                               "c" "c"], :auto)
        @test ds.x1 !== ds.x2

        ds = Dataset(x1=1:3, x2=4:6)
        ds[!, selector] .= categorical(["a"])
        @test ds == Dataset(["a" "a"
                               "a" "a"
                               "a" "a"], :auto)
        @test ds.x1.val isa CategoricalVector
        @test ds.x2.val isa CategoricalVector
        @test ds.x1 !== ds.x2

        ds = Dataset(x1=1:3, x2=4:6)
        ds[!, selector] .= Dataset(["a" "b"], :auto)
        @test ds == Dataset(["a" "b"
                               "a" "b"
                               "a" "b"], :auto)
        @test ds.x1 !== ds.x2

        ds = Dataset(x1=1:3, x2=4:6)
        ds[!, selector] .= Dataset(["a" "d"
                                      "b" "e"
                                      "c" "f"], :auto)
        @test ds == Dataset(["a" "d"
                               "b" "e"
                               "c" "f"], :auto)
        @test ds.x1 !== ds.x2

        ds = Dataset(x1=1:3, x2=4:6)
        ds[!, selector] .= ["a" "d"
                            "b" "e"
                            "c" "f"]
        @test ds == Dataset(["a" "d"
                               "b" "e"
                               "c" "f"], :auto)
        @test ds.x1 !== ds.x2

        ds = Dataset(x1=1:3, x2=4:6, x3=1)
        ds[!, selector] .= "a"
        @test ds == Dataset(["a" "a" 1
                               "a" "a" 1
                               "a" "a" 1], :auto)
        @test ds.x1 !== ds.x2

        ds = Dataset(x1=1:3, x2=4:6, x3=1)
        ds[!, selector] .= Ref((a=1, b=2))
        @test ds[:, 1:2] == Dataset(fill((a=1, b=2), 3, 2), :auto)
        @test ds[:, 3] == [1, 1, 1]
        @test ds.x1 !== ds.x2

        ds = Dataset(x1=1:3, x2=4:6, x3=1)
        ds[!, selector] .= ["a" "b"]
        @test ds == Dataset(["a" "b" 1
                               "a" "b" 1
                               "a" "b" 1], :auto)
        @test ds.x1 !== ds.x2

        ds = Dataset(x1=1:3, x2=4:6, x3=1)
        ds[!, selector] .= ["a", "b", "c"]
        @test ds == Dataset(["a" "a" 1
                               "b" "b" 1
                               "c" "c" 1], :auto)
        @test ds.x1 !== ds.x2

        ds = Dataset(x1=1:3, x2=4:6, x3=1)
        ds[!, selector] .= categorical(["a"])
        @test ds == Dataset(["a" "a" 1
                               "a" "a" 1
                               "a" "a" 1], :auto)
        @test ds.x1.val isa CategoricalVector
        @test ds.x2.val isa CategoricalVector
        @test ds.x1 !== ds.x2

        ds = Dataset(x1=1:3, x2=4:6, x3=1)
        ds[!, selector] .= Dataset(["a" "b"], :auto)
        @test ds == Dataset(["a" "b" 1
                               "a" "b" 1
                               "a" "b" 1], :auto)
        @test ds.x1 !== ds.x2

        ds = Dataset(x1=1:3, x2=4:6, x3=1)
        ds[!, selector] .= Dataset(["a" "d"
                                      "b" "e"
                                      "c" "f"], :auto)
        @test ds == Dataset(["a" "d" 1
                               "b" "e" 1
                               "c" "f" 1], :auto)
        @test ds.x1 !== ds.x2

        ds = Dataset(x1=1:3, x2=4:6, x3=1)
        ds[!, selector] .= ["a" "d"
                            "b" "e"
                            "c" "f"]
        @test ds == Dataset(["a" "d" 1
                               "b" "e" 1
                               "c" "f" 1], :auto)
        @test ds.x1 !== ds.x2
    end

    ds = Dataset(x1=1:3, x2=4:6)
    @test_throws ArgumentError ds[!, [:x1, :x3]] .= "a"
end

@testset "broadcasting over heterogenous columns" begin
    ds = Dataset(x=[1, 1.0, big(1), "1"])
    f_identity(x) = x
    @test ds == f_identity.(ds)
end

@testset "@views on ds[!, col]" begin
    ds = Dataset(ones(3, 4), :auto)
    @views ds[!, 1] .+= 1
    @test ds[!, 1] == [2.0, 2.0, 2.0]
    @views ds[:, 2] .= ds[!, 4] .+ ds[!, 3]
    @test ds[!, 2] == [2.0, 2.0, 2.0]

    # make sure we do not mess with maybeview
    @test @views typeof(ds[!, 1:2]) <: SubDataset
end

@testset "broadcasting of ds[:, col] = value" begin
    ds = Dataset(ones(3, 4), :auto)
    z = ["a", "b", "c"]
    ds[:, :z] .= z
    @test ds.z == z
    @test ds.z !== z
    @test_throws ArgumentError ds[:, 6] .= z
    @test_throws MethodError ds[:, 1] .= z

    ds = Dataset(ones(3, 4), :auto)
    z = "abc"
    ds[:, :z] .= z
    @test ds.z == fill("abc", 3)
    @test_throws ArgumentError ds[:, 6] .= z
    @test_throws MethodError ds[:, 1] .= z

    ds = Dataset(ones(3, 4), :auto)
    z = fill("abc", 1, 1, 2)
    @test_throws DimensionMismatch ds[:, :z] .= z

    ds = Dataset(ones(3, 4), :auto)
    z = ["a", "b", "c"]
    ds[:, "z"] .= z
    @test ds.z == z
    @test ds.z !== z

    ds = Dataset(ones(3, 4), :auto)
    z = "abc"
    ds[:, "z"] .= z
    @test ds.z == fill("abc", 3)

    ds = Dataset(ones(3, 4), :auto)
    z = fill("abc", 1, 1, 2)
    @test_throws DimensionMismatch ds[:, "z"] .= z
end

@testset "broadcasting of getproperty" begin
    if isdefined(Base, :dotgetproperty)
        ds = Dataset(a=1:4)
        ds.b .= 1
        ds.c .= 4:-1:1
        # TODO: enable this in the future when the deprecation period is finished
        # ds.a .= 'a':'d'
        # @test ds.a isa Vector{Char}
        # @test ds == Dataset(a='a':'d', b=1, c=4:-1:1)
        # dsv = view(ds, 2:3, 2:3)
        # @test_throws ArgumentError dsv.b .= 0
    end
end


text:
# common utils for applications -- staging area towards RoME itself

using Dates
using DelimitedFiles
using RoME
using DSP
using ArgParse

# NOTE, moved to RoME
# @enum HandlerStateMachine HSMReady HSMHandling HSMOverlapHandling HSMBlocking
# @enum SolverStateMachine SSMReady SSMConsumingSolvables SSMSolving

# Accumulate delta X values using FGOS
# TODO, refactor into RoME
function devAccumulateOdoPose2(DX::Array{Float64,2},
                               X0::Vector{Float64}=zeros(3);
                               P0=1e-3*Matrix(LinearAlgebra.I, 3,3),
                               Qc=1e-6*Matrix(LinearAlgebra.I, 3,3),
                               dt::Float64=1.0,
                               timestamp::DateTime=now() )
  #
  # entries are rows with columns dx,dy,dtheta
  @assert size(DX,2) == 3
  mpp = MutablePose2Pose2Gaussian(MvNormal(X0, P0), timestamp=timestamp )
  nXYT = zeros(size(DX,1), 3)
  for i in 1:size(DX,1)
    RoME.accumulateDiscreteLocalFrame!(mpp,DX[i, :],Qc,dt)
    nXYT[i,:] .= mpp.Zij.μ
  end

  return nXYT
end

function loadResultsDRT(fg)
  drt_data = readdlm(joinLogPath(fg, "DRT.csv"), ',')

  DRTT = DateTime.(drt_data[:,1])
  XX = Float64.(drt_data[:,4])
  YY = Float64.(drt_data[:,5])

  # filter the signals for hard jumps between drt transitions
  responsetype = Lowpass(1.0; fs=50)
  designmethod = FIRWindow(hanning(64))

  # remove those near zero
  mask = 40.0 .< (XX.^2 + YY.^2)

  TTm = DRTT[mask]
  XXm = XX[mask]
  YYm = YY[mask]
  XXf = filt(digitalfilter(responsetype, designmethod), XXm)
  YYf = filt(digitalfilter(responsetype, designmethod), YYm)

  return drt_data, TTm, XXm, YYm, XXf, YYf
end

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--plotSeriesBeliefs"
            help = "Glob fg_* archives and draw belief frames as many as is requested, default 0 is for all"
            arg_type = Int
            default = -1
        "--plotSeriesBeliefsNew"
            help = "Glob fg_* archives and draw belief frames as many as is requested, default 0 is for all with background"
            arg_type = Int
            default = -1
        "--reportPoses"
            help = "Generate report on interpose pose factors"
            action = :store_true
        "--reportRanges"
            help = "Generate report on range factors"
            action = :store_true
        "--drawAllRanges"
            help = "Draw ranges as separate images"
            action = :store_true
        "--skip"
            help = "Skip existing images"
            action = :store_true
        "--drawFG"
            help = "Draw factor graph to PDF"
            action = :store_true
        "--reportDir"
            help = "which folder in which to produce results."
            required = false
        "--kappa_odo"
            help = "Scale the odometry covariance"
            arg_type = Float64
            default = 1.0
        "--stride_range"
            help = "Every how many poses to try add a range measurement"
            arg_type = Int
            default = 4
        "--fixedlag"
            help = "FIFO fixed-lag solve length"
            arg_type = Int
            default = 30
        "--magStdDeg"
            help = "Magnetometer standard deviation"
            arg_type = Float64
            default = 5.0
        "--iters", "-i"
            help = "LCM messages to handle"
            arg_type = Int
            default = 10000
        "--speed", "-s"
            help = "Target playback speed for LCMLog"
            arg_type = Float64
            default = 0.2
        "--dbg"
            help = "debug flag"
            action = :store_true
        "--genResults"
            help = "Generate output results"
            action = :store_true
        "--warmupjit"
            help = "Warm up JIT during startup"
            action = :store_true
        "--savePlotting"
            help = "Store factor graph after each pose"
            action = :store_true
        "--odoGyroBias"
            help = "Use gyro biased channel for odometry"
            action = :store_true
        "--limitfixeddown"
            help = "Limit numerical computations for recycled and marginalized cliques during down solve."
            action = :store_true
        "--recordTrees"
            help = "Store copies of the Bayes tree for later animation, sets the rate in sleep(1/rate)"
            arg_type = Float64
            default = -1.0
        "--scaleOdoX"
            help = "Scale the delta X odo parameters."
            arg_type = Float64
            default = 1.0
    end

    return parse_args(s)
end

#

ext: jl


text:
abstract type AbstractStateSpaceModel{T<:Real} end

struct StateSpaceModel{T} <: AbstractStateSpaceModel{T}

    # Process transition matrix, control matrix, and noise covariance
    A::Function
    B::Function
    V::Function

    # Observation matrix, feed-forward matrix, and noise covariance
    C::Function
    D::Function
    W::Function

    # Inital state and error covariance conditions
    x1::Vector{T}
    P1::Union{Matrix{T}, SparseMatrixCSC{T}}

    # Model dimensions
    nx::Int
    ny::Int
    nu::Int

    function StateSpaceModel{T}(A::Function, B::Function, V::Function,
                                C::Function, D::Function, W::Function,
                                x1::Vector{T}, P1::Union{Matrix{T}, SparseMatrixCSC{T}}) where {T}

        ispossemidef(x::AbstractMatrix) = issymmetric(x) && (eigmin(full(x)) >= 0)
        @assert ispossemidef(V(1))
        @assert ispossemidef(W(1))
        @assert ispossemidef(P1)

        nx, ny, nu = confirm_matrix_sizes(A(1), B(1), V(1), C(1), D(1), W(1), x1, P1)
        new(A, B, V, C, D, W, x1, P1, nx, ny, nu)
    end
end

# Time-dependent constructor
StateSpaceModel(A::Function, V::Function, C::Function, W::Function,
                       x1::Vector{T}, P1::Union{Matrix{T}, SparseMatrixCSC{T}};
                       B::Function=_->zeros(size(V(1), 1), 1),
                       D::Function=_->zeros(size(W(1), 1), 1)) where {T} =
	  StateSpaceModel{T}(A, B, V, C, D, W, x1, P1)

# Time-independent constructor
StateSpaceModel(A::Matrix{T}, V::Matrix{T}, C::Matrix{T}, W::Matrix{T},
                       x1::Vector{T}, P1::Union{Matrix{T}, SparseMatrixCSC{T}};
                       B::Matrix{T}=zeros(size(A, 1), 1),
                       D::Matrix{T}=zeros(size(C, 1), 1)) where {T} =
	  StateSpaceModel{T}(_->A, _->B, _->V, _->C, _->D, _->W, x1, P1)

function show(io::IO, mod::StateSpaceModel{T}) where T
    dx, dy = mod.nx, mod.ny
    println("StateSpaceModel{$T}, $dx-D process x $dy-D observations")
    println("Process evolution matrix A:")
    show(mod.A(1))
    println("\n\nControl input matrix B:")
    show(mod.B(1))
    println("\n\nProcess error covariance V:")
    show(mod.V(1))
    println("\n\nObservation matrix C:")
    show(mod.C(1))
    println("\n\nFeed-forward matrix D:")
    show(mod.D(1))
    println("\n\nObseration error covariance W:")
    show(mod.W(1))
end

# Time-independent parametrized matrix type
# TODO: A more general ParametrizedArray (allowing vectors) could be preferable
struct ParametrizedMatrix{T}

    f::AbstractVector{T}
    D::Union{Matrix{T}, SparseMatrixCSC{T}}
    np::Int
    dims::Tuple{Int, Int}

    function ParametrizedMatrix{T}(f::AbstractVector{T},
                D::Union{Matrix{T}, SparseMatrixCSC{T, Int}}, dims::Tuple{Int, Int}) where {T}
        @assert length(f) == size(D, 1)
        @assert length(f) == dims[1] * dims[2]
        new(f, D, size(D, 2), dims)
    end #ParametrizedMatrix

end #ParametrizedMatrix

ParametrizedMatrix(f::Vector{T}, D::Matrix{T}, dims::Tuple{Int, Int}) where {T} = ParametrizedMatrix{T}(f, D, dims)

function ParametrizedMatrix(f::SparseVector{T}, D::SparseMatrixCSC{T}, dims::Tuple{Int, Int}) where T
    ParametrizedMatrix{T}(f, D, dims)
end #ParametrizedMatrix

function show(io::IO, cpm::ParametrizedMatrix{T}) where T

    function combinestrelems(a, b)
        if (a != "") & (b != "")
          a * " + " * b
        elseif (a == "") & (b == "")
          "0"
        else
          a * b
        end #if
    end #combinestrelems

    nx, ny = cpm.dims
    println(nx, "x", ny, " ParametrizedMatrix{$T}")
    conststring = map(x -> x == 0 ? "" : string(x), cpm.f)
    paramstring = fill("", length(cpm.f))
    paramelems = fill("", cpm.np)

    for i = 1:size(cpm.D,1)
        for j = 1:cpm.np
            if  cpm.D[i,j] == 0
                paramelems[j] = ""
            elseif cpm.D[i,j] == 1
                paramelems[j] = string(Char(96+j))
            else
                paramelems[j] = string(cpm.D[i,j],Char(96+j))
            end #if
        end #for
        paramstring[i] = isempty(paramelems) ? "0" : reduce(combinestrelems, paramelems)
    end #for

    finalstrings = map(combinestrelems, conststring, paramstring)

    reshape(finalstrings, cpm.dims) |> showcompact
end

Base.length(pm::ParametrizedMatrix) = pm.dims[1] * pm.dims[2]
Base.size(pm::ParametrizedMatrix) = pm.dims
Base.size(pm::ParametrizedMatrix, dim::Int) = pm.dims[dim]

function (pm::ParametrizedMatrix{T})(params::Vector{T}) where {T}
    @assert pm.np == length(params)
    return pm.np == 0 ?
        reshape(pm.f, pm.dims) :
        reshape(pm.f + pm.D * sparse(params), pm.dims)
end #call

parametrize_full(X::Matrix{T}) where {T} =
        ParametrizedMatrix{T}(zeros(length(X)), eye(length(X)), size(X))

function parametrize_diag(x::Vector{T}; sparse=false) where T
    n = length(x)
    f = sparse ? spzeros(n^2, 1) : zeros(n^2)
    D = sparse ? spzeros(n^2, n) : zeros(n^2, n)
    D[[1 + (i-1)*(n^2 + n + 1) for i in 1:n]] = 1
    return ParametrizedMatrix{T}(f, D, (n,n))
end #parametrize_diag

parametrize_none(X::Matrix{T}) where {T} =
        ParametrizedMatrix{T}(vec(X), zeros(length(X), 0), size(X))

parametrize_none(X::SparseMatrixCSC{T}) where {T} =
        ParametrizedMatrix{T}(vec(X), spzeros(length(X), 0), size(X))

# Time-independent parametrized state space model
struct ParametrizedSSM{T} <: AbstractStateSpaceModel{T}

    # Transition equation and noise covariance

    A1::Function
    A2::ParametrizedMatrix{T}
    A3::Function

    B1::Function
    B2::ParametrizedMatrix{T}

    G::Function
    Q::ParametrizedMatrix{T}

    # Observation equation and noise covariance

    C1::Function
    C2::ParametrizedMatrix{T}
    C3::Function

    D1::Function
    D2::ParametrizedMatrix{T}

    H::Function
    R::ParametrizedMatrix{T}

    # Initial state and error covariance
    x1::ParametrizedMatrix{T}
    J::AbstractMatrix
    S::ParametrizedMatrix{T}

    nx::Int
    ny::Int
    nu::Int
    nq::Int
    nr::Int
    ns::Int

    function ParametrizedSSM{T}(A1::Function, A2::ParametrizedMatrix{T}, A3::Function,
                                B1::Function, B2::ParametrizedMatrix{T},
                                G::Function, Q::ParametrizedMatrix{T},
                                C1::Function, C2::ParametrizedMatrix{T}, C3::Function,
                                D1::Function, D2::ParametrizedMatrix{T},
                                H::Function, R::ParametrizedMatrix{T},
                                x1::ParametrizedMatrix{T}, J::AbstractMatrix, S::ParametrizedMatrix{T}) where {T}
        nx, ny, nu, nq, nr, ns = confirm_matrix_sizes(A1(1), A2, A3(1), B1(1), B2, G(1), Q,
                                              C1(1), C2, C3(1), D1(1), D2, H(1), R, x1, J, S)
        new(A1, A2, A3, B1, B2, G, Q, C1, C2, C3, D1, D2, H, R, x1, J, S, nx, ny, nu, nq, nr, ns)
    end

end #ParametrizedSSM

ParametrizedSSM(A2::ParametrizedMatrix{T}, Q::ParametrizedMatrix{T},
                       C2::ParametrizedMatrix{T}, R::ParametrizedMatrix{T},
                       x1::ParametrizedMatrix{T}, S::ParametrizedMatrix{T};
                       A1::Function=_->speye(size(A2,1)), A3::Function=_->speye(size(A2,2)),
                       B1::Function=_->speye(size(x1,1)),
                       B2::ParametrizedMatrix{T}=parametrize_none(spzeros(size(B1(1),2), 1)),
                       G::Function=_->speye(size(x1,1)),
                       C1::Function=_->speye(size(C2, 1)), C3::Function=_->speye(size(C2,2)),
                       D1::Function=_->speye(size(C1(1),1)),
                       D2::ParametrizedMatrix{T}=parametrize_none(spzeros(size(C1(1),1), 1)),
                       H::Function=_->speye(size(C1(1),1)), J::AbstractMatrix=speye(size(x1, 1))) where {T} =
          ParametrizedSSM{T}(A1, A2, A3, B1, B2, G, Q, C1, C2, C3, D1, D2, H, R, x1, J, S)

# State space model parameters
struct SSMParameters{T}

    # Process transition and noise covariance
    A::Vector{T}
    B::Vector{T}
    Q::Vector{T}

    # Observation and noise covariance
    C::Vector{T}
    D::Vector{T}
    R::Vector{T}

    # Initial state and error covariance
    x1::Vector{T}
    S::Vector{T}

end #SSMParameters

SSMParameters(__::T; A::Vector{T}=T[], B::Vector{T}=T[], Q::Vector{T}=T[],
               C::Vector{T}=T[], D::Vector{T}=T[], R::Vector{T}=T[],
               x1::Vector{T}=T[], S::Vector{T}=T[]) where {T} =
              SSMParameters{T}(A, B, Q, C, D, R, x1, S)

function (m::ParametrizedSSM{T})(p::SSMParameters{T}) where {T}
    A2, B2, Q = m.A2(p.A), m.B2(p.B), m.Q(p.Q)
    C2, D2, R = m.C2(p.C), m.D2(p.D), m.R(p.R)
    A(t) = m.A1(t) * A2 * m.A3(t)
    B(t) = m.B1(t) * B2
    V(t) = m.G(t) * Q * m.G(t)'
    C(t) = m.C1(t) * C2 * m.C3(t)
    D(t) = m.D1(t) * D2
    W(t) = m.H(t) * R * m.H(t)'
    x1 = vec(m.x1(p.x1))
    P1 = m.J * m.S(p.S) * m.J'
    return StateSpaceModel(A, V, C, W, x1, P1, B=B, D=D)
end #call

function confirm_matrix_sizes(F::AbstractMatrix, B::AbstractMatrix, V::AbstractMatrix,
                              G::AbstractMatrix, D::AbstractMatrix, W::AbstractMatrix,
                              x1::Vector, P1::AbstractMatrix)

    nx = size(F, 1)
    nu = size(B, 2)
    ny = size(G, 1)

    @assert size(F) == (nx, nx)
    @assert size(B) == (nx, nu)
    @assert size(V) == (nx, nx)

    @assert size(G) == (ny, nx)
    @assert size(D) == (ny, nu)
    @assert size(W) == (ny, ny)

    @assert size(x1) == (nx,) || size(x1) == (nx, 1)
    @assert size(P1) == (nx, nx)

    return nx, ny, nu

end #confirm_matrix_sizes

function confirm_matrix_sizes(A1::AbstractMatrix, A2::ParametrizedMatrix{T}, A3::AbstractMatrix,
                           B1::AbstractMatrix, B2::ParametrizedMatrix{T},
                           G::AbstractMatrix, Q::ParametrizedMatrix{T},
                           C1::AbstractMatrix, C2::ParametrizedMatrix{T}, C3::AbstractMatrix,
                           D1::AbstractMatrix, D2::ParametrizedMatrix{T},
                           H::AbstractMatrix, R::ParametrizedMatrix{T},
                           x1::ParametrizedMatrix{T}, J::AbstractMatrix, S::ParametrizedMatrix{T}) where T

    @assert size(B2, 2) == size(D2, 2)

    nx = size(A1, 1)
    ny = size(C1, 1)
    nu = size(B2, 2)

    na1, na2  = size(A1, 2), size(A2, 2)
    nb        = size(B1, 2)
    nc1, nc2  = size(C1, 2), size(C2, 2)
    nd        = size(D1, 2)

    nq = size(Q, 1)
    nr = size(R, 1)
    ns = size(S, 1)

    @assert size(A1) == (nx, na1)
    @assert size(A2) == (na1, na2)
    @assert size(A3) == (na2, nx)

    @assert size(B1) == (nx, nb)
    @assert size(B2) == (nb, nu)

    @assert size(G)  == (nx, nq)
    @assert size(Q)  == (nq, nq)

    @assert size(C1) == (ny, nc1)
    @assert size(C2) == (nc1, nc2)
    @assert size(C3) == (nc2, nx)

    @assert size(D1) == (ny, nd)
    @assert size(D2) == (nd, nu)

    @assert size(H)  == (ny, nr)
    @assert size(R)  == (nr, nr)

    @assert length(x1)  == nx
    @assert size(J)     == (nx, ns)
    @assert size(S)     == (ns, ns)

    return nx, ny, nu, nq, nr, ns

end #confirm_matrix_sizes

"
Generates a realization of a state space model.

Arguments:
  model : StateSpaceModel
  Model defining the process

  n : Int
  Number of steps to simulate.

  u : Array (optional)
  n x nu array of exogenous inputs
"
function simulate(model::StateSpaceModel{T}, n::Int; u::Array{T}=zeros(n, model.nu)) where T
    @assert size(u, 1) == n
    @assert size(u, 2) == model.nu

    # create empty arrays to hold the state and observed series
    x = zeros(model.nx, n)
    y = zeros(model.ny, n)
    u = u'

    # Cholesky decompositions of the covariance matrices, for generating
    # random noise
    V_chol(t) = all(model.V(t) .== 0) ? model.V(t) : chol(Hermitian(model.V(t), :L))
    W_chol(t) = all(model.W(t) .== 0) ? model.W(t) : chol(Hermitian(model.W(t), :L))

    # Generate the series
    x[:, 1] = model.x1
    y[:, 1] = model.C(1) * model.x1 + model.D(1) * u[:, 1] + W_chol(1) * randn(model.ny)
    for i=2:n
        x[:, i] = model.A(i-1) * x[:, i-1] + model.B(i-1) * u[:, i-1] + V_chol(i-1) * randn(model.nx)
        y[:, i] = model.C(i)   * x[:, i]   + model.D(i)   * u[:, i]   + W_chol(i) * randn(model.ny)
    end

    return x', y'
end


text:
using YaoBase, YaoArrayRegister, YaoBlocks, LightGraphs
using YaoBase: @interface
export GraphReg, nqubits, nactive, addbits!, add_edge!

struct GraphReg{B, T, Graph <: AbstractGraph} <: AbstractRegister{B, T}
    graphs::NTuple{B, Graph}

    function GraphReg(::Type{T}, graphs::NTuple{B, G}) where {B, T, G <: AbstractGraph}
        # check if each graph has the same number of vertices
        n = first(graphs)
        for k in 2:B
            n == nv(graphs[k]) || error("expect all input graph has the same number of vertices, got $(nv(graphs[k]))")
        end
        new{B, T, G}(graphs)
    end
end

"""
    GraphReg([T=ComplexF64], graphs...)

Create a `GraphReg` with a set of graphs. Graphs should be created as a concrete
type of [`AbstractGraph`](@ref) in [LightGraphs.jl](https://github.com/JuliaGraphs/LightGraphs.jl).
"""
GraphReg(graphs::NTuple{B, G}) where {B, G} = GraphReg(ComplexF64, graphs)

GraphReg(::Type{T}, graphs::AbstractGraph...) where T = GraphReg(T, graphs)
GraphReg(graphs::AbstractGraph...) = GraphReg(ComplexF64, graphs)

# generate a random graph state
"""
    GraphReg{B}([T=ComplexF64], nv[, ne=0])
    GraphReg([T=ComplexF64], nv[, ne=0])

Create a `GraphReg` of random graphs with batch size `B` and number of vertices `nv`
number of edges `ne`.
"""
GraphReg(::Type{T}, nv::Integer, ne::Integer=0) where T = GraphReg{1}(T, nv, ne)
GraphReg(nv::Integer, ne::Integer=0) = GraphReg{1}(nv, ne)

GraphReg{B}(::Type{T}, nv::Integer, ne::Integer=0) where {T, B} = GraphReg(ntuple(k->SimpleGraph(nv, ne), B))
GraphReg{B}(nv::Integer, ne::Integer=0) where B = GraphReg{B}(ComplexF64, nv, ne)

YaoBase.nqubits(r::GraphReg) = nv(r.graphs[1])
YaoBase.nactive(r::GraphReg) = nqubits(r)

function YaoBase.addbits!(r::GraphReg, n::Int)
    for _ in 1:n
        foreach(add_vertex!, r.graphs)
    end
    return r
end

"""
    czcircuit(n, graph)

Return the Control-Z circuit to construct graph state from ``|+⟩^V``
"""
@interface function czcircuit(n, graph::AbstractGraph)
    circuit = chain(n)
    for E in edges(graph)
        # CZ
        push!(circuit, control(E.dst, E.src=>Z))
    end
    return circuit
end

# convert to normal array state
function YaoArrayRegister.ArrayReg(r::GraphReg{1, T}) where T
    # CZ * H |0>
    n = nqubits(r)
    out = zero_state(T, n)
    # H |0> => |+>
    out |> repeat(n, H)
    # CZ |+> => G
    out |> czcircuit(n, r.graphs[1])
    return out
end

function YaoArrayRegister.ArrayReg(r::GraphReg{B, T}) where {B, T}
    n = nqubits(r)
    out = zero_state(T, n; nbatch=B)
    out |> repeat(n, H)

    for k in 1:B
        viewbatch(out, k) |> czcircuit(n, r.graphs[k])
    end
    return out
end

"""
    add_edge!(regsiter, x, y)

Add edge to a `GraphReg` with given vertices `x` and `y`.

!!! note

    For register with batch size > 1, this edge will be added to all the graphs
    stored in the register.
"""
LightGraphs.add_edge!(r::GraphReg{1}, x, y) = add_edge!(r.graphs[1], x, y)

function LightGraphs.add_edge!(r::GraphReg, x, y)
    for each in r.graphs
        add_edge!(each, x, y)
    end
    return r
end

max_stars_repo_name: QuantumBFS/YaoGraphs

max_issues_repo_name: QuantumBFS/YaoGraphs

max_forks_repo_name: QuantumBFS/YaoGraphs


text:
struct Rates
    interior::Matrix{Float64}
    lwr::Vector{Float64}
    upr::Vector{Float64}
end

"""
Indexing for Rates by cell, k, then phase, i.
"""
function Base.getindex(s::Rates, k::Int, i::Int)
    if k==0
        r = s.lwr[i]
    elseif k==size(s.interior,2)+1
        r = s.upr[i]
    else
        r = s.interior[i,k]
    end
    return r
end
Base.getindex(s::Rates, i::AbstractVector{Int}, k::AbstractVector{Int}) =
    [s[ii,kk] for ii in i, kk in k]
Base.getindex(s::Rates, i::AbstractVector{Int}, k::Int) =
    [s[ii,k] for ii in i]
Base.getindex(s::Rates, i::Int, k::AbstractVector{Int}) =
    [s[i,kk] for kk in k]

struct FluidFluidQueue{T<:DiscretisedFluidQueue}# where N<:AbstractVector{Float64}}
    dq::T
    rates::Rates
    function FluidFluidQueue{T}(dq::T,rates::Rates) where T<:DiscretisedFluidQueue
        if !(n_intervals(dq)==size(rates.interior,2))
            throw(DomainError("size(rates.interior,2) needs to match the number of intervals in the DiscretisedFluidQueue"))
        elseif !(n_phases(dq)==size(rates.interior,1))
            throw(DomainError("size(rates.interior,1) needs to match the number of phases in the DiscretisedFluidQueue"))
        end
        if !(sum(DiscretisedFluidQueues._has_left_boundary(dq.model.S))==length(rates.lwr))
            throw(DomainError("rates.lwr must have length equal to the number of phases with
                the a left point mass, sum(DiscretisedFluidQueues._has_left_boundary(B.dq.model.S))"))
        end
        if !(sum(DiscretisedFluidQueues._has_right_boundary(dq.model.S))==length(rates.upr))
            throw(DomainError("rates.upr must be have length equal to the number of phases with
                the a left point mass, sum(DiscretisedFluidQueues._has_left_boundary(B.dq.model.S))"))
        end
        return new{T}(dq,rates)
    end
end
FluidFluidQueue(dq::DiscretisedFluidQueue,rates::Rates)=FluidFluidQueue{typeof(dq)}(dq,rates)

function FluidFluidQueue(dq::DiscretisedFluidQueue,interior::Matrix{Float64},lwr::Vector{Float64},upr::Vector{Float64})
    return FluidFluidQueue(dq,Rates(interior,lwr,upr))
end

function getrates(ffq::FluidFluidQueue,i::Int,k::Int)
    if k==0
        if rates(ffq.dq,i)<=0.0
            idx = sum(DiscretisedFluidQueues._has_left_boundary(ffq.dq.model.S[1:i]))
            r = ffq.rates.lwr[idx]
        else
            r = 0.0
        end
    elseif k==size(ffq.rates.interior,2)+1
        if rates(ffq.dq,i)<=0.0
            idx = sum(DiscretisedFluidQueues._has_right_boundary(ffq.dq.model.S[1:i]))
            r = ffq.rates.upr[idx]
        else
            r = 0.0
        end
    else
        r = ffq.rates.interior[i,k]
    end
    return r
end
getrates(ffq::FluidFluidQueue,i::AbstractVector{Int},k::AbstractVector{Int}) =
    [getrates(ffq::FluidFluidQueue,ii,kk) for ii in i, kk in  k]
getrates(ffq::FluidFluidQueue,i::Int,k::AbstractVector{Int}) =
    [getrates(ffq::FluidFluidQueue,i,kk) for kk in  k]
getrates(ffq::FluidFluidQueue,i::AbstractVector{Int},k::Int) =
    [getrates(ffq::FluidFluidQueue,ii,k) for ii in i]

function DiscretisedFluidQueues.augment_model(ffq::FluidFluidQueue)
    aug_model = augment_model(ffq.dq.model)
    # aug_rates = zeros(n_phases(augment_model),n_intervals(ffq.dq))
    interior_rates = zeros(0,n_intervals(ffq.dq))
    for i in 1:n_phases(ffq.dq.model)
        if rates(ffq.dq)[i]==0.0
            interior_rates = [interior_rates; repeat(ffq.rates.interior[[i],:],2)]
        else
            interior_rates = [interior_rates; ffq.rates.interior[[i],:]]
        end
    end
    aug_dq = DiscretisedFluidQueue(aug_model,ffq.dq.mesh)
    return FluidFluidQueue(aug_dq,interior_rates,ffq.rates.lwr,ffq.rates.upr)
end


text:
export elitesel, elitesel!

# TODO: Implement for `MetaPopulation`.

@doc """elitesel(p::Population, q::Int64=1)

Create a new population through elite selection, which keeps only the `q`
fittest members of the population and copies them to create the new population.

If `q > 1`, the elite individuals will be copied proportionally based on their
fitnesses. So, as an example, assume the fittest individual has fitness 0.75,
the second fittest individual has fitness 0.5, and the population size is 100.
The fittest individual will be copied 60 times and the other will be copied 40
times to make the new population.

Each elite individual will always show up in the resulting population at least
once, even when its fitness is so low that its share in the population would
round to zero.
"""
function elitesel(p::Population, q::Int64=1)
  np = Population(p)
  elitesel!(np, q)
  return np
end

@doc """elitesel!(p::Population, q::Int64=1)

Conduct elite selection in-place.
"""
function elitesel!(p::Population, q::Int64=1)
  elites = sort(p.genotypes, by=(g) -> fitness(g))[end-q+1:end]
  fits = [fitness(g) for g = elites]
  counts = [round(Int64, f / sum(fits) * popsize(p) |> ceil) for f = fits]
  # We need to adjust counts to achieve the same population size, we do this by
  # penalizing the most-fit individual. This is a naive approach, but it should
  # work fine for most use cases. It would be better to ensure that the
  # most-fit individual always ends up more common than the next most-fit, etc.
  while sum(counts) > popsize(p)
    counts[end] -= 1
  end
  next = 1
  for i = 1:length(elites)
    for _ = 1:counts[i]
      p.genotypes[next] = Genotype(elites[i])
      next += 1
    end
  end
end


text:
const func_to_jamsd = Dict(
:ceil => (1,1),
:floor => (2,1),
:round => (3,1),
:fmod => (4,1),
:trunc => (5,1),
:sign => (6,1),
:min => (7,2),
:max => (8,2),
# :sqr => (9,1),
:exp => (10,1),
:log => (11,1),
:log10 => (12,1),
:sqrt => (13,1),
:abs => (14,1),
:cos => (15,1),
:sin => (16,1),
:atan => (17,1),
:erf => (18,1),
# (dunfm),
# (dnorm),
:^ => (21,2), #power
# (jdate),
# (jtime),
# (jstart),
# (jnow),
# (error),
# (gyear),
# (gmonth),
# (gday),
# (gdow),
# (gleap),
# (ghour),
# (gminute),
# (gsecond),
# (curseed),
# (timest),
# (timeco),
# (timeex),
# (timecl),
:frac => (40,1),
# (errorl),
# (heaps),
# (fact),
# (unfmi),
# (pi),
# (ncpf),
# (ncpcm),
# (entropy),
# (sigmoid),
:log2 => (50,1),
# (boolnot),
# (booland),
# (boolor),
# (boolxor),
# (boolimp),
# (booleqv),
:(==)   => (57,2),  # (relopeq),
:(>)    => (58,2),  # (relopgt),
:(>=)   => (59,2),  # (relopge),
:(<)    => (60,2),  # (reloplt),
:(<=)   => (61,2),  # (relople),
:(!=)   => (62,2),  # (relopne),
:ifelse => (63,3),  # (ifthen),
#   &_rpower => (64,1),
# (edist),
# (div),
# (div0),
# (sllog10),
# (sqlog10),
# (slexp),
# (sqexp),
# (slrec),
# (sqrec),
# (cvpower),
# (vcpower),
# (centropy),
# (gmillisec),
# (maxerror),
# (timeel),
# (gamma),
# (loggamma),
# (beta),
# (logbeta),
# (gammareg),
# (betareg),
:sinh => (86,1),
:cosh => (87,1),
:tanh => (88,1),
# (mathlastrc),
# (mathlastec),
# (mathoval),
# (signpower),
# (handle),
# (ncpvusin),
# (ncpvupow),
# (binomial),
# (rehandle),
# (gamsver),
# (delhandle),
:tan => (100,1),
:acos => (101,1),
:asin => (102,1),
:atan2 => (103,1),
# (sleep),
# (heapf),
# (cohandle),
# (gamsrel),
# (poly),
# (licensestatus),
# (licenselevel),
# (heaplimit),
# (dummy),
)

const sense_to_jamsd = Dict(
:Min => 0,
:Max => 1,
)

const relation_to_jamsd = Dict(
:multiple => -1,
:<=       => 2,
:>=       => 1,
:(==)     => 0,
'='       => 0,
)

const quad_relation_sense = Dict(
'=' => '=',
'<' => :<=,
'>' => :>=
)

const arithm_ops = Dict(
:+ => 2,
:sum => 2,
:- => 3,
:* => 4,
:/ => 5,
)

const reverse_relation = Dict(
:0 => 0,
:1 => 2,
:2 => 1,
:4 => 4,
)

const nary_functions = Set{Symbol}([
:min,
:max,
# :count,
# :numberof,
# :numberofs,
# :and_n,
# :or_n,
# :alldiff,
])


text:
using DelimitedFiles

function solve(input)
    raw = readdlm(input,'\n', Int32)
    n_recipes=raw[1]
    @show n_recipes
    @time @show hotchocolate(n_recipes)
    @time @show hotchocolate_b(string(n_recipes))
end

function hotchocolate(n_recipes)
    scores=[3,7]
    elves=[1,2]
    for iteration in 1:n_recipes+10
        append!(scores,reverse(digits(sum(scores[elves]))))
        elves=((elves+(scores[elves])).%length(scores)).+1
        #@show scores
    end
    ten_scores = join(scores[n_recipes+1:n_recipes+10])
    return ten_scores
end

@assert hotchocolate(9)=="5158916779"
@assert hotchocolate(5)=="0124515891"
@assert hotchocolate(18)=="9251071085"
@assert hotchocolate(2018)=="5941429882"

function hotchocolate_b(input_str)
    scores=[3,7]
    elves=[1,2]
    while !occursin(input_str,join(view(scores,max(1,length(scores)-10):length(scores))))
        new_score=reverse(digits(sum(scores[elves])))
        append!(scores,new_score)
        elves=((elves+(scores[elves])).%length(scores)).+1
    end
    sequence = map(c->parse(Int8,c),collect(input_str))
    result = findfirst(input_str,join(scores))[1]-1
    return result
end

@assert hotchocolate_b("51589")==9
@assert hotchocolate_b("01245")==5
@assert hotchocolate_b("92510")==18
@assert hotchocolate_b("59414")==2018

solve("input.txt")