text:
############################################################################################
# CREATING
############################################################################################

facts("Creating a VirtualArray") do
    context("normal case") do
        # set up
        num = rand(1:1000)
        len = rand(1:100)

        a = collect(num:num+len)
        b = collect(num:num+len)
        expected = cat(1,a,b)

        test = VirtualArray{Int, 1}(a,b)
        test_2 = VirtualArray{Int, 1}(1,a,b)

        @fact test.parents[1] --> a
        @fact test.parents[2] --> b
        @fact test --> expected
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> test_2
        @fact eachindex(test) --> eachindex(expected)
    end
    context("one parent") do
        a = collect(1:9)
        test = VirtualArray{Int, 1}(a)
        test_2 = VirtualArray{Int, 1}(1, a)
        expected = cat(1,a)

        @fact test.parents[1] --> a
        @fact length(test.parents) --> 1
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> test_2
        @fact eachindex(test) --> eachindex(expected)
    end
    context("multiple parent") do
        # set up
        num = rand(1:1000)
        num_parents = rand(3:10)
        len = rand(1:100)

        parents = []
        for i in 1:num_parents
            push!(parents, collect(num:num+len-1))
        end

        expected = cat(1, parents...)
        test = VirtualArray{Int, 1}(parents...)
        test_2 = VirtualArray{Int, 1}(1, parents...)

        @fact test.parents --> parents
        @fact length(test.parents) --> num_parents
        @fact test --> expected
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> test_2
        @fact eachindex(test) --> eachindex(expected)
    end
    context("1 2 dimensional parent") do
        # set up
        len = rand(2:4)

        a = rand(len,len)

        expected = cat(1, a)
        test = VirtualArray{Float64, 2}(a)

        @fact test.parents[1] --> a
        @fact length(test.parents) --> 1
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> expected
        @fact eachindex(test) --> eachindex(expected)
    end
    context("2 dimensional parents") do
        # set up
        len = rand(2:4)

        a = rand(len,len)
        b = rand(len,len)

        expected = cat(1, a, b)
        test = VirtualArray{Float64, 2}(a, b)

        @fact test.parents[1] --> a
        @fact test.parents[2] --> b
        @fact length(test.parents) --> 2
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> expected
        @fact eachindex(test) --> eachindex(expected)
    end
    context("multiple 2 dimensional parents") do
        # set up
        num_parents = rand(3:10)
        len = rand(1:100)

        parents = []
        for i in 1:num_parents
            push!(parents, rand(len,len))
        end

        expected = cat(1, parents...)
        test = VirtualArray{Float64, 2}(parents...)

        @fact length(test.parents) --> num_parents
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> expected
        @fact eachindex(test) --> eachindex(expected)
    end
    context("1 multi dimensional parents") do
        # set up
        # keep these numbers small because we can run out of memory or get very slow tests
        num_dims = rand(3:6) # no larger than 6
        len = rand(2:5) # no larger than 5

        dims = zeros(Int, num_dims) + len

        a = rand(dims...)

        expected = cat(1, a)
        test = VirtualArray{Float64, num_dims}(a)

        @fact test.parents[1] --> a
        @fact length(test.parents) --> 1
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> expected
        @fact eachindex(test) --> eachindex(expected)
    end
    context("2 multi dimensional parents") do
        # set up
        # keep these numbers small because we can run out of memory or get very slow tests
        num_dims = rand(3:6) # no larger than 6
        len = rand(2:5) # no larger than 5

        dims = zeros(Int, num_dims) + len

        a = rand(dims...)
        b = rand(dims...)

        expected = cat(1, a, b)
        test = VirtualArray{Float64, num_dims}(a, b)

        @fact test.parents[1] --> a
        @fact test.parents[2] --> b
        @fact length(test.parents) --> 2
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> expected
        @fact eachindex(test) --> eachindex(expected)
    end
    context("multi multi dimensional parents") do
        # set up
        # keep these numbers small because we can run out of memory or get very slow tests
        num_parents = rand(3:10) # no larger than 10
        num_dims = rand(3:6) # no larger than 6
        len = rand(2:5) # no larger than 5

        dims = zeros(Int, num_dims) + len

        parents = []
        for i in 1:num_parents
            push!(parents, rand(dims...))
        end

        expected = cat(1, parents...)
        test = VirtualArray{Float64, num_dims}(parents...)

        @fact test.parents --> parents
        @fact length(test.parents) --> num_parents
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> expected
        @fact eachindex(test) --> eachindex(expected)
    end
end

facts("Creating VirtualArrays and expanding on different dimensions") do
    context("1 2 dimensional parent") do
        # set up
        len = rand(2:4)
        expanded_dim = 2

        a = rand(len,len)

        expected = cat(expanded_dim, a)
        test = VirtualArray{Float64, 2}(expanded_dim, a)

        @fact test.parents[1] --> a
        @fact length(test.parents) --> 1
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> expected
        @fact eachindex(test) --> eachindex(expected)
    end
    context("2 dimensional parents") do
        # set up
        len = rand(2:4)
        expanded_dim = 2

        a = rand(len,len)
        b = rand(len,len)

        expected = cat(expanded_dim, a, b)
        test = VirtualArray{Float64, 2}(expanded_dim, a, b)

        @fact test.parents[1] --> a
        @fact test.parents[2] --> b
        @fact length(test.parents) --> 2
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> expected
        @fact eachindex(test) --> eachindex(expected)
    end
    context("multiple 2 dimensional parents") do
        # set up
        num_parents = rand(3:10)
        len = rand(1:100)
        expanded_dim = 2

        parents = []
        for i in 1:num_parents
            push!(parents, rand(len,len))
        end

        expected = cat(expanded_dim, parents...)
        test = VirtualArray{Float64, 2}(expanded_dim, parents...)

        @fact length(test.parents) --> num_parents
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> expected
        @fact eachindex(test) --> eachindex(expected)
    end
    context("1 multi dimensional parents") do
        # set up
        # keep these numbers small because we can run out of memory or get very slow tests
        num_dims = rand(3:6) # no larger than 6
        len = rand(2:5) # no larger than 5
        expanded_dim = rand(3:num_dims)

        dims = zeros(Int, num_dims) + len

        a = rand(dims...)

        expected = cat(expanded_dim, a)
        test = VirtualArray{Float64, num_dims}(expanded_dim, a)

        @fact test.parents[1] --> a
        @fact length(test.parents) --> 1
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> expected
        @fact eachindex(test) --> eachindex(expected)
    end
    context("2 multi dimensional parents") do
        # set up
        # keep these numbers small because we can run out of memory or get very slow tests
        num_dims = rand(3:6) # no larger than 6
        len = rand(2:5) # no larger than 5
        expanded_dim = rand(3:num_dims)

        dims = zeros(Int, num_dims) + len

        a = rand(dims...)
        b = rand(dims...)

        expected = cat(expanded_dim, a, b)
        test = VirtualArray{Float64, num_dims}(expanded_dim, a, b)

        @fact test.parents[1] --> a
        @fact test.parents[2] --> b
        @fact length(test.parents) --> 2
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> expected
        @fact eachindex(test) --> eachindex(expected)
    end
    context("multi multi dimensional parents") do
        # set up
        # keep these numbers small because we can run out of memory or get very slow tests
        num_parents = rand(3:10) # no larger than 10
        num_dims = rand(3:6) # no larger than 6
        len = rand(2:5) # no larger than 5
        expanded_dim = rand(3:num_dims)

        dims = zeros(Int, num_dims) + len

        parents = []
        for i in 1:num_parents
            push!(parents, rand(dims...))
        end

        expected = cat(expanded_dim, parents...)
        test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

        @fact test.parents --> parents
        @fact length(test.parents) --> num_parents
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> expected
        @fact eachindex(test) --> eachindex(expected)
    end

    context("creating an 1 d VirtualArray with 1 d parents of different lengths") do
        # set up
        num_dims = 1
        len = rand(2:2:20) # no larger than 5
        len_2 = rand(3:2:30) # no larger than 5
        expanded_dim = 1

        a = rand(len)
        b = rand(len_2)
        parents = Array[a, b]

        expected = cat(expanded_dim, parents...)
        test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

        @fact test.parents --> parents
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> expected
        @fact eachindex(test) --> eachindex(expected)
    end
    context("creating an 2 d VirtualArray with 2 d parents of different lengths") do
        # set up
        num_dims = 2
        len = rand(2:2:20) # no larger than 5
        len_2 = rand(3:2:21) # no larger than 5
        expanded_dim = 2

        a = rand(len, len)
        b = rand(len, len_2)
        parents = Array[a, b]

        expected = cat(expanded_dim, parents...)
        test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

        @fact test.parents --> parents
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> expected
        @fact eachindex(test) --> eachindex(expected)
    end
    context("creating an 2 d VirtualArray with 2 d parents of different heights") do
        # set up
        num_dims = 2
        len = rand(2:2:20) # no larger than 5
        len_2 = rand(3:2:21) # no larger than 5
        expanded_dim = 1

        a = rand(len, len)
        b = rand(len_2, len)
        parents = Array[a, b]

        expected = cat(expanded_dim, parents...)
        test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

        @fact test.parents --> parents
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> expected
        @fact eachindex(test) --> eachindex(expected)
    end

    context("creating an 2 d VirtualArray with 1 d parents") do
        # set up
        num_dims = 2
        len = rand(2:10) # no larger than 5
        expanded_dim = 2

        a = rand(len)
        b = rand(len)
        parents = Array[a, b]

        expected = cat(expanded_dim, parents...)
        test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

        @fact test.parents --> parents
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> expected
        @fact eachindex(test) --> eachindex(expected)
    end
    context("creating an 3> d VirtualArray with 1 d parents") do
        # set up
        num_dims = rand(3:10)
        len = rand(2:10) # no larger than 5
        expanded_dim = num_dims

        a = rand(len)
        b = rand(len)
        parents = Array[a, b]

        expected = cat(expanded_dim, parents...)
        test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

        @fact test.parents --> parents
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> expected
        @fact eachindex(test) --> eachindex(expected)
    end
end

facts("Creating VirtualArrays oddly") do
    context("storing float array into a real VirtualArray") do
        # set up
        num_dims = 1
        len = rand(2:10)
        expanded_dim = 1

        a = rand(Float64, len)
        b = rand(Float64, len)
        parents = Array[a, b]

        expected = cat(expanded_dim, parents...)
        test = VirtualArray{Real, num_dims}(expanded_dim, parents...)

        @fact test.parents --> parents
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact test --> expected
        @fact eachindex(test) --> eachindex(expected)
    end
end

############################################################################################
# MODIFYING
############################################################################################

facts("Modifying values in a VirtualArray with 1 d arrays") do
    context("normal case changing one VirtualArray element in the first parent") do

        # set up
        num = rand(1:1000)
        len = rand(1:100)
        num_picked = rand(1:1000)
        index_picked = rand(1:len)

        a = collect(num:num+len)
        b = collect(num:num+len)

        expected = cat(1, a, b)
        test = VirtualArray{Int, 1}(a,b)

        test[index_picked] = num_picked
        expected[index_picked] = num_picked

        @fact test[index_picked] --> num_picked
        @fact a[index_picked] --> num_picked
        @fact b[index_picked] --> num + index_picked - 1
        @fact test --> expected

    end
    context("normal case changing one parent element in the first parent") do

        # set up
        num = rand(1:1000)
        len = rand(1:100)
        num_picked = rand(1:1000)
        index_picked = rand(1:len)

        a = collect(num:num+len)
        b = collect(num:num+len)

        expected = cat(1, a, b)
        test = VirtualArray{Int, 1}(a,b)

        a[index_picked] = num_picked
        expected = cat(1, a, b)

        @fact test[index_picked] --> num_picked
        @fact a[index_picked] --> num_picked
        @fact b[index_picked] --> num + index_picked - 1
        @fact test --> expected

    end
    context("normal case changing one VirtualArray element in the any parent") do

        # set up
        num = rand(1:1000)
        num_parents = rand(3:10)
        len = rand(1:10)
        change_p = rand(3:num_parents)
        change_i = rand(1:len)
        change_to = rand(1:10)

        parents = []
        for i in 1:num_parents
            push!(parents, collect(num:num+len-1))
        end

        expected = cat(1, parents...)
        test = VirtualArray{Int, 1}(parents...)

        test[(change_p-1)*len+change_i] = change_to
        expected[(change_p-1)*len+change_i] = change_to

        @fact test[(change_p-1)*len+change_i] --> change_to
        @fact test.parents[change_p][change_i] --> change_to
        @fact test --> expected

    end
    context("normal case changing one element in the any parent") do

        # set up
        num = rand(1:1000)
        num_parents = rand(3:10)
        len = rand(1:100)
        change_p = rand(3:num_parents)
        change_i = rand(1:len)
        change_to = rand(1:10)

        parents = []
        for i in 1:num_parents
            push!(parents, collect(num:num+len-1))
        end

        expected = cat(1, parents...)
        test = VirtualArray{Int, 1}(parents...)

        parents[change_p][change_i] = change_to
        expected = cat(1, parents...)

        @fact test.parents[change_p][change_i] --> change_to
        @fact test[(change_p-1)*len+change_i] --> change_to
        @fact test --> expected

    end
    context("changing a value when the parent is a subarray") do

        # set up
        len = rand(1:100)
        change_to = rand(1:1000)
        change_i = rand(1:len)
        parent = rand(Int, len)
        sub_array = sub(parent, 1:len)

        test = VirtualArray{Int, 1}(sub_array)

        test[change_i] = change_to

        @fact test.parents[1][change_i] --> change_to
        @fact parent[change_i] --> change_to
        @fact sub_array[change_i] --> change_to
        @fact test --> parent
        @fact test --> sub_array

    end
end

facts("Modifying values in a VirtualArray with 2 d arrays") do
    context("normal case changing one VirtualArray element in the first parent") do

        # set up
        len = rand(1:100)
        num_picked = rand(1:1000)
        index_picked = [rand(1:len),rand(1:len)]

        a = rand(len,len)
        b = rand(len,len)

        expected = cat(1, a, b)
        test = VirtualArray{Float64, 2}(a,b)

        test[index_picked...] = num_picked
        expected[index_picked...] = num_picked

        @fact test[index_picked...] --> num_picked
        @fact test[index_picked...] --> num_picked
        @fact a[index_picked...] --> num_picked
        @fact test --> expected

    end
    context("normal case changing one parent element in the first parent") do

        # set up
        len = rand(1:100)
        num_picked = rand(1:1000)
        index_picked = [rand(1:len),rand(1:len)]

        a = rand(len,len)
        b = rand(len,len)

        expected = cat(1, a, b)
        test = VirtualArray{Float64, 2}(a,b)

        a[index_picked...] = num_picked
        expected = cat(1, a, b)

        @fact test[index_picked...] --> num_picked
        @fact a[index_picked...] --> num_picked
        @fact test --> expected

    end
    context("normal case changing one VirtualArray element in the any parent") do

        # set up
        num_parents = rand(3:10)
        len = rand(1:10)
        change_p = rand(3:num_parents)
        change_i = [rand(1:len),rand(1:len)]
        change_to = rand(1:10)
        combined_i = [(change_p-1)*len+change_i[1],change_i[2]]

        parents = []
        for i in 1:num_parents
            push!(parents, rand(len,len))
        end

        expected = cat(1, parents...)
        test = VirtualArray{Float64, 2}(parents...)

        test[combined_i...] = change_to
        expected[combined_i...] = change_to

        @fact test[combined_i...] --> change_to
        @fact test.parents[change_p][change_i...] --> change_to
        @fact test --> expected

    end
    context("normal case changing one element in the any parent") do

        # set up
        num_parents = rand(3:10)
        len = rand(1:100)
        change_p = rand(3:num_parents)
        change_i = [rand(1:len),rand(1:len)]
        change_to = rand(1:10)
        combined_i = [(change_p-1)*len+change_i[1],change_i[2]]

        parents = []
        for i in 1:num_parents
            push!(parents, rand(len,len))
        end

        expected = cat(1, parents...)
        test = VirtualArray{Float64, 2}(parents...)

        parents[change_p][change_i...] = change_to
        expected = cat(1, parents...)

        @fact test.parents[change_p][change_i...] --> change_to
        @fact test[combined_i...] --> change_to
        @fact test --> expected

    end
end

facts("Modifying values in a VirtualArray with mulit d arrays") do
    context("normal case changing one VirtualArray element in the first parent") do

        # set up
        # keep these numbers small because we can run out of memory or get very slow tests
        num_parents = rand(3:10) # no larger than 10
        num_dims = rand(3:6) # no larger than 6
        len = rand(2:5) # no larger than 5

        dims = zeros(Int, num_dims) + len

        num_picked = rand(1:1000)
        index_picked = []
        for i in 1:num_dims
            push!(index_picked, rand(1:len))
        end

        a = rand(dims...)
        b = rand(dims...)

        expected = cat(1, a, b)
        test = VirtualArray{Float64, num_dims}(a,b)

        test[index_picked...] = num_picked
        expected[index_picked...] = num_picked

        @fact test[index_picked...] --> num_picked
        @fact a[index_picked...] --> num_picked
        @fact test --> expected

    end
    context("normal case changing one parent element in the first parent") do

        # set up
        # keep these numbers small because we can run out of memory or get very slow tests
        num_parents = rand(3:10) # no larger than 10
        num_dims = rand(3:6) # no larger than 6
        len = rand(2:5) # no larger than 5

        dims = zeros(Int, num_dims) + len

        num_picked = rand(1:1000)
        index_picked = []
        for i in 1:num_dims
            push!(index_picked, rand(1:len))
        end

        a = rand(dims...)
        b = rand(dims...)

        expected = cat(1, a, b)
        test = VirtualArray{Float64, num_dims}(a,b)

        a[index_picked...] = num_picked
        expected = cat(1, a, b)

        @fact test[index_picked...] --> num_picked
        @fact a[index_picked...] --> num_picked
        @fact test --> expected

    end
    context("normal case changing one VirtualArray element in the any parent") do

        # set up
        # keep these numbers small because we can run out of memory or get very slow tests
        num_parents = rand(3:10) # no larger than 10
        num_dims = rand(3:6) # no larger than 6
        len = rand(2:5) # no larger than 5

        dims = zeros(Int, num_dims) + len

        parents = []
        for i in 1:num_parents
            push!(parents, rand(dims...))
        end
        change_i = []
        for i in 1:num_dims
            push!(change_i, rand(1:len))
        end
        change_p = rand(3:num_parents)
        change_to = rand(1:10)
        combined_i = [(change_p-1)*len+change_i[1],change_i[2:end]...]

        expected = cat(1, parents...)
        test = VirtualArray{Float64, num_dims}(parents...)

        test[combined_i...] = change_to
        expected[combined_i...] = change_to

        @fact test[combined_i...] --> change_to
        @fact test.parents[change_p][change_i...] --> change_to
        @fact test --> expected

    end
    context("normal case changing one element in the any parent") do

        # set up
        # keep these numbers small because we can run out of memory or get very slow tests
        num_parents = rand(3:10) # no larger than 10
        num_dims = rand(3:6) # no larger than 6
        len = rand(2:5) # no larger than 5

        dims = zeros(Int, num_dims) + len

        parents = []
        for i in 1:num_parents
            push!(parents, rand(dims...))
        end
        change_i = []
        for i in 1:num_dims
            push!(change_i, rand(1:len))
        end
        change_p = rand(3:num_parents)
        change_to = rand(1:10)
        combined_i = [(change_p-1)*len+change_i[1],change_i[2:end]...]

        expected = cat(1, parents...)
        test = VirtualArray{Float64, num_dims}(parents...)

        parents[change_p][change_i...] = change_to
        expected = cat(1, parents...)

        @fact test.parents[change_p][change_i...] --> change_to
        @fact test[combined_i...] --> change_to
        @fact test --> expected

    end
end

facts("Modifying values in a VirtualArray with mulit d arrays and expand in different dimensions") do
    context("normal case changing one VirtualArray element in the first parent") do

        # set up
        # keep these numbers small because we can run out of memory or get very slow tests
        num_parents = rand(3:10) # no larger than 10
        num_dims = rand(3:6) # no larger than 6
        len = rand(2:5) # no larger than 5
        expanded_dim = rand(3:num_dims)

        dims = zeros(Int, num_dims) + len

        num_picked = rand(1:1000)
        index_picked = []
        for i in 1:num_dims
            push!(index_picked, rand(1:len))
        end

        a = rand(dims...)
        b = rand(dims...)

        expected = cat(expanded_dim, a, b)
        test = VirtualArray{Float64, num_dims}(expanded_dim, a,b)

        test[index_picked...] = num_picked
        expected[index_picked...] = num_picked

        @fact test[index_picked...] --> num_picked
        @fact a[index_picked...] --> num_picked
        @fact test --> expected

    end
    context("normal case changing one parent element in the first parent") do

        # set up
        # keep these numbers small because we can run out of memory or get very slow tests
        num_parents = rand(3:10) # no larger than 10
        num_dims = rand(3:6) # no larger than 6
        len = rand(2:5) # no larger than 5
        expanded_dim = rand(3:num_dims)

        dims = zeros(Int, num_dims) + len

        num_picked = rand(1:1000)
        index_picked = []
        for i in 1:num_dims
            push!(index_picked, rand(1:len))
        end

        a = rand(dims...)
        b = rand(dims...)

        expected = cat(expanded_dim, a, b)
        test = VirtualArray{Float64, num_dims}(expanded_dim, a,b)

        a[index_picked...] = num_picked
        expected = cat(expanded_dim, a, b)

        @fact test[index_picked...] --> num_picked
        @fact a[index_picked...] --> num_picked
        @fact test --> expected

    end
    context("normal case changing one VirtualArray element in the any parent") do

        # set up
        # keep these numbers small because we can run out of memory or get very slow tests
        num_parents = rand(3:10) # no larger than 10
        num_dims = rand(3:6) # no larger than 6
        len = rand(2:5) # no larger than 5
        expanded_dim = rand(3:num_dims)

        dims = zeros(Int, num_dims) + len

        parents = []
        for i in 1:num_parents
            push!(parents, rand(dims...))
        end
        change_i = []
        for i in 1:num_dims
            push!(change_i, rand(1:len))
        end
        change_p = rand(3:num_parents)
        change_to = rand(1:10)
        combined_i = [
                change_i[1:expanded_dim-1]...,
                (change_p-1)*len+change_i[expanded_dim],
                change_i[expanded_dim+1:end]...
            ]

        expected = cat(expanded_dim, parents...)
        test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

        test[combined_i...] = change_to
        expected[combined_i...] = change_to

        @fact test[combined_i...] --> change_to
        @fact test.parents[change_p][change_i...] --> change_to
        @fact test --> expected

    end
    context("normal case changing one element in the any parent") do

        # set up
        # keep these numbers small because we can run out of memory or get very slow tests
        num_parents = rand(3:10) # no larger than 10
        num_dims = rand(3:6) # no larger than 6
        len = rand(2:5) # no larger than 5
        expanded_dim = rand(3:num_dims)

        dims = zeros(Int, num_dims) + len

        parents = []
        for i in 1:num_parents
            push!(parents, rand(dims...))
        end
        change_i = []
        for i in 1:num_dims
            push!(change_i, rand(1:len))
        end
        change_p = rand(3:num_parents)
        change_to = rand(1:10)
        combined_i = [
                change_i[1:expanded_dim-1]...,
                (change_p-1)*len+change_i[expanded_dim],
                change_i[expanded_dim+1:end]...
            ]

        expected = cat(expanded_dim, parents...)
        test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

        parents[change_p][change_i...] = change_to
        expected = cat(expanded_dim, parents...)

        @fact test.parents[change_p][change_i...] --> change_to
        @fact test[combined_i...] --> change_to
        @fact test --> expected

    end
end

facts("Accessing VirtualArray values oddly") do
    context("accessing the element with a bunch of 1's padded on") do

        # set up
        len = rand(1:50)
        num_picked = rand(1:1000)
        index_picked = [rand(1:len),rand(1:len)]
        oness = ones(Int, rand(1:20))

        a = rand(len,len)
        b = rand(len,len)

        expected = cat(1, a, b)
        test = VirtualArray{Float64, 2}(a,b)

        test[index_picked...] = num_picked
        expected[index_picked...] = num_picked

        @fact test[index_picked...] --> num_picked
        @fact test[index_picked..., oness...] --> num_picked
        @fact a[index_picked...] --> num_picked
        @fact a[index_picked..., oness...] --> num_picked
        @fact test --> expected

    end
    context("setting the element with a bunch of 1's padded on") do

        # set up
        len = rand(1:50)
        num_picked = rand(1:1000)
        index_picked = [rand(1:len),rand(1:len)]
        oness = ones(Int, rand(1:20))

        a = rand(len,len)
        b = rand(len,len)

        expected = cat(1, a, b)
        test = VirtualArray{Float64, 2}(a,b)

        test[index_picked..., oness...] = num_picked
        expected[index_picked..., oness...] = num_picked

        @fact test[index_picked...] --> num_picked
        @fact test[index_picked..., oness...] --> num_picked
        @fact a[index_picked...] --> num_picked
        @fact a[index_picked..., oness...] --> num_picked
        @fact test --> expected

    end
    context("accessing a multi d VirtualArray like a 1 d array when expanded on dimension 1") do

        # set up
        # keep these numbers small because we can run out of memory or get very slow tests
        num_parents = rand(3:10) # no larger than 10
        num_dims = rand(3:6) # no larger than 6
        len = rand(2:5) # no larger than 5
        expanded_dim = 1

        dims = zeros(Int, num_dims) + len

        parents = []
        for i in 1:num_parents
            push!(parents, rand(dims...))
        end

        expected = cat(expanded_dim, parents...)
        test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

        change_to = rand(1:10)
        change_i = rand(1:length(expected))

        test[change_i] = change_to
        expected[change_i] = change_to

        @fact test --> expected

    end
    context("accessing a multi d VirtualArray like a 1 d array when expanded on dimension 2") do

        # set up
        # keep these numbers small because we can run out of memory or get very slow tests
        num_parents = rand(3:10) # no larger than 10
        num_dims = rand(3:6) # no larger than 6
        len = rand(2:5) # no larger than 5
        expanded_dim = 2

        dims = zeros(Int, num_dims) + len

        parents = []
        for i in 1:num_parents
            push!(parents, rand(dims...))
        end

        expected = cat(expanded_dim, parents...)
        test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

        change_to = rand(1:10)
        change_i = rand(1:length(expected))

        test[change_i] = change_to
        expected[change_i] = change_to

        @fact test --> expected

    end
    context("accessing a multi d VirtualArray like a 1 d array when expanded on dimension above 2") do

        # set up
        # keep these numbers small because we can run out of memory or get very slow tests
        num_parents = rand(3:10) # no larger than 10
        num_dims = rand(3:6) # no larger than 6
        len = rand(2:5) # no larger than 5
        expanded_dim = rand(3:num_dims)

        dims = zeros(Int, num_dims) + len

        parents = []
        for i in 1:num_parents
            push!(parents, rand(dims...))
        end

        expected = cat(expanded_dim, parents...)
        test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

        change_to = rand(1:10)
        change_i = rand(1:length(expected))

        test[change_i] = change_to
        expected[change_i] = change_to

        @fact test --> expected

    end
end

facts("Modifying values in a 1 d VirtualArray with 1 d arrays of different lengths") do
    context("modifying a value in the first parent from VirtualArray") do
        # set up
        num_dims = 1
        len = rand(2:2:20) # no larger than 5
        len_2 = rand(3:2:30) # no larger than 5
        expanded_dim = 1
        num_picked = rand(1:1000)
        index_picked = rand(1:len)

        a = rand(len)
        b = rand(len_2)
        parents = Array[a, b]

        expected = cat(expanded_dim, parents...)
        test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

        test[index_picked] = num_picked
        expected[index_picked] = num_picked

        @fact test[index_picked] --> num_picked
        @fact a[index_picked] --> num_picked
        @fact test --> expected
    end
    context("modifying a value in the second parent from VirtualArray") do
        # set up
        num_dims = 1
        len = rand(2:2:20) # no larger than 5
        len_2 = rand(3:2:30) # no larger than 5
        expanded_dim = 1
        num_picked = rand(1:1000)
        index_picked = rand(1:len_2)

        a = rand(len)
        b = rand(len_2)
        parents = Array[a, b]

        expected = cat(expanded_dim, parents...)
        test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

        test[index_picked + len] = num_picked
        expected[index_picked + len] = num_picked

        @fact test[index_picked + len] --> num_picked
        @fact b[index_picked] --> num_picked
        @fact test --> expected
    end
    context("modifying a value in the first parent from the parent") do
        # set up
        num_dims = 1
        len = rand(2:2:20) # no larger than 5
        len_2 = rand(3:2:30) # no larger than 5
        expanded_dim = 1
        num_picked = rand(1:1000)
        index_picked = rand(1:len)

        a = rand(len)
        b = rand(len_2)
        parents = Array[a, b]

        test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

        a[index_picked] = num_picked
        expected = cat(expanded_dim, parents...)

        @fact test[index_picked] --> num_picked
        @fact a[index_picked] --> num_picked
        @fact test --> expected
    end
    context("modifying a value in the second parent from the parent") do
        # set up
        num_dims = 1
        len = rand(2:2:20) # no larger than 5
        len_2 = rand(3:2:30) # no larger than 5
        expanded_dim = 1
        num_picked = rand(1:1000)
        index_picked = rand(1:len_2)

        a = rand(len)
        b = rand(len_2)
        parents = Array[a, b]

        expected = cat(expanded_dim, parents...)
        test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

        b[index_picked] = num_picked
        expected = cat(expanded_dim, parents...)

        @fact test[index_picked + len] --> num_picked
        @fact b[index_picked] --> num_picked
        @fact test --> expected
    end
end

facts("modifying values of m d VirtualArrays with n d parents") do
    context("modifying the value in the first parent from VirtualArray") do
        # set up
        num_dims = 2
        len = rand(2:10) # no larger than 5
        expanded_dim = 2
        num_picked = rand(1:1000)
        index_picked = rand(1:len)

        a = rand(len)
        b = rand(len)
        parents = Array[a, b]

        expected = cat(expanded_dim, parents...)
        test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

        test[index_picked, 1] = num_picked
        expected[index_picked, 1] = num_picked

        @fact test[index_picked, 1] --> num_picked
        @fact a[index_picked] --> num_picked
        @fact test --> expected
    end
    context("modifying the value in the second parent from VirtualArray") do
        # set up
        num_dims = 2
        len = rand(2:10) # no larger than 5
        expanded_dim = 2
        num_picked = rand(1:1000)
        index_picked = rand(1:len)

        a = rand(len)
        b = rand(len)
        parents = Array[a, b]

        expected = cat(expanded_dim, parents...)
        test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

        test[index_picked, 2] = num_picked
        expected[index_picked, 2] = num_picked

        @fact test[index_picked, 2] --> num_picked
        @fact b[index_picked] --> num_picked
        @fact test --> expected
    end
    context("modifying the value in the first parent from the parent") do
        # set up
        num_dims = 2
        len = rand(2:10) # no larger than 5
        expanded_dim = 2
        num_picked = rand(1:1000)
        index_picked = rand(1:len)

        a = rand(len)
        b = rand(len)
        parents = Array[a, b]

        expected = cat(expanded_dim, parents...)
        test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

        a[index_picked] = num_picked
        expected = cat(expanded_dim, parents...)

        @fact test[index_picked, 1] --> num_picked
        @fact a[index_picked] --> num_picked
        @fact test --> expected
    end
    context("modifying the value in the second parent from the parent") do
        # set up
        num_dims = 2
        len = rand(2:10) # no larger than 5
        expanded_dim = 2
        num_picked = rand(1:1000)
        index_picked = rand(1:len)

        a = rand(len)
        b = rand(len)
        parents = Array[a, b]

        test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

        b[index_picked] = num_picked
        expected = cat(expanded_dim, parents...)

        @fact test[index_picked, 2] --> num_picked
        @fact b[index_picked] --> num_picked
        @fact test --> expected
    end
end

############################################################################################
# ERRORS
############################################################################################

facts("Errors while using VirtualArray") do
    context("out of bounds indexing on 1 1d array") do
        len = rand(1:100)
        a = rand(len)
        test = VirtualArray{Float64, 1}(a)
        @fact_throws BoundsError test[-1]
        @fact_throws BoundsError test[0]
        @fact_throws BoundsError test[len+1]
    end
    context("tring to create a VirtualArray of the wrong type") do
        len = rand(1:100)
        a = rand(len)
        @fact_throws MethodError test = VirtualArray{AbstractString, 1}(a)
        @fact_throws InexactError test = VirtualArray{Int, 1}(a)
    end
    context("out of bounds indexing on multiple 1d arrays") do
        # set up
        num_parents = rand(3:10)
        len = rand(1:100)

        parents = []
        for i in 1:num_parents
            push!(parents, rand(len))
        end
        test = VirtualArray{Float64, 1}(parents...)
        @fact_throws BoundsError test[-1]
        @fact_throws BoundsError test[len*(num_parents)+1]
    end
    context("out of bounds setting on 1 1d array") do
        len = rand(1:100)
        a = rand(len)
        test = VirtualArray{Float64, 1}(a)
        @fact_throws BoundsError test[-1] = 1
        @fact_throws BoundsError test[len+1] = 1
    end
    context("out of bounds setting on multiple 1d arrays") do
        # set up
        num_parents = rand(3:10)
        len = rand(1:100)

        parents = []
        for i in 1:num_parents
            push!(parents, rand(len))
        end
        test = VirtualArray{Float64, 1}(parents...)
        @fact_throws BoundsError test[-1] = 1
        @fact_throws BoundsError test[len*(num_parents)+1] = 1
    end
    context("throw the right type of error when accessing the array wrong") do

        # set up
        a = rand(3,3,3)
        b = rand(3,3,3)

        c = cat(1,a,b)
        test = VirtualArray{Float64, 3}(1, a, b)

        @fact_throws BoundsError c[6,4,1]
        @fact_throws BoundsError test[6,4,1]

        c_error = 1
        try
            c[6,4,1]
        catch e
            c_error = e
        end

        test_error = 1
        try
            test[6,4,1]
        catch e
            test_error = e
        end

        @fact test_error.a --> c_error.a
        @fact test_error.i --> c_error.i

    end
    context("trying to make an 2 d array of 1 d parents of different sizes") do
        # set up
        num_dims = 2
        len = rand(2:2:20) # no larger than 5
        len_2 = rand(3:2:30) # no larger than 5
        expanded_dim = 2

        a = rand(len)
        b = rand(len_2)
        parents = Array[a, b]

        @fact_throws DimensionMismatch expected = cat(expanded_dim, parents...)
        @fact_throws DimensionMismatch test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

        c_error = 1
        try
            expected = cat(expanded_dim, parents...)
        catch e
            c_error = e
        end

        test_error = 1
        try
            test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)
        catch e
            test_error = e
        end

        @fact contains(c_error.msg, test_error.msg) --> true
    end
    context("creating an 2 d VirtualArray with 2 d parents of different lengths on the wrong dimension") do
        # set up
        num_dims = 2
        len = rand(2:2:20) # no larger than 5
        len_2 = rand(3:2:21) # no larger than 5
        expanded_dim = 1 # to work, this should be 2

        a = rand(len, len)
        b = rand(len, len_2)
        parents = Array[a, b]

        @fact_throws DimensionMismatch expected = cat(expanded_dim, parents...)
        @fact_throws DimensionMismatch test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)
    end
end

############################################################################################
# MEMORY ALLOCATION
############################################################################################

facts("Check the memory usage of VirtualArray") do
    context("simple test") do
        a = collect(1:1)

        # Make sure everything is compiled before testing
        VirtualArray{Int, 1}(1,a)

        expected = @allocated VirtualArray{Int, 1}(1,a)

        test = @allocated VirtualArray{Int, 1}(1,a)
        @fact test --> expected

        a = collect(1:10)
        test = @allocated VirtualArray{Int, 1}(1,a)
        @fact test --> expected

        a = collect(1:10000)
        test = @allocated VirtualArray{Int, 1}(1,a)
        @fact test --> expected

        a = collect(1:10000000)
        test = @allocated VirtualArray{Int, 1}(1,a)
        @fact test --> expected
    end
end

@testset "VirtualArrays Tests" begin
    @testset "range getindex" begin
        @testset "getting a range value in one 1 d parent" begin
            # set up
            num = rand(1:1000)
            num_parents = rand(3:10)
            len = rand(1:100)
            change_p = rand(1:num_parents)
            change_start = rand(1:len)
            change_end = rand(change_start:len)
            change_range = change_start:change_end

            parents = []
            for i in 1:num_parents
                push!(parents, collect(num:num+len-1))
            end

            expected = cat(1, parents...)
            test = VirtualArray{Int, 1}(parents...)

            @test test.parents[change_p][change_range] == test[(change_p-1)*len+change_range]
            @test test[(change_p-1)*len+change_range] == expected[(change_p-1)*len+change_range]
            @test test == expected
        end
        @testset "getting a range value in multi 1 d parent" begin
            # set up
            num = rand(1:1000)
            num_parents = rand(3:10)
            len = rand(1:100)
            change_p_start = rand(1:num_parents-1)
            change_p_end = rand(change_p_start+1:num_parents)
            change_start = rand(1:len)
            change_end = rand(1:len)
            change_range = (change_p_start-1)*len+change_start:(change_p_end-1)*len+change_end

            parents = []
            for i in 1:num_parents
                push!(parents, collect(num:num+len-1))
            end

            expected = cat(1, parents...)
            test = VirtualArray{Int, 1}(parents...)

            @test test[change_range] == expected[change_range]
            @test test == expected
        end
        @testset "getting a range value in multi multi d parent" begin

            # set up
            # keep these numbers small because we can run out of memory or get very slow tests
            num_parents = rand(3:10) # no larger than 10
            num_dims = rand(3:6) # no larger than 6
            len = rand(2:5) # no larger than 5
            expanded_dim = rand(1:num_dims)

            dims = zeros(Int, num_dims) + len

            parents = []
            for i in 1:num_parents
                push!(parents, rand(dims...))
            end
            change_to = rand(1:10)
            change_p_start = rand(1:num_parents-1)
            change_p_end = rand(change_p_start+1:num_parents)
            change_start = rand(1:len)
            change_end = rand(1:len)
            change_range_expanded_dim = (change_p_start-1)*len+change_start:(change_p_end-1)*len+change_end

            change_range = []
            for i in 1:num_dims
                if i != expanded_dim
                    change_start = rand(1:len)
                    change_end = rand(change_start:len)
                    push!(change_range, change_start:change_end)
                else
                    push!(change_range, change_range_expanded_dim)
                end
            end

            expected = cat(expanded_dim, parents...)
            test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

            @test test[change_range...] == expected[change_range...]
            @test test == expected

        end
        @testset "getting a range value in multi M d parent like an N d range" begin

            # set up
            # keep these numbers small because we can run out of memory or get very slow tests
            num_parents = rand(3:10) # no larger than 10
            num_dims = rand(3:6) # no larger than 6
            len = rand(2:5) # no larger than 5
            expanded_dim = rand(1:num_dims)

            dims = zeros(Int, num_dims) + len

            parents = []
            for i in 1:num_parents
                push!(parents, rand(dims...))
            end
            change_to = rand(1:10)
            change_p_start = rand(1:num_parents-1)
            change_p_end = rand(change_p_start+1:num_parents)
            change_start = rand(1:len)
            change_end = rand(1:len)
            change_range_expanded_dim = (change_p_start-1)*len+change_start:(change_p_end-1)*len+change_end

            range_dim = rand(1:num_dims-1) # there's a chance we don't have an index on the expanded_dim

            change_range = []
            for i in 1:range_dim-1
                if i != expanded_dim
                    change_start = rand(1:len)
                    change_end = rand(change_start:len)
                    push!(change_range, change_start:change_end)
                else
                    push!(change_range, change_range_expanded_dim)
                end
            end

            expected = cat(expanded_dim, parents...)
            test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

            end_left = 1
            for i in range_dim:num_dims
                end_left *= size(test)[i]
            end

            push!(change_range, 1:end_left)

            @test test[change_range...] == expected[change_range...]
            @test test == expected

        end
        @testset "getting a range value in multi M d parent like an N d range where N is < expanded dim and expanded dim < M" begin

            # set up
            # keep these numbers small because we can run out of memory or get very slow tests
            num_parents = rand(3:10) # no larger than 10
            num_dims = rand(3:6) # no larger than 6
            len = rand(2:5) # no larger than 5
            expanded_dim = rand(2:num_dims-1)

            dims = zeros(Int, num_dims) + len

            parents = []
            for i in 1:num_parents
                push!(parents, rand(dims...))
            end
            change_to = rand(1:10)
            change_p_start = rand(1:num_parents-1)
            change_p_end = rand(change_p_start+1:num_parents)
            change_start = rand(1:len)
            change_end = rand(1:len)
            change_range_expanded_dim = (change_p_start-1)*len+change_start:(change_p_end-1)*len+change_end

            range_dim = rand(1:expanded_dim-1)

            change_range = []
            for i in 1:range_dim-1
                if i != expanded_dim
                    change_start = rand(1:len)
                    change_end = rand(change_start:len)
                    push!(change_range, change_start:change_end)
                else
                    push!(change_range, change_range_expanded_dim)
                end
            end

            expected = cat(expanded_dim, parents...)
            test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

            end_left = 1
            for i in range_dim:num_dims
                end_left *= size(test)[i]
            end

            push!(change_range, 1:end_left)

            @test test[change_range...] == expected[change_range...]
            @test test == expected

        end
        @testset "getting a range value in multi M d parent like an N d range where N is < expanded dim, expanded dim < M, and length of last dim > number of dim" begin

            # set up
            # keep these numbers small because we can run out of memory or get very slow tests
            num_parents = rand(3:10) # no larger than 10
            num_dims = rand(4:6) # no larger than 6
            expanded_dim = rand(2:num_dims-2)
            len = rand(expanded_dim+1:5) # no larger than 5

            dims = zeros(Int, num_dims) + len

            parents = []
            for i in 1:num_parents
                push!(parents, rand(dims...))
            end
            change_to = rand(1:10)
            change_p_start = rand(1:num_parents-1)
            change_p_end = rand(change_p_start+1:num_parents)
            change_start = rand(1:len)
            change_end = rand(1:len)
            change_range_expanded_dim = (change_p_start-1)*len+change_start:(change_p_end-1)*len+change_end

            range_dim = rand(1:expanded_dim-1)

            change_range = []
            for i in 1:range_dim-1
                if i != expanded_dim
                    change_start = rand(1:len)
                    change_end = rand(change_start:len)
                    push!(change_range, change_start:change_end)
                else
                    push!(change_range, change_range_expanded_dim)
                end
            end

            expected = cat(expanded_dim, parents...)
            test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

            end_left = 1
            for i in range_dim:num_dims
                end_left *= size(test)[i]
            end

            push!(change_range, 1:end_left)

            @test test[change_range...] == expected[change_range...]
            @test test == expected

        end
        @testset "getting a range value in multi M d parent like an N d range where N is < expanded dim and expanded dim = M" begin

            # set up
            # keep these numbers small because we can run out of memory or get very slow tests
            num_parents = rand(3:10) # no larger than 10
            num_dims = rand(3:6) # no larger than 6
            len = rand(2:5) # no larger than 5
            expanded_dim = num_dims

            dims = zeros(Int, num_dims) + len

            parents = []
            for i in 1:num_parents
                push!(parents, rand(dims...))
            end
            change_to = rand(1:10)
            change_p_start = rand(1:num_parents-1)
            change_p_end = rand(change_p_start+1:num_parents)
            change_start = rand(1:len)
            change_end = rand(1:len)
            change_range_expanded_dim = (change_p_start-1)*len+change_start:(change_p_end-1)*len+change_end

            range_dim = rand(1:expanded_dim-1)

            change_range = []
            for i in 1:range_dim-1
                if i != expanded_dim
                    change_start = rand(1:len)
                    change_end = rand(change_start:len)
                    push!(change_range, change_start:change_end)
                else
                    push!(change_range, change_range_expanded_dim)
                end
            end

            expected = cat(expanded_dim, parents...)
            test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

            end_left = 1
            for i in range_dim:num_dims
                end_left *= size(test)[i]
            end

            push!(change_range, 1:end_left)

            @test test[change_range...] == expected[change_range...]
            @test test == expected

        end
        @testset "getting a range value in multi M d parent like an N d range where N = expanded dim" begin

            # set up
            # keep these numbers small because we can run out of memory or get very slow tests
            num_parents = rand(3:10) # no larger than 10
            num_dims = rand(3:6) # no larger than 6
            len = rand(2:5) # no larger than 5
            expanded_dim = rand(1:num_dims-1)

            dims = zeros(Int, num_dims) + len

            parents = []
            for i in 1:num_parents
                push!(parents, rand(dims...))
            end
            change_to = rand(1:10)
            change_p_start = rand(1:num_parents-1)
            change_p_end = rand(change_p_start+1:num_parents)
            change_start = rand(1:len)
            change_end = rand(1:len)
            change_range_expanded_dim = (change_p_start-1)*len+change_start:(change_p_end-1)*len+change_end

            range_dim = expanded_dim

            change_range = []
            for i in 1:range_dim-1
                if i != expanded_dim
                    change_start = rand(1:len)
                    change_end = rand(change_start:len)
                    push!(change_range, change_start:change_end)
                else
                    push!(change_range, change_range_expanded_dim)
                end
            end

            expected = cat(expanded_dim, parents...)
            test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

            end_left = 1
            for i in range_dim:num_dims
                end_left *= size(test)[i]
            end

            push!(change_range, 1:end_left)

            @test test[change_range...] == expected[change_range...]
            @test test == expected

        end
        @testset "getting a range value in multi M d parent like an N d range where N > expanded dim" begin

            # set up
            # keep these numbers small because we can run out of memory or get very slow tests
            num_parents = rand(3:10) # no larger than 10
            num_dims = rand(3:6) # no larger than 6
            len = rand(2:5) # no larger than 5
            expanded_dim = rand(1:num_dims-2)

            dims = zeros(Int, num_dims) + len

            parents = []
            for i in 1:num_parents
                push!(parents, rand(dims...))
            end
            change_to = rand(1:10)
            change_p_start = rand(1:num_parents-1)
            change_p_end = rand(change_p_start+1:num_parents)
            change_start = rand(1:len)
            change_end = rand(1:len)
            change_range_expanded_dim = (change_p_start-1)*len+change_start:(change_p_end-1)*len+change_end

            range_dim = rand(expanded_dim+1:num_dims-1) # there's a chance we don't have an index on the expanded_dim

            change_range = []
            for i in 1:range_dim-1
                if i != expanded_dim
                    change_start = rand(1:len)
                    change_end = rand(change_start:len)
                    push!(change_range, change_start:change_end)
                else
                    push!(change_range, change_range_expanded_dim)
                end
            end

            expected = cat(expanded_dim, parents...)
            test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

            end_left = 1
            for i in range_dim:num_dims
                end_left *= size(test)[i]
            end

            push!(change_range, 1:end_left)

            @test test[change_range...] == expected[change_range...]
            @test test == expected

        end
        @testset "getting a range value with nothing in the square brackets" begin
            # set up
            # keep these numbers small because we can run out of memory or get very slow tests
            num_parents = rand(3:10) # no larger than 10
            num_dims = rand(3:6) # no larger than 6
            len = rand(2:5) # no larger than 5
            expanded_dim = rand(1:num_dims)

            dims = zeros(Int, num_dims) + len

            parents = []
            for i in 1:num_parents
                push!(parents, rand(dims...))
            end

            expected = cat(expanded_dim, parents...)
            test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

            @test test[] == expected[]
            @test test[] == test[1]
            @test test == expected
        end
    end
    @testset "range setindex" begin
        @testset "setting a range value in one 1 d parent" begin
            # set up
            num = rand(1:1000)
            num_parents = rand(3:10)
            len = rand(1:100)
            change_p = rand(1:num_parents)
            change_start = rand(1:len)
            change_end = rand(change_start:len)
            change_range = change_start:change_end
            change_to = rand(1:10)

            parents = []
            for i in 1:num_parents
                push!(parents, collect(num:num+len-1))
            end

            test = VirtualArray{Int, 1}(parents...)
            test[change_range] = change_to
            expected = cat(1, parents...)

            @test test.parents[change_p][change_range] == test[(change_p-1)*len+change_range]
            @test test[(change_p-1)*len+change_range] == expected[(change_p-1)*len+change_range]
            @test test == expected
        end
        @testset "setting a range value in multi 1 d parent" begin
            # set up
            num = rand(1:1000)
            num_parents = rand(3:10)
            len = rand(1:100)
            change_p_start = rand(1:num_parents-1)
            change_p_end = rand(change_p_start+1:num_parents)
            change_start = rand(1:len)
            change_end = rand(1:len)
            change_range = (change_p_start-1)*len+change_start:(change_p_end-1)*len+change_end
            change_to = rand(1:10)

            parents = []
            for i in 1:num_parents
                push!(parents, collect(num:num+len-1))
            end

            test = VirtualArray{Int, 1}(parents...)
            test[change_range] = change_to
            expected = cat(1, parents...)

            @test test[change_range] == expected[change_range]
            @test test == expected
        end
        @testset "setting a range value in multi multi d parent" begin

            # set up
            # keep these numbers small because we can run out of memory or get very slow tests
            num_parents = rand(3:10) # no larger than 10
            num_dims = rand(3:6) # no larger than 6
            len = rand(2:5) # no larger than 5
            expanded_dim = rand(1:num_dims)

            dims = zeros(Int, num_dims) + len

            parents = []
            for i in 1:num_parents
                push!(parents, rand(dims...))
            end
            change_to = rand(1:10)
            change_p_start = rand(1:num_parents-1)
            change_p_end = rand(change_p_start+1:num_parents)
            change_start = rand(1:len)
            change_end = rand(1:len)
            change_range_expanded_dim = (change_p_start-1)*len+change_start:(change_p_end-1)*len+change_end

            change_range = []
            for i in 1:num_dims
                if i != expanded_dim
                    change_start = rand(1:len)
                    change_end = rand(change_start:len)
                    push!(change_range, change_start:change_end)
                else
                    push!(change_range, change_range_expanded_dim)
                end
            end
            change_to = rand(1:10)

            test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)
            test[change_range...] = change_to
            expected = cat(expanded_dim, parents...)

            @test test[change_range...] == expected[change_range...]
            @test test == expected

        end
        @testset "setting a range value in multi M d parent like an N d range " begin

            # set up
            # keep these numbers small because we can run out of memory or get very slow tests
            num_parents = rand(3:10) # no larger than 10
            num_dims = rand(3:6) # no larger than 6
            len = rand(2:5) # no larger than 5
            expanded_dim = rand(1:num_dims)

            dims = zeros(Int, num_dims) + len

            parents = []
            for i in 1:num_parents
                push!(parents, rand(dims...))
            end
            change_to = rand(1:10)
            change_p_start = rand(1:num_parents-1)
            change_p_end = rand(change_p_start+1:num_parents)
            change_start = rand(1:len)
            change_end = rand(1:len)
            change_range_expanded_dim = (change_p_start-1)*len+change_start:(change_p_end-1)*len+change_end

            range_dim = rand(1:num_dims-1) # there's a chance we don't have an index on the expanded_dim

            change_range = []
            for i in 1:range_dim-1
                if i != expanded_dim
                    change_start = rand(1:len)
                    change_end = rand(change_start:len)
                    push!(change_range, change_start:change_end)
                else
                    push!(change_range, change_range_expanded_dim)
                end
            end

            expected = cat(expanded_dim, parents...)
            test = VirtualArray{Float64, num_dims}(expanded_dim, parents...)

            end_left = 1
            for i in range_dim:num_dims
                end_left *= size(test)[i]
            end
            push!(change_range, 1:end_left)

            change_to = rand(1:10)
            test[change_range...] = change_to
            expected = cat(expanded_dim, parents...)

            @test test[change_range...] == expected[change_range...]
            @test test == expected

        end
    end
end


text:
#=
predictors.jl; This file implements methods to compare encoded strings to references and to make predictions from this comparison.
=#

### Cosine distance predictor
"""
This function makes predictions for the given encodings by finding the closest class encoding in terms of cosine distance.\n
Input:
    - a vector of encoded test sequences.\n
    - a dictionary of encoded reference sequences with their class a key.\n
Output: a vector of class predictions in the same order as the input test encodings.\n

DISCLAIMER: Currently every test encoding given to this function will be classified to the closest class, even when the true class is not in the train encodings.
=> ALLE KLASSES EVEN VER VAN EXAMPLE => UNCLASSIFIED
"""
function cosine_predict(test_encodings::Vector, train_encodings::Dict)
    predictions = Vector{String}()
    for hv in test_encodings
        closest_distance = Inf
        match = missing
        for (k, v) in train_encodings
            dist = cosine_dist(hv, v)
            if (dist < closest_distance)
                closest_distance = dist
                match = k
            end
        end
        push!(predictions, match)
    end
    return predictions
end

### Naive bayes predictor

ext: jl


text:
# ------------------------------- #
#  -- Logistic Regression (L2) -- #
# ------------------------------- #
# ->  f(θ)  = -∑logσ(yXθ) + λ|θ|₂²/2
# -> ∇f(θ)  =  X'(y(w-1)) + λθ
# -> ∇²f(θ) =  X' Diag(w(1-w)) X + λI
# NOTE:
# * w = σ(yXθ)
# * yᵢ ∈ {±1} so that y² = 1
# * -σ(-x) == (σ(x)-1)
# NOTE: https://github.com/JuliaAI/MLJLinearModels.jl/issues/104
# ---------------------------------------------------------

function fgh!(glr::GLR{LogisticLoss,<:L2R}, X, y, scratch)
    J    = objective(glr) # GLR objective (loss+penalty)
    n, p = size(X)
    λ    = getscale(glr.penalty)
    if glr.fit_intercept
        (f, g, H, θ) -> begin
            Xθ = scratch.n
            apply_X!(Xθ, X, θ)                       # -- Xθ = apply_X(X, θ)
            # precompute σ(yXθ)
            w    = scratch.n2
            w   .= σ.(Xθ .* y)                       # -- w  = σ.(Xθ .* y)
            g === nothing || begin
                t  = scratch.n3
                t .= y .* w .- y                     # -- t = y .* (w .- 1.0)
                apply_Xt!(g, X, t)                   # -- g = X't
                g .+= λ .* θ
                glr.penalize_intercept || (g[end] -= λ * θ[end])
            end
            H === nothing || begin
                # NOTE: we could try to be clever to reduce the allocations for
                # ΛX but computing the full hessian allocates a lot anyway so
                # probably not really worth it
                t    = scratch.n3
                t   .= w .- w.^2                      # σ(yXθ)(1-σ(yXθ))
                a    = 1:p
                ΛX   = t .* X                         # !! big allocs
                mul!(view(H, a, a), X', ΛX)           # -- H[1:p,1:p] = X'ΛX
                ΛXt1   = view(scratch.p, a)
                ΛXt1 .*= 0
                @inbounds for i in a, j in 1:n
                    ΛXt1[i] += ΛX[j, i]               # -- (ΛX)'1
                end
                @inbounds for i in a
                    H[i, end] = H[end, i] = ΛXt1[i]   # -- H[:,p+1] = (ΛX)'1
                end
                H[end, end] = sum(t)                  # -- 1'Λ1'
                add_λI!(H, λ, glr.penalize_intercept) # -- H = -X'ΛX + λI
            end
            f === nothing || return J(y, Xθ, view_θ(glr, θ))
        end
    else
        # see comments above, same computations just no additional things for
        # fit_intercept
        (f, g, H, θ) -> begin
            Xθ = scratch.n
            apply_X!(Xθ, X, θ)
            w    = scratch.n2
            w   .= σ.(y .* Xθ)
            g === nothing || begin
                t  = scratch.n3
                t .= y .* w .- y
                apply_Xt!(g, X, t)
                g .+= λ .* θ
            end
            H === nothing || begin
                t  = scratch.n3
                t .= w .- w.^2
                mul!(H, X', t .* X)
                add_λI!(H, λ)
            end
            f === nothing || return J(y, Xθ, θ)
        end
    end
end

function Hv!(glr::GLR{LogisticLoss,<:L2R}, X, y, scratch)
    n, p = size(X)
    λ    = getscale(glr.penalty)
    if glr.fit_intercept
        # H = [X 1]'Λ[X 1] + λ I
        # rows a 1:p = [X'ΛX + λI | X'Λ1]
        # row  e end = [1'ΛX      | sum(a)+λ]
        (Hv, θ, v) -> begin
            Xθ = scratch.n
            apply_X!(Xθ, X, θ)                       # -- Xθ = apply_X(X, θ)
            w   = scratch.n2
            w  .= σ.(Xθ .* y)                        # -- w  = σ.(Xθ .* y)
            w .-= w.^2                               # -- w  = w(1-w)
            # view on the first p rows
            a    = 1:p
            Hvₐ  = view(Hv, a)
            vₐ   = view(v,  a)
            XtΛ1 = view(scratch.p, 1:p)
            mul!(XtΛ1, X', w)                        # -- X'Λ1; O(np)
            vₑ   = v[end]
            # update for the first p rows -- (X'X + λI)v[1:p] + (X'1)v[end]
            Xvₐ  = scratch.n
            mul!(Xvₐ, X, vₐ)
            Xvₐ .*=  w                               # --  ΛXvₐ
            mul!(Hvₐ, X', Xvₐ)                       # -- (X'ΛX)vₐ
            Hvₐ .+= λ .* vₐ .+ XtΛ1 .* vₑ
            # update for the last row -- (X'1)'v + n v[end]
            Hv[end] = dot(XtΛ1, vₐ) +
                        (sum(w) + λ_if_penalize_intercept(glr, λ)) * vₑ
        end
    else
        (Hv, θ, v) -> begin
            Xθ = scratch.n
            apply_X!(Xθ, X, θ)
            w   = scratch.n2
            w  .= σ.(Xθ .* y)                # -- σ(yXθ)
            w .-= w.^2
            Xv  = scratch.n
            mul!(Xv, X, v)
            Xv .*= scratch.n2               # -- ΛXv
            mul!(Hv, X', Xv)                # -- X'ΛXv
            Hv .+= λ .* v
        end
    end
end

# ----------------------------------- #
#  -- L1/Elnet Logistic Regression -- #
# ----------------------------------- #
# ->  J(θ)  = f(θ) + r(θ)
# ->  f(θ)  = LL + λ|θ|₂²  // smooth (LL = LogisticLoss)
# ->  r(θ)  = γ|θ|₁        // non-smooth with prox
# -> ∇f(θ)  = ∇LL + λθ
# -> ∇²f(θ) = ∇²LL + λI
# -> prox_r = soft-thresh
# ---------------------------------------------------------

function smooth_fg!(glr::GLR{LogisticLoss,<:ENR}, X, y, scratch)
    smooth = get_smooth(glr)
    (g, θ) -> fgh!(smooth, X, y, scratch)(0.0, g, nothing, θ)
end

# ---------------------------------- #
#  -- Multinomial Regression (L2) -- #
# ---------------------------------- #
# ->  c is the number of classes, θ has dims p * c
# ->  P = X * θ
# -> Zᵢ = ∑ exp(Pᵢ)
# -> Λ  = Diagonal(-Z)
# ->  f(θ)   = ∑(log Zᵢ - P[i, y[i]]) +  λ|θ|₂²
# -> ∇f(θ)   = reshape(X'ΛM, c * p)
# -> ∇²f(θ)v = via R operator
# NOTE:
# * yᵢ ∈ {1, 2, ..., c}
# ---------------------------------------------------------

function fg!(glr::GLR{<:MultinomialLoss,<:L2R}, X, y, scratch)
    n, p = size(X)
    c    = getc(glr, y)
    λ    = getscale(glr.penalty)
    (f, g, θ) -> begin
        P  = scratch.nc
        apply_X!(P, X, θ, c, scratch)                # O(npc) store n * c
        M  = scratch.nc2
        M .= exp.(P)                                 # O(npc) store n * c
        g === nothing || begin
            ΛM  = scratch.nc3
            ΛM .= M ./ sum(M, dims=2)                # O(nc)  store n * c
            Q   = scratch.nc4
            @inbounds for i = 1:n, j=1:c
                Q[i, j] = ifelse(y[i] == j, 1.0, 0.0)
            end
            ∑ΛM = sum(ΛM, dims=1)
            ∑Q  = sum(Q, dims=1)
            R   = ΛM
            R .-= Q
            G   = scratch.pc
            if glr.fit_intercept
                mul!(view(G, 1:p, :), X', R)
                @inbounds for k in 1:c
                    G[end, k] = ∑ΛM[k] - ∑Q[k]
                end
            else
                mul!(G, X', R)
            end
            g  .= reshape(G, (p + Int(glr.fit_intercept)) * c)
            g .+= λ .* θ
            glr.fit_intercept &&
                (glr.penalize_intercept || (g[end] -= λ * θ[end]))
        end
        f === nothing || begin
            # we re-use pre-computations here, see also MultinomialLoss
            # ms = maximum(P, dims=2)
            # ss = sum(M ./ exp.(ms), dims=2)
            ms   = maximum(P, dims=2)
            ems  = scratch.n
            @inbounds for i in 1:n
                ems[i] = exp(ms[i])
            end
            ΛM  = scratch.nc2  # note that _NC is already linked to P
            ΛM .= M ./ ems
            ss  = sum(ΛM, dims=2)
            t   = 0.0
            @inbounds for i in eachindex(y)
                t += log(ss[i]) + ms[i] - P[i, y[i]]
            end
            return sum(t) + glr.penalty(view_θ(glr, θ))
        end
    end
end

function Hv!(glr::GLR{<:MultinomialLoss,<:L2R}, X, y, scratch)
    p = size(X, 2)
    λ = getscale(glr.penalty)
    c = getc(glr, y)
    # NOTE:
    # * ideally P and Q should be recuperated from gradient computations (fghv!)
    # * assumption that c is small so that storing matrices of size n * c is
    # not too bad; if c is large and allocations should be minimized, all these
    # computations can be done per class with views over (c-1)p+1:cp; it will
    # allocate less but is likely slower; maybe in the future we could have a
    # keyword indicating which one the user wants to use.
    (Hv, θ, v) -> begin
        P  = apply_X(X, θ, c)     # P_ik = <x_i, θ_k> // dims n * c; O(npc)
        Q  = apply_X(X, v, c)     # Q_ik = <x_i, v_k>    // dims n * c; O(npc)
        M  = exp.(P)              # M_ik = exp<x_i, w_k> // dims n * c;
        MQ = M .* Q               #                      // dims n * c; O(nc)
        ρ  = 1 ./ sum(M, dims=2)  # ρ_i = 1/Z_i = 1/∑_k exp<x_i, w_k>
        κ  = sum(MQ, dims=2)      # κ_i  = ∑_k exp<x_i, w_k><x_i, v_k>
        γ  = κ .* ρ.^2            # γ_i  = κ_i / Z_i^2
        # computation of Hv
        U      = (ρ .* MQ) .- (γ .* M)                  # // dims n * c; O(nc)
        Hv_mat = X' * U                                 # // dims n * c; O(npc)
        if glr.fit_intercept
            Hv .= reshape(vcat(Hv_mat, sum(U, dims=1)), (p+1)*c)
        else
            Hv .= reshape(Hv_mat, p * c)
        end
        Hv .+= λ .* v
        glr.fit_intercept && (glr.penalize_intercept || (Hv[end] -= λ * v[end]))
    end
end

# -------------------------------------- #
#  -- L1/Elnet Multinomial Regression -- #
# -------------------------------------- #
# ->  J(θ)  = f(θ) + r(θ)
# ->  f(θ)  = MN + λ|θ|₂²  // smooth (MN = MultinomialLoss)
# ->  r(θ)  = γ|θ|₁        // non-smooth with prox
# -> ∇f(θ)  = ∇MN + λθ
# -> ∇²f(θ) = ∇²MN + λI
# -> prox_r = soft-thresh
# ---------------------------------------------------------

function smooth_fg!(glr::GLR{<:MultinomialLoss,<:ENR}, X, y, scratch)
    smooth = get_smooth(glr)
    (g, θ) -> fg!(smooth, X, y, scratch)(0.0, g, θ)
end


text:
immutable ObsAdaptiveSRFilter{S} <: Updater
    pomdp::POMDP{S}
    resample::Any
    max_frac_replaced::Float64
    rng::AbstractRNG
end

function POMDPs.initialize_belief{S}(up::ObsAdaptiveSRFilter{S}, d::Any)
    sc = resample(up.resample, d, up.rng)
    return ParticleCollection{Tuple{S,Float64}}([(p,0.0) for p in particles(sc)])
end
POMDPs.update(up::ObsAdaptiveSRFilter, b, a, o) = update(up, initialize_belief(up, b), a, o)

function POMDPs.update{S}(up::ObsAdaptiveSRFilter{S}, b::ParticleFilters.ParticleCollection{Tuple{S,Float64}}, a, o)
    if n_particles(b) > 2*up.resample.n
        b = resample(up.resample, b, up.rng)
    end

    ps = particles(b)
    pm = Array{Tuple{S,Float64}}(0)
    wm = Array{Float64}(0)
    sizehint!(pm, n_particles(b))
    sizehint!(wm, n_particles(b))
    all_terminal = true
    for i in 1:n_particles(b)
        sr = ps[i]
        s = first(sr)
        if !isterminal(up.pomdp, s)
            all_terminal = false
            sp, r = generate_sr(up.pomdp, s, a, up.rng)
            push!(pm, (sp,r))
            od = observation(up.pomdp, s, a, sp)
            push!(wm, pdf(od, o))
        end
    end
    ws = sum(wm)
    if all_terminal || sum(wm) == 0.0
        # warn("All states in particle collection were terminal.")
        return initialize_belief(up, reset_distribution(up.pomdp, a, o))
    end

    pc = resample(up.resample, WeightedParticleBelief{Tuple{S,Float64}}(pm, wm, ws, nothing), up.rng)
    ps = particles(pc)
    # for i in 1:length(ps)
    #     ps[i] += 0.001*randn(up.rng, 2)
    # end

    mpw = max_possible_weight(up.pomdp, a, o)
    frac_replaced = up.max_frac_replaced*max(0.0, 1.0 - maximum(wm)/mpw)
    n_replaced = floor(Int, frac_replaced*length(ps))
    is = randperm(up.rng, length(ps))[1:n_replaced]
    for i in is
        ps[i] = (new_particle(up.pomdp, a, o, up.rng), 0.0)
    end
    return pc
end


text:
# Length
# ------

resolve(box::AbsoluteBox, x::Length{:mm}) = x
resolve(box::AbsoluteBox, x::Length{:w})  = width(box).value * x.value * mm
resolve(box::AbsoluteBox, x::Length{:h})  = height(box).value * x.value * mm
resolve(box::AbsoluteBox, x::Length{:d})  = depth(box).value * x.value * mm

# Operations
# ----------

resolve(box::AbsoluteBox, x::Neg) = -resolve(box, x.a)
resolve(box::AbsoluteBox, x::Add) = resolve(box, x.a) + resolve(box, x.b)
resolve(box::AbsoluteBox, x::Mul) = resolve(box, x.a) * x.b
resolve(box::AbsoluteBox, x::Div) = resolve(box, x.a) / x.b
resolve(box::AbsoluteBox, x::Min) = min(resolve(box, x.a), resolve(box, x.b))
resolve(box::AbsoluteBox, x::Max) = max(resolve(box, x.a), resolve(box, x.b))
resolve(box::AbsoluteBox, p::Vec) = map(x -> resolve(box, x), p) + box.x0
resolve(outer::AbsoluteBox, box::BoundingBox) =
    BoundingBox(resolve(outer, box.x0), map(x -> resolve(outer, x), box.a))


text:
# please instantiate ./test/Manifest.toml for this

using JLD2
import Statistics

# ========================================
# Asymmetric Fuzz Test Fixtures
# ========================================

# format: [(values, window_size, expected_medians)]

fixtures = Vector{Tuple{Vector{Float64}, Integer, Vector{Float64}}}(undef, 0)

function naive_asymmetric_median(input, window_size)
    if window_size > length(input)
        window_size = length(input)
    end
    growing_phase_inds = [1:k for k in 1:window_size]
    rolling_phase_inds = [k:k + window_size - 1 for k in 2:(length(input) - window_size + 1)]
    shrinking_phase_inds = [length(input) - k + 1:length(input) for k in window_size - 1:-1:1]
    phase_inds = [growing_phase_inds; rolling_phase_inds; shrinking_phase_inds]
    output = [Statistics.median(input[inds]) for inds in phase_inds]
    return output
end

for i in 1:100
    N = rand(1:50)
    window_size = rand(1:60)
    values = rand(N)
    expected_medians = naive_asymmetric_median(values, window_size)
    push!(fixtures, (values, window_size, expected_medians))
end

filename = dirname(@__FILE__) * "/asymmetric.jld2"
@save filename fixtures

# Again with full Int32 range

fixtures = Vector{Tuple{Vector{Int32}, Integer, Vector{Float64}}}(undef, 0)

for i in 1:100
    N = rand(1:50)
    window_size = rand(1:60)
    values = rand(Int32, N)
    expected_medians = naive_asymmetric_median(values, window_size)
    push!(fixtures, (values, window_size, expected_medians))
end

filename = dirname(@__FILE__) * "/asymmetric_int.jld2"
@save filename fixtures

max_stars_repo_name: Firionus/FastRunningMedian

max_issues_repo_name: Firionus/FastRunningMedian


text:
type Dimension
	name::Symbol # the column name in dataframe
	accessor::ASCIIString #var paymentsByTotal = payments.dimension(function(d) { return d.total; });
	bin_width::Float64 # discretization width, NaN if unused

	Dimension(name::Symbol, accessor::ASCIIString, bin_width::Float64=NaN) = new(name, accessor, bin_width)
end

"""
    round_to_nearest_half_order_of_magnitude

Rounds a number to the nearest half order of magnitude, {...0.1,0.5,1,5,10,50,100,500...}
"""
function round_to_nearest_half_order_of_magnitude(w::Real)
  a = round(Int, log(w)/log(10))
  x_mid = 10.0^a
  x_lo = x_mid*0.5
  x_hi = 0.5*10.0^(a+1)
  i = indmin([abs(x_mid-w), abs(x_lo-w), abs(x_hi-w)])
  i == 1 ? x_mid : i ==2 ? x_lo : x_hi
end

"""
	infer_dimension

Constructs a Dimension suitable for the type in arr.
"""
function infer_dimension{I<:Integer}(arr::AbstractDataArray{I}, name::Symbol)
	accessor = @sprintf("function(d){return d.%s; }", name)
	Dimension(name, accessor)
end
function infer_dimension{F<:AbstractFloat}(arr::AbstractDataArray{F}, name::Symbol, desired_bincount::Int=10)

    lo,hi = extrema(arr)
    bin_width = round_to_nearest_half_order_of_magnitude((hi-lo)/desired_bincount)

	accessor = @sprintf("function(d){return Math.round(d.%s / %f)*%f; }", name, bin_width, bin_width)
	Dimension(name, accessor, bin_width)
end
function infer_dimension{S<:AbstractString}(arr::AbstractDataArray{S}, name::Symbol)
	accessor = @sprintf("function(d){return d.%s; }", name)
	Dimension(name, accessor)
end

Base.write(io::IO, dim::Dimension) = print(io, "var ", dim.name, " = cf.dimension(", dim.accessor, ");")


text:
@doc raw"""
```
Z,p,q = retrieve_surface(algorithm::Photometric, img1::AbstractArray, img2::AbstractArray, img3::AbstractArray, illumination_direction1::Vector{T} where T <: Real, illumination_direction2::Vector{T} where T <: Real, illumination_direction3::Vector{T} where T <: Real; integration_scheme::AbstractIntegrationScheme = Horn())
```
Attempts to produce the gradients and integrate them to retiereve a height map
using three-point photometric stereo.

This algorithm solves a system of linear equation formed by varying the lighting
conditions on an object while holding the viewing angle constant. Note that
the ilumination direction id defined negative of most other algorithms and as
such inouts may need to be corrected (see example).
# Output
Returns an M by N array (matching dimensions of original image) of Float `Z`
that represents the reconstructed height at the point and the gradients in
M by N arrays of Float `p` and `q`.
# Details
Let ñ₁, ñ₂ and ñ₃ define the ilumination direction let N be a 3 X 3 matrix
whos rows are formed by the row vectors of the illumination direction as per:
```math
N=\begin{bmatrix}&n_{11}&n_{12}&n_{13}\\&n_{21}&n_{22}&n_{23}\\&n_{31}&n_{32}&n_{33}\end{bmatrix}
```
and I₁, I₂ and I₃ be the intensity values at a point (x,y) and let Ĩ=[I₁,I₂,I₃]′
be the column vector formed by these values.

The surface normal at the point (x,y) can now be directly calculated using:
```math
ñ=N^{-1}Ĩ
```
where n is the surface normal.
# Arguments
The function arguments are described in more detail below.
##  `img1`
An `AbstractArray` storing the grayscale value of each pixel within
the range [0,1].
##  `img2
An `AbstractArray` storing the grayscale value of each pixel within
the range [0,1].
##  `img3`
An `AbstractArray` storing the grayscale value of each pixel within
the range [0,1].
## `illumination_direction1`
A `Vector{T} where T <: Real` that specifies the first illumination vector to be used by the
algorithm. The `illumination_direction` should be a vector of the form [x,y,z]
where x,y,z are int he range [0,1].
## `illumination_direction2`
A `Vector{T} where T <: Real` that specifies the second illumination vector to be used by the
algorithm. The `illumination_direction` should be a vector of the form [x,y,z]
where x,y,z are int he range [0,1].
## `illumination_direction3`
A `Vector{T} where T <: Real` that specifies the third illumination vector to be used by the
algorithm. The `illumination_direction` should be a vector of the form [x,y,z]
where x,y,z are int he range [0,1].
## `integration_scheme`
A keywork argument of type `AbstractIntegrationScheme` which specifies which integration
scheme is to be used to reconstruct the surface from the calculated normals.
```julia
using Images, Makie, ShapeFromShading

#generate synthetic images
n₁ = [0,0,1.0]
n₂ = [1.0,0.2,1.0]
n₃ = [0.2,0.9,1.0]
img1, img2, img3 = generate_photometric(n₁, n₂, n₃, 1, shape = SynthSphere(), r = 5)

#calculate the heightmap
Z,p,q = retrieve_surface(Photometric(), img1, img2, img3, -n₁, -n₂, -n₃, integration_scheme = Durou())

#normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)

#display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:2
surface(r, r, Z)
```
# Reference
1. R. Woodham, "Photometric Method For Determining Surface Orientation From Multiple Images", Optical Engineering, vol. 19, no. 1, 1980. [doi:10.1117/12.7972479](https://doi.org/10.1117/12.7972479)
"""
function retrieve_surface(algorithm::Photometric, img1::AbstractArray, img2::AbstractArray, img3::AbstractArray, illumination_direction1::Vector{T} where T <: Real, illumination_direction2::Vector{T} where T <: Real, illumination_direction3::Vector{T} where T <: Real; integration_scheme::AbstractIntegrationScheme = Horn())
    #setup illumination matrix
    N = zeros(Float64, 3, 3)
    for i = 1:3
        N[1,i] = illumination_direction1[i]
        N[2,i] = illumination_direction2[i]
        N[3,i] = illumination_direction3[i]
    end
    Ninv = inv(N)

    #setup image and gradient arrays
    I₁ = Float64.(img1)
    I₂ = Float64.(img2)
    I₃ = Float64.(img3)
    p = zeros(Float64, size(I₁))
    q = zeros(Float64, size(I₁))

    #calculate gradients
    for i in CartesianIndices(I₁)
        I = [I₁[i], I₂[i], I₃[i]]
        n = Ninv * I
        if n[3] == 0.0
            n = [0.0,0.0,-1.0]
        end
        p[i] = -n[1] / n[3]
        q[i] = -n[2] / n[3]
    end

    #reconstruct surface
    Z = convert_gradient(integration_scheme, p, q)
    return Z, p, q
end


text:
#EXAMPLE: broadcasting

using MPI

MPI.Init()
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)

if rank == 0
    data = [7.0,8.0,9.0,10.0]
else
    data = Float64[]
end

data = MPI.bcast(data, 0, comm)
@show rank, data

MPI.Finalize()


text:
@static if isdefined(Base, Symbol("@constprop"))
  using Base: @constprop
else
  macro constprop(_, ex); esc(ex); end
end

@generated function static_sizeof(::Type{T}) where {T}
  st = Base.allocatedinline(T) ? sizeof(T) : sizeof(Int)
  Expr(:block, Expr(:meta, :inline), StaticInt(st))
end

smax(a::StaticInt, b::StaticInt) = ifelse(gt(a, b), a, b)
smin(a::StaticInt, b::StaticInt) = ifelse(lt(a, b), a, b)

_pick_vector_width_float16(::StaticInt{RS}, ::True) where {RS} = StaticInt{RS}() ÷ StaticInt{2}()
_pick_vector_width_float16(::StaticInt{RS}, ::False) where {RS} = StaticInt{RS}() ÷ StaticInt{4}()
pick_vector_width(::Type{Float16}) = _pick_vector_width_float16(register_size(Float32), fast_half())
pick_vector_width(::Type{T}) where {T} = register_size(T) ÷ static_sizeof(T)
@inline @constprop :aggressive function _pick_vector_width(min_W, max_W, ::Type{T}, ::Type{S}, args::Vararg{Any,K}) where {K,S,T}
  _max_W = smin(max_W, pick_vector_width(T))
  _pick_vector_width(min_W, _max_W, S, args...)
end
@inline @constprop :aggressive function _pick_vector_width(min_W, max_W, ::Type{T}) where {T}
  _max_W = smin(max_W, pick_vector_width(T))
  smax(min_W, _max_W)
end
@inline @constprop :aggressive function pick_vector_width(::Type{T}, ::Type{S}, args::Vararg{Any,K}) where {T,S,K}
  _pick_vector_width(One(), register_size(), T, S, args...)
end
@inline @constprop :aggressive function pick_vector_width(::Union{Val{P},StaticInt{P}}, ::Type{T}, ::Type{S}, args::Vararg{Any,K}) where {P,T,S,K}
  _pick_vector_width(One(), smin(register_size(), nextpow2(StaticInt{P}())), T, S, args...)
end
@inline @constprop :aggressive function pick_vector_width(::Union{Val{P},StaticInt{P}}, ::Type{T}) where {P,T}
  _pick_vector_width(One(), smin(register_size(), nextpow2(StaticInt{P}())), T)
end
@inline function pick_vector_width_shift(args::Vararg{Any,K}) where {K}
  W = pick_vector_width(args...)
  W, intlog2(W)
end


text:
module FastSolversForWeightedTV

# Modules
using LinearAlgebra, AbstractLinearOperators, CUDA, NNlib, Flux, Roots

# Constant utils
const RealOrComplex{T<:Real} = Union{T,Complex{T}}

# Types and utils
include("./abstract_functional_types.jl")
include("./type_utils.jl")

# Optimization solvers
include("./optimization_solvers.jl")

# Differentiable function examples
include("./differentiable_examples.jl")

# Gradient operator
include("./gradient_operator.jl")

# Vector-field projection operator
include("./weighting_operator.jl")

# Proximable function examples
include("./proximable_examples.jl")

end

max_stars_repo_name: grizzuti/FastSolversForWeightedTV

max_issues_repo_name: grizzuti/FastSolversForWeightedTV

max_forks_repo_name: grizzuti/FastSolversForWeightedTV


text:
"""
`kfr(R, A)`

Kim-Forger activator-repressor inhibition for K = 0.
"""
kfr(R, A) = R < A ? 1 - R / A : 0.0

"""
`kfr(R, A, K)`

Kim-Forger activator-repressor inhibition.
"""
kfr(R, A, K) = (A - R - K + sqrt((A - R - K)^2 + 4*A*K)) / (2*A)

"""
`create_events(X)`

Create event matrix from specifications `X` defined as `X = [(:LL, 2*24),
(:LD, 5, 12, 12), (:DD, 5*24), (:LL, 5*24)]`
"""
function create_events(X)

    # Time counter
    current_time = 0.0

    # Initialize matrix of empty events
    events = Matrix{Float64}(undef, 0, 2)

    for x in X

        # Get code of the event
        code = Symbol(x[1])

        if code == :LD  # light-dark period

            # Iterate days
            for _ = 1:x[2]

                if !isempty(events) && events[end, 2] == current_time
                    # Extend an already existing ligt period
                    events[end, 2] = current_time + x[3]
                else
                    # Add a new light-dark cycle
                    events = [events; [current_time current_time + x[3]]]
                end

                # Increment the time counter by the duration of the added day
                current_time = current_time + x[3] + x[4]

            end

        elseif code == :DD  # constant darkness

            # Increment the time counter by the duration of the dark period
            duration = x[2]
            current_time += duration

        elseif code == :LL  # constant light

            if !isempty(events) && events[end, 2] == current_time
                # Extend an already existing ligt period
                events[end, 2] = current_time + x[2]
            else
                # Add a new light period
                events = [events; [current_time current_time + x[2]]]
            end

            # Increment the time counter the duration of the light period
            current_time += x[2]

        else  # unknown code
            error("Unknown code!")
        end
    end

    return events

end

"""
`plot_events(events; ax=gca(), ylims=[], color="#ffbfbfff", zorder=0)`

Plot events.
"""
function plot_events(events; ax=gca(), ylims=[], color="#ffbfbfff", zorder=0)

    if isempty(ylims)
        # Get ylims based on the plotted lines of the Axes

        lines = ax.get_lines()
        if isempty(lines)
            # If there are no lines, default ylims = [0, 1]
            B = 0.0
            T = 1.0
        else
            # If there are lines, ylims = [min(lines), max(lines)]
            B = Inf
            T = -Inf
            for line in lines
                y = line.get_ydata()
                B = min(B, minimum(y))
                T = max(T, maximum(y))
            end
        end
    else
        B = ylims[1]
        T = ylims[2]
    end

    n = size(events, 1)
    for i = 1:n
        # Print a box for each event
        L = events[i, 1]
        R = events[i, 2]
        ax.fill_between([L, R], [B, B], [T, T], linewidth=0.0, color=color,
            zorder=zorder)
    end

    return ax

end

"""
`estimate_period(t, x; kwargs...)`

Estimate period of a signal based on its autocorrelation function.

**Argument**
- `t`: Time vector.
- `x`: Data vector.

**Keyword Arguments**
- `n_lags`: Number of lags used for autocorrelation. Default is 1000.
- `show_plots`: Visualize the result. Default is false.

**Returns**
- `period`: Estimated period (the position of the highest).
- `peak`: Height of the highest peak.
"""
function estimate_period(t, x; n_lags=1000, show_plots=false)

    # If x is constat, return NaNs
    if all(x[1] .== x)
        return (NaN, NaN)
    end

    # Normalize input
    x = zscore(x)

    # Calculate autocorrelation
    n = min(n_lags, length(x))  # number of lags to calculate
    R = autocor(x, 0:(n-1))  # autocorrelation function

    # Find peaks of the autocorrelation function
    pr = findpeaks(R, (0:(n-1)) .* mean(diff(t)), sortstr="descend")

    # If there are no peaks, return NaNs
    if length(pr) == 0
        return (NaN, NaN)
    end

    # Get the highest peaks (peak height and location = period)
    peak = peakheights(pr)[1]
    period = peaklocations(pr)[1]

    # Visualize the results
    if show_plots
        fig, axs = subplots(2)
        axs[1].plot(t, x, color="black")
        axs[1].set_title("Time Series")
        axs[2].plot((0:(n-1)) .* mean(diff(t)), R, color="black")
        axs[2].plot([period, period], [0, peak], color="red")
        axs[2].set_title("peak = $(round(peak, digits=2)), " *
            "period = $(round(period, digits=2))")
        fig.tight_layout()
    end

    return (period, peak)

end

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
module TestHalve

import SplittablesBase
using Base.Iterators: partition, product
using OffsetArrays: OffsetArray

raw_examples = """
1:10
1:11
[1:10;]
[1:11;]
reshape(1:6, 2, 3)
adjoint(reshape(1:6, 2, 3))
transpose(reshape(1:6, 2, 3))
permutedims(reshape(1:6, 2, 3))
permutedims(reshape(1:24, 2, 3, 4), (1, 2, 3))
permutedims(reshape(1:24, 2, 3, 4), (1, 3, 2))
permutedims(reshape(1:24, 2, 3, 4), (2, 1, 3))
permutedims(reshape(1:24, 2, 3, 4), (2, 3, 1))
permutedims(reshape(1:24, 2, 3, 4), (3, 1, 2))
permutedims(reshape(1:24, 2, 3, 4), (3, 2, 1))
OffsetArray(1:15, -1)
OffsetArray(reshape(1:15, 3, 5), -1:1, 0:4)
pairs(1:10)
pairs(1:11)
pairs(reshape(1:6, 2, 3))
pairs(adjoint(reshape(1:6, 2, 3)))
"abcde"
"abcdef"
"αβγδϵ"
"αβγδϵζ"
"あいうえお"
"あいうえおか"
"aいυe𝒐"
"aいυe𝒐か"
eachindex("abcde")
eachindex("abcdef")
eachindex("αβγδϵ")
eachindex("αβγδϵζ")
eachindex("あいうえお")
eachindex("あいうえおか")
eachindex("aいυe𝒐")
eachindex("aいυe𝒐か")
(1, 2, 3, 4, 5)
(1, 2, 3, 4, 5, 6)
(a=1, b=2, c=3, d=4, e=5)
(a=1, b=2, c=3, d=4, e=5, f=6)
zip(1:3, (11, 22, 33), (a=1, b=2, c=3))
(x^2 for x in 1:5)
(x^2 for x in 1:6)
(x^2 for x in 1:5 if isodd(x))
(x^2 for x in 1:6 if isodd(x))
Iterators.filter(isodd, 1:5)
Iterators.filter(isodd, 1:6)
Iterators.flatten([1:3, 4:5])
product(1:3, 1:4)
product(1:3, 1:4, 1:5)
product(pairs(1:10), pairs(1:11))
product(1:3, pairs(1:11), (1, 2), (a=1, b=2))
partition(1:10, 1)
partition(1:10, 2)
partition(1:10, 3)
partition(1:10, 4)
enumerate([11, 22, 33, 44])
enumerate([11, 22, 33, 44, 55])
enumerate(Iterators.product(1:3, 'a':'b'))
zip(1:3, partition(1:10, 4))
zip("αβγ")
zip(enumerate([11, 22, 33, 44]), 'a':'d')
zip(enumerate([11, 22, 33, 44, 55]), 'a':'e')
zip(Iterators.product(1:3, 'a':'b'), ones(3, 2))
Iterators.reverse(1:10)
Iterators.reverse(zip(1:3, 'a':'c'))
skipmissing([1, 2, missing, 4])
skipmissing([1, 2, missing, 4, 5])
"""

# An array of `(label = ..., data = ...)`
examples = map(split(raw_examples, "\n", keepempty = false)) do code
    (label = code, data = Base.include_string(@__MODULE__, code))
end

SplittablesBase.Testing.test_ordered(examples)

end  # module


text:
module NonIterableCellMaps

using Gridap
using Gridap.Helpers

export UnimplementedMap
export NonIterableCellMap
import Base: iterate
import Gridap: evaluate
import Gridap: return_size

struct UnimplementedMap{S,M,T,N} <: Map{S,M,T,N} end

function evaluate!(
  this::UnimplementedMap{S,M,T,N},
  points::AbstractArray{<:S,M},
  v::AbstractArray{T,N}) where {S,M,T,N}
  @notimplemented
end

function return_size(::UnimplementedMap{S,M,T,N},::NTuple{M,Int}) where {S,M,T,N}
  @notimplemented
end

abstract type NonIterableCellMap{S,M,T,N} <: IterCellValue{UnimplementedMap{S,M,T,N}} end

function iterate(::NonIterableCellMap)
  _error()
end

function iterate(::NonIterableCellMap,state)
  _error()
end

function _error()
  error("Iteration is intentionally disabled for this type.")
end

end # module


text:
# You need to do to the data folder and generate Stackoverflow_data.mat
# before this will work
s = time()
M = matread("../data/Stackoverflow_data.mat")
LabelMatrix = M["LabelMat"]
LabelNames = M["LabelNames"]
H = M["H"]
Tconds = M["ClusConds"]
Tsizes = M["ClusSizes"]
order = round.(Int64,vec(sum(H,dims=2)))
d = vec(sum(H,dims=1))
volA = sum(d)
m,n = size(H)
Ht = sparse(H')
toc = time()-s
println("Done loading things into memory in $toc seconds.")

max_stars_repo_name: nveldt/SparseCardDSFM

max_issues_repo_name: nveldt/SparseCardDSFM

max_forks_repo_name: nveldt/SparseCardDSFM


text:
using SnoopCompileBot

snoop_bot(BotConfig("SnoopCompileWithMocking"));


text:
using LinearAlgebra
using Plots
gr()

"""
Global Newton Fractal for Solving Systems of Nonlinear Equations
@author Oscar Veliz
"""

"""
    F(X)

System of nonlinear equations on vector X ([x;y]).
Overwrite for your own testing, currently set to:
x^2 - y - 1
x - y^2 + 1

For a system version of z^3-1 use F(X) = [X[1]^3-3*X[1]*X[2]^2-1;3*X[1]^2*X[2]-X[2]^3]

# Example
```julia-repl
julia> F([1.0; 2.0])
2-element Array{Float64,1}:
 -2.0
 -2.0
```
"""
F(X) = [X[1]^2-X[2]-1;X[1]-X[2]^2+1]

"""
    J(X)

Jacobian matrix for F given vector X ([x;y]).
Overwrite for your own testing, currently set to:
[2x -1;1 -2y]

For a system version of z^3-1 use J(X) = [3*X[1]^2-3*X[2]^2 -6*X[1]*X[2];6*X[1]*X[2] 3*X[1]^2-3*X[2]^2]

# Example
```julia-repl
julia> J([1.0; 2.0])
2×2 Array{Float64,2}:
 2.0  -1.0
 1.0  -4.0
```
"""
J(X) = [2*X[1] -1;1 -2*X[2]]

"""
    Newton(X, eps=10^-6)

Newton's Method for the system of nonlinear equations in [`F`](@ref). Returns the number of iterations.
# Example
```julia-repl
julia> Newton([1.0;2.0])
5
```
"""
function Newton(X, eps::Float64=10^-6)
    i = 1
    FX = F(X)
    while norm(FX) > eps && i < 50
        X = X - J(X) \ FX
        FX = F(X)
        i = i + 1
    end
    i
end

"""
    GlobalNewton(X, eps=10^-6)

Global Newton's Method for the system of nonlinear equations in [`F`](@ref). Returns the number of iterations.
# Example
```julia-repl
julia> GlobalNewton([1.0;2.0])

```
"""
function GlobalNewton(X, eps::Float64=10^-6)
    FX = F(X)
    nF = norm(FX)
    i = 1
    JX = J(X)
    delta = 1000
    while norm(delta) > eps && nF > eps
        a = 1.0
        delta = JX \ FX
        T = X - a * delta
        FT = F(T)
        nFT = norm(FT)
        while nFT > nF
            a = a * 0.5
            T = X - a * delta
            FT = F(T)
            nFT = norm(FT)
        end
        delta = X - T
        X = T
        FX = FT
        nF = nFT
        JX = J(X)
        i = i + 1
    end
    i
end

"""
    Fractal(useGlobal=true, eps=10^-3)

Function for creating fractals. Call with false for normal newton fractal. Uses [`F`](@ref) and [`J`](@ref).
Saves the file to `globalnewton.svg` or `newton.svg` depending on `useGlobal`.
# Example
```julia-repl
julia> Fractal()
```
"""
function Fractal(useGlobal::Bool=true, step::Float64=.001)
    domain = -3.6:step:3.6
    range = -2.0:step:2.0
    data = rand(length(range), length(domain))
    i = 1
    for x in domain
        j = 1
        for y in range
            X = [x;y]
            try
                data[j,i] = if useGlobal GlobalNewton(X) else Newton(X) end
            catch
                data[j,i] = 50
            end
            j = j + 1
        end
        if i > length(domain)
            i = 1
        else
            i = i + 1
        end
    end
    heatmap(domain,
    range, data,
    c = cgrad([:black, :red, :yellow, :white]),
    xlabel = "x", ylabel = "y",
    title = if useGlobal "Global Newton Fractal" else "Newton Fractal" end,
    fmt = :svg)
    savefig(if useGlobal "globalnewton.svg" else "newton.svg" end)
end

"""
Main
"""

Fractal()

max_stars_repo_name: deboradeben/numerical-veliz

max_issues_repo_name: KulakovaEA/numerical-veliz

max_forks_repo_name: KulakovaEA/numerical-veliz