text:
# Social Security Mortality
maleMort = [
    0.00699,
    0.000447,
    0.000301,
    0.000233,
    0.000177,
    0.000161,
    0.00015,
    0.000139,
    0.000123,
    0.000105,
    0.000091,
    0.000096,
    0.000135,
    0.000217,
    0.000332,
    0.000456,
    0.000579,
    0.000709,
    0.000843,
    0.000977,
    0.001118,
    0.00125,
    0.001342,
    0.001382,
    0.001382,
    0.00137,
    0.001364,
    0.001362,
    0.001373,
    0.001393,
    0.001419,
    0.001445,
    0.001478,
    0.001519,
    0.001569,
    0.001631,
    0.001709,
    0.001807,
    0.001927,
    0.00207,
    0.002234,
    0.00242,
    0.002628,
    0.00286,
    0.003117,
    0.003396,
    0.003703,
    0.004051,
    0.004444,
    0.004878,
    0.005347,
    0.005838,
    0.006337,
    0.006837,
    0.007347,
    0.007905,
    0.008508,
    0.009116,
    0.009723,
    0.010354,
    0.011046,
    0.011835,
    0.012728,
    0.013743,
    0.014885,
    0.016182,
    0.017612,
    0.019138,
    0.020752,
    0.022497,
    0.024488,
    0.026747,
    0.029212,
    0.031885,
    0.034832,
    0.038217,
    0.042059,
    0.046261,
    0.050826,
    0.055865,
    0.06162,
    0.068153,
    0.075349,
    0.08323,
    0.091933,
    0.101625,
    0.112448,
    0.124502,
    0.137837,
    0.152458,
    0.168352,
    0.185486,
    0.203817,
    0.223298,
    0.243867,
    0.264277,
    0.284168,
    0.303164,
    0.320876,
    0.336919,
    0.353765,
    0.371454,
    0.390026,
    0.409528,
    0.430004,
    0.451504,
    0.474079,
    0.497783,
    0.522673,
    0.548806,
    0.576246,
    0.605059,
    0.635312,
    0.667077,
    0.700431,
    0.735453,
    0.772225,
    0.810837,
    0.851378,
    0.893947,
]

femaleMort = [
    0.005728,
    0.000373,
    0.000241,
    0.000186,
    0.00015,
    0.000133,
    0.000121,
    0.000112,
    0.000104,
    0.000098,
    0.000094,
    0.000098,
    0.000114,
    0.000143,
    0.000183,
    0.000229,
    0.000274,
    0.000314,
    0.000347,
    0.000374,
    0.000402,
    0.000431,
    0.000458,
    0.000482,
    0.000504,
    0.000527,
    0.000551,
    0.000575,
    0.000602,
    0.00063,
    0.000662,
    0.000699,
    0.000739,
    0.00078,
    0.000827,
    0.000879,
    0.000943,
    0.00102,
    0.001114,
    0.001224,
    0.001345,
    0.001477,
    0.001624,
    0.001789,
    0.001968,
    0.002161,
    0.002364,
    0.002578,
    0.0028,
    0.003032,
    0.003289,
    0.003559,
    0.003819,
    0.004059,
    0.004296,
    0.004556,
    0.004862,
    0.005222,
    0.005646,
    0.006136,
    0.006696,
    0.007315,
    0.007976,
    0.008676,
    0.009435,
    0.010298,
    0.011281,
    0.01237,
    0.013572,
    0.014908,
    0.01644,
    0.018162,
    0.020019,
    0.022003,
    0.024173,
    0.026706,
    0.029603,
    0.032718,
    0.036034,
    0.039683,
    0.043899,
    0.048807,
    0.054374,
    0.060661,
    0.067751,
    0.075729,
    0.084673,
    0.094645,
    0.105694,
    0.117853,
    0.131146,
    0.145585,
    0.161175,
    0.17791,
    0.195774,
    0.213849,
    0.231865,
    0.249525,
    0.266514,
    0.282504,
    0.299455,
    0.317422,
    0.336467,
    0.356655,
    0.378055,
    0.400738,
    0.424782,
    0.450269,
    0.477285,
    0.505922,
    0.536278,
    0.568454,
    0.602561,
    0.638715,
    0.677038,
    0.71766,
    0.76072,
    0.806363,
    0.851378,
    0.893947,
]

bowers_illus_life_vec = [20.4217,1.3431,1.2237,1.1239,1.0421,0.9770,0.9269,0.8904,0.8660,0.8522,0.8475,0.8504,0.8594,0.8730,0.8898,0.9081,0.9282,0.9502,0.9744,1.0009,1.0299,1.0618,1.0967,1.1350,1.1770,1.2330,1.2735,1.3288,1.3895,1.4560,1.5289,1.6089,1.6965,1.7927,1.8980,2.0136,2.1402,2.2791,2.4313,2.5982,2.7812,2.9818,3.2017,3.4427,3.7070,3.9966,4.3141,4.6621,5.0436,5.4617,5.9199,6.4221,6.9724,7.5755,8.2364,8.9605,9.7538,10.6230,11.5752,12.6181,13.7604,15.0114,16.3813,17.8812,19.5231,21.3203,23.2871,25.4391,27.7932,30.3680,33.1833,36.2608,39.6240,43.2982,47.3108,51.6911,56.4708,61.6840,67.3671,73.5589,80.3009,87.6369,95.6134,104.2794,113.6860,123.8867,134.9367,146.8926,159.8121,173.7533,188.7738,204.9298,222.2749,240.8589,260.7257,281.9122,304.4456,328.3410,353.5993,380.2041,408.1188,437.2837,467.6133,498.9935,531.2793,564.2937,597.8266,631.6360,665.4495,698.9685,731.8742,] ./ 1000

bowers_illus_life = mortality_vector(
    [decrement(MortalityTables.Makeham(a=0.05/1000,b=log(10^0.04),c=0.7/1000),age,age+1) for age in 13:111],
    start_age=13
    )

# AMLCR 2nd ed, §4.3
AMLCR_std = mortality_vector(
    [decrement(MortalityTables.Makeham(a=2.7*10^(-6),b=log(1.124),c=0.00022),age,age+1) for age in 13:120],
    start_age=13
    )


text:
module CyclicDiagonal

using SpecialMatrices, LinearAlgebra

export iterDCD

# Iteratively find an approximate factorization of A as a product D1 * C * D2
# where D1 and D2 are diagonal matrices, and C is a circulant matrix
function iterDCD(A::Matrix, max_iter = 1e5, tol = 1e-5)
	s_col = ones(size(A,1))
	s_row = Array{Float64}(undef, 0, 2)
	mu_row = Array{Float64}(undef, 0, 2)
	mu_col = Array{Float64}(undef, 0, 2)

	d = Inf
	i = 0
	while i < max_iter && d > tol
		mu_row, s_row = find_mu_and_s( A .* (s_col.^-1)', max_iter, tol)
		mu_col, s_col = find_mu_and_s( ((s_row.^-1) .* A)', max_iter, tol )
		d = norm( [ mu_col[1]; mu_col[end:-1:2] ]  - mu_row )
		i = i + 1
	end

	D1 = Diagonal( s_row[:] )
	C =  Circulant( mu_col[:] )
	D2 = Diagonal( s_col[:] )

	D1, C, D2
end

# Find left side diagonal scaling matrix with diagonal s and circulant diagonal matrix factorization approximations
function find_mu_and_s(A, max_iter , tol)

    uncirc!(A)
    mu = sum(A, dims=1)'/size(A,1)
    mu_old = fill(Inf, size(mu))
	s = ones(size(A,1))

    i = 0
    while i < max_iter && norm( mu_old - mu ) > tol

        s = A*mu/(mu'*mu)
        mu_old = mu
        mu = sum(A ./ s, dims=1)'/size(A,1)

        i = i + 1
    end

	mu, s
end

# Circulate rows of A
function circ!(A)
    for i in 1:size(A,1)
        A[i,:] = circshift(A[i,:], (i-1))
    end
	A
end

# Uncirculate rows of A
function uncirc!(A)
    for i in 1:size(A,1)
        A[i,:] = circshift(A[i,:], (1-i))
    end
	A
end

end # module

ext: jl


text:
function upload_iris_to_s3(bucket_name)
    # iris_download_path = "https://raw.githubusercontent.com/banyan-team/banyan-julia/v0.1.3/BanyanDataFrames/test/res/iris.csv"
    iris_download_path = "https://raw.githubusercontent.com/banyan-team/banyan-julia/v0.1.3/BanyanDataFrames/test/res/iris.csv"
    iris_local_path = download(iris_download_path)
    df = CSV.read(iris_local_path, DataFrames.DataFrame)
    # Duplicate df six times and change the species names
    species_list = df[:, :species]
    df = reduce(vcat, [df, df, df, df, df, df])
    for i = 4:18
        species_list = append!(species_list, Base.fill("species_$(i)", 50))
    end
    df[:, :species] = species_list
    CSV.write("iris.csv", df)
    verify_file_in_s3(bucket_name, "iris_large.csv", p"iris.csv")
    verify_file_in_s3(
        bucket_name,
        "iris_species_info.csv",
        "https://raw.githubusercontent.com/banyan-team/banyan-julia/v0.1.3/BanyanDataFrames/test/res/iris_species_info.csv",
    )
end

# TODO: Add test for setting multiple columns with a matrix

@testset "Basic data analytics on a small dataset" begin
    run_with_job("Indexing small dataset") do job
        bucket = get_cluster_s3_bucket_name(get_cluster_name())
        # upload_iris_to_s3(bucket)
        iris = read_csv("s3://$(bucket)/iris_large.csv")

        # Select rows that have a thin petal (i.e., length > 4 * width). Select only the petal columns.
        res = collect(
            sort(
                iris[
                    map(
                        (pl, pw) -> pl > 4 * pw,
                        iris[:, :petal_length],
                        iris[:, :petal_width],
                    ),
                    :,
                ][
                    :,
                    [3, 4],
                ],
            ),
        )
        @test round(sum(res[:, :petal_width])) == 61.0
        @test res[1, 1] == 1.0
        @test res[36, 2] == 0.2
        @test res[44, 1] == 1.3
        res_row_1 = res[256, [1, 2]]
        @test collect(res_row_1) == [1.9, 0.4]
        res_row_2 = res[nrow(res), [1, 2]]
        @test collect(res_row_2) == [4.1, 1.0]
    end

    run_with_job("Filtering small dataset") do job
        bucket = get_cluster_s3_bucket_name(get_cluster_name())
        # upload_iris_to_s3(bucket)
        iris = read_csv("s3://$(bucket)/iris_large.csv")

        # Filter
        iris_filtered =
            sort(filter(row -> (row.sepal_length > 5.0 && row.sepal_width < 3.0), iris))
        iris_filtered_collect = collect(iris_filtered)
        @test nrow(iris_filtered) == 306
        @test iris_filtered_collect[1, :sepal_width] == 2.5
        @test iris_filtered_collect[46, :petal_length] == 3.9
        @test iris_filtered_collect[51, :sepal_length] == 5.6
        @test collect(iris_filtered_collect[298, :]) == [7.7, 2.6, 6.9, 2.3, "species_6"]

        # Unique
        @test nrow(unique(iris[:, [:species]])) == 18
        unique_iris = sort(unique(iris))
        @test nrow(unique_iris) == 882
        @test nrow(
            unique(iris[:, [:sepal_length, :sepal_width, :petal_length, :petal_width]]),
        ) == 147

        # Nonunique before aggregation
        nonunique_idxs = nonunique(iris)
        nonunique_iris = collect(sort(iris[nonunique_idxs, :]))
        @show collect(iris[nonunique_idxs, :])
        @show nonunique_iris
        @show iris
        @show typeof(iris)
        @show collect(iris)
        @show collect(nonunique_idxs)
        @test collect(first(nonunique_iris)) == [4.9, 3.1, 1.5, 0.1, "setosa"]
        @test collect(nonunique_iris[10, :]) == [4.9, 3.1, 1.5, 0.1, "species_4"]
        @test collect(last(nonunique_iris)) == [5.8, 2.7, 5.1, 1.9, "virginica"]
        nonunique_species = collect(nonunique(iris, :species))
        @test sum(nonunique_species) == 882

        # Nonunique after aggregation
        @test collect(sum(nonunique_idxs)) == 18
        # @show collect(sum(nonunique_idxs))
        # @show sample(nonunique_idxs)
        # @show collect(nonunique_idxs)
        nonunique_iris = collect(sort(iris[nonunique_idxs, :]))
        @show collect(iris[nonunique_idxs, :])
        @show nonunique_iris
        @show collect(iris)
        @show collect(nonunique_idxs)
        @test collect(first(nonunique_iris)) == [4.9, 3.1, 1.5, 0.1, "setosa"]
        @test collect(nonunique_iris[10, :]) == [4.9, 3.1, 1.5, 0.1, "species_4"]
        @test collect(last(nonunique_iris)) == [5.8, 2.7, 5.1, 1.9, "virginica"]
        nonunique_species = collect(nonunique(iris, :species))
        @test sum(nonunique_species) == 882

        # Nonunique after _another_ aggregation
        @test collect(sum(nonunique_idxs)) == 18
        # @show collect(sum(nonunique_idxs))
        # @show sample(nonunique_idxs)
        # @show collect(nonunique_idxs)
        nonunique_iris = collect(sort(iris[nonunique_idxs, :]))
        @show collect(iris[nonunique_idxs, :])
        @show nonunique_iris
        @show collect(iris)
        @show collect(nonunique_idxs)
        @test collect(first(nonunique_iris)) == [4.9, 3.1, 1.5, 0.1, "setosa"]
        @test collect(nonunique_iris[10, :]) == [4.9, 3.1, 1.5, 0.1, "species_4"]
        @test collect(last(nonunique_iris)) == [5.8, 2.7, 5.1, 1.9, "virginica"]
        nonunique_species = collect(nonunique(iris, :species))
        @test sum(nonunique_species) == 882

        # Dropmissing, allowmissing, disallowmissing
        @test nrow(dropmissing(iris)) == 900
        iris_allowmissing = allowmissing(iris)
        #setindex!(iris_allowmissing, missing, 1, 1)
        #setindex!(iris_allowmissing, missing, 123, 2)
        #setindex!(iris_allowmissing, missing, 150, 5)
        #iris_cleaned = dropmissing(iris_allowmissing)
        #iris_cleaned = disallowmissing(iris_cleaned)
        #@test nrow(iris_cleaned) == 147
        #@test collect(iris_cleaned[123, :])[:sepal_length] == 6.7
        #iris_cleaned = collect(iris_cleaned)
        #@test first(iris_cleaned)[:sepal_width] == 3.0
        #@test last(iris_cleaned)[:petal_width] == 2.3
    end

    run_with_job("Sorting small dataset") do job
        bucket = get_cluster_s3_bucket_name(get_cluster_name())
        # upload_iris_to_s3(bucket)
        iris = read_csv("s3://$(bucket)/iris_large.csv")

        # Sort forward and backward
        iris_sorted = sort(iris, :sepal_width)
        iris_sorted_reverse = sort(iris, :sepal_width, rev = true)

        iris_sorted = collect(iris_sorted)
        iris_sorted_reverse = collect(iris_sorted_reverse)

        @test iris_sorted[[1, 142, nrow(iris_sorted)], [:sepal_width, :petal_length]] ==
              DataFrames.DataFrame(:sepal_width => [2.0, 2.6, 4.4], :petal_length => [3.5, 4.0, 1.5])
        @test iris_sorted_reverse[
            [1, 142, nrow(iris_sorted)],
            [:sepal_width, :petal_length],
        ] == DataFrames.DataFrame(:sepal_width => [4.4, 3.5, 2.0], :petal_length => [1.5, 1.3, 3.5])
    end

    run_with_job("Join and group-by-aggregate small dataset") do job
        # Join and groupby aggregration
        bucket = get_cluster_s3_bucket_name(get_cluster_name())
        # upload_iris_to_s3(bucket)
        iris = read_csv("s3://$(bucket)/iris_large.csv")
        species_info = read_csv("s3://$(bucket)/iris_species_info.csv")
        iris_new = innerjoin(iris, species_info, on = :species)
        @test ncol(iris_new) == 6

        gdf_species = groupby(iris_new, :species)
        gdf_region = groupby(iris_new, :region)

        gr_cols = collect(groupcols(gdf_region))
        gs_cols = collect(groupcols(gdf_species))
        vr_cols = Set(collect(valuecols(gdf_region)))

        @test gr_cols == [:region]
        @test gs_cols == [:species]
        @test vr_cols ==
              Set([:sepal_length, :sepal_width, :petal_length, :petal_width, :species])

        lengths_species =
            collect(sort(combine(gdf_species, :petal_length => mean), :petal_length_mean))
        counts_species = collect(sort(combine(gdf_region, nrow), :region))
        lengths_region =
            collect(sort(combine(gdf_region, :petal_length => mean), :petal_length_mean))
        counts_region = collect(sort(combine(gdf_region, nrow), :region))

        @test lengths_species[:, :petal_length_mean] ==
              lengths_region[:, :petal_length_mean] ==
              [
                  1.464,
                  1.464,
                  1.464,
                  1.464,
                  1.464,
                  1.464,
                  4.26,
                  4.26,
                  4.26,
                  4.26,
                  4.26,
                  4.26,
                  5.552,
                  5.552,
                  5.552,
                  5.552,
                  5.552,
                  5.552,
              ]
        @test counts_species[:, :nrow] == counts_region[:, :nrow] == Base.fill(50, 18)
    end

    run_with_job("Group and map over small dataset") do job
        bucket = get_cluster_s3_bucket_name(get_cluster_name())
        # upload_iris_to_s3(bucket)
        iris = read_csv("s3://$(bucket)/iris_large.csv")

        # Prepare some larger data. Groupby by rounded petal length into 54 groups.
        iris = innerjoin(iris, iris[:, [:species]], on = :species)
        setindex!(iris, map(pl -> round(pl), iris[:, :petal_length]), :, :petal_length_rounded)
        gdf = groupby(iris, [:species, :petal_length_rounded])
        @test length(gdf) == 54

        # Select
        iris_select = sort(
            select(gdf, :, [:petal_length] => (pl) -> pl .- mean(pl)),
            :petal_length_function,
        )
        petal_length_function_sum = round(collect(reduce(-, iris_select[:, :petal_length_function])), digits=2)
        @test petal_length_function_sum == -52.71

        #@test round.(
        #    collect(
        #        iris_select[:, :petal_length_function][[1, 7488, 34992, nrow(iris_select)]],
        #    ),
        #    digits = 3,
        #) == [-0.639, -0.186, 0.115, 0.592]

        # Transform
        iris_new = transform(gdf, :species => x -> "iris-" .* x)
        species_names = collect(sort(unique(iris_new[:, "species_function"])))
        @test species_names == [
            "iris-setosa",
            "iris-species_10",
            "iris-species_11",
            "iris-species_12",
            "iris-species_13",
            "iris-species_14",
            "iris-species_15",
            "iris-species_16",
            "iris-species_17",
            "iris-species_18",
            "iris-species_4",
            "iris-species_5",
            "iris-species_6",
            "iris-species_7",
            "iris-species_8",
            "iris-species_9",
            "iris-versicolor",
            "iris-virginica",
        ]

        # Split-Apply-Combine
        iris_mins =
            collect(sort(combine(gdf, :petal_length => minimum), :petal_length_minimum))
        @test iris_mins[:, :petal_length_minimum] == [
            1.0,
            1.0,
            1.0,
            1.0,
            1.0,
            1.0,
            1.5,
            1.5,
            1.5,
            1.5,
            1.5,
            1.5,
            3.0,
            3.0,
            3.0,
            3.0,
            3.0,
            3.0,
            3.5,
            3.5,
            3.5,
            3.5,
            3.5,
            3.5,
            4.5,
            4.5,
            4.5,
            4.5,
            4.5,
            4.5,
            4.6,
            4.6,
            4.6,
            4.6,
            4.6,
            4.6,
            4.8,
            4.8,
            4.8,
            4.8,
            4.8,
            4.8,
            5.5,
            5.5,
            5.5,
            5.5,
            5.5,
            5.5,
            6.6,
            6.6,
            6.6,
            6.6,
            6.6,
            6.6,
        ]

        # Subset
        long_petal_iris = sort(
            combine(
                groupby(subset(gdf, :petal_length => pl -> pl .>= mean(pl)), :species),
                nrow,
            ),
        )
        long_petal_iris_nrow = collect(long_petal_iris)[:, :nrow]
        long_petal_iris_species = collect(long_petal_iris)[:, :species]
        @test long_petal_iris_nrow == [
            1250,
            1250,
            1250,
            1200,
            1250,
            1250,
            1200,
            1250,
            1250,
            1200,
            1250,
            1250,
            1200,
            1250,
            1250,
            1200,
            1250,
            1200,
        ]
        @test long_petal_iris_species == [
            "setosa",
            "species_10",
            "species_11",
            "species_12",
            "species_13",
            "species_14",
            "species_15",
            "species_16",
            "species_17",
            "species_18",
            "species_4",
            "species_5",
            "species_6",
            "species_7",
            "species_8",
            "species_9",
            "versicolor",
            "virginica",
        ]
    end

    run_with_job("Simple group-by aggregation on small dataset") do job
        bucket = get_cluster_s3_bucket_name(get_cluster_name())
        # upload_iris_to_s3(bucket)
        iris = read_csv("s3://$(bucket)/iris_large.csv")
        gdf = groupby(iris, :species)
        lengths = collect(sort(combine(gdf, :petal_length => mean)))
        counts = collect(combine(gdf, nrow))

        @show lengths
        @show counts

        @test length(gdf) == 18
        @test lengths[:, :petal_length_mean] == [
            1.464,
            1.464,
            4.26,
            5.552,
            1.464,
            4.26,
            5.552,
            1.464,
            4.26,
            5.552,
            1.464,
            4.26,
            5.552,
            1.464,
            4.26,
            5.552,
            4.26,
            5.552,
        ]
        @test counts[:, :nrow] == Base.fill(50, 18)
    end

    run_with_job("Inner Join") do job
        bucket = get_cluster_s3_bucket_name(get_cluster_name())
        # upload_iris_to_s3(bucket)
        iris = read_csv("s3://$(bucket)/iris_large.csv")
        species_info = read_csv("s3://$(bucket)/iris_species_info.csv")

        # Rename columns of species_info
        iris_species_info =
            rename(species_info, :species => :species_info, :region => :region_info)

        # Join iris and species_info on the :species/:species_info column
        iris_joined = innerjoin(
            iris,
            iris_species_info,
            on = :species => :species_info,
            renamecols = uppercase => uppercase,
        )
        @test Set(names(iris_joined)) == Set([
            "SEPAL_WIDTH",
            "PETAL_WIDTH",
            "REGION_INFO",
            "PETAL_LENGTH",
            "SEPAL_LENGTH",
            "species",
        ])
        @test nrow(iris_joined) == 900

        # Join iris on multiple columns
        iris[:, :petal_length_rounded] = map(pl -> round(pl), iris[:, :petal_length])
        species_info[:, :pl_rounded] =
            [1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4]
        iris_joined = sort(
            innerjoin(
                iris,
                species_info,
                on = [:species, :petal_length_rounded => :pl_rounded],
                :sepal_length,
            ),
        )
        @test Set(names(iris_joined)) == Set([
            "sepal_length",
            "sepal_width",
            "petal_length",
            "petal_width",
            "species",
            "petal_length_rounded",
            "region",
        ])
        @test nrow(iris_joined) == 146
        res = collect(iris_joined)
        @test collect(res[1, :]) == [4.3, 3.0, 1.1, 0.1, "setosa", 1.0, "Arctic"]
        @test collect(res[121, :]) == [6.6, 2.9, 4.6, 1.3, "species_5", 5.0, "region_5"]
        @test collect(res[146, :]) == [7.9, 3.8, 6.4, 2.0, "species_6", 6.0, "region_6"]
    end
end

@testset "Complex data analytics on a small dataset" begin
    run_with_job("Multiple evaluations together (test 1) on small dataset") do job
        bucket = get_cluster_s3_bucket_name(get_cluster_name())
        # upload_iris_to_s3(bucket)
        iris = read_csv("s3://$(bucket)/iris_large.csv")

        # Compute the min-max normalized petal_length (within each species)
        # for each flower and sort by this value. Select those with normalized
        # petal length greater than 0.5.

        gdf = groupby(iris, :species)

        # Method 1
        result_1 = transform(
            gdf,
            :petal_length => x -> (x .- minimum(x)) ./ (maximum(x) - minimum(x)),
        )
        result_1 = rename(result_1, :petal_length_function => :petal_length_normalized)
        result_1 = filter(row -> row.petal_length_normalized < row.petal_length, result_1)
        result_1 = sort(result_1, :petal_length_normalized)

        # Method 2
        min_max = combine(gdf, :petal_length => minimum, :petal_length => maximum)
        iris_new = innerjoin(iris, min_max, on = :species)
        iris_new[:, :petal_length_normalized] = map(
            (l, min, max) -> (pl - min) / (max - min),
            iris_new[:, :petal_length],
            iris_new[:, :petal_length_minimum],
            iris_new[:, :petal_length_maximum],
        )
        result_2 = filter(row -> row.petal_length_normalized < row.petal_length, result_2)
        result_2 = sort(result_2, :petal_length_normalized)

        @test nrow(result_1) == nrow(result_2) == 462

        result_1 = collect(result_1)
        result_2 = collect(result_2)

        @test isapprox(result_1[7, :petal_length_normalized], 0.52380, atol = 1e-4)
        @test isapprox(result_2[7, :petal_length_normalized], 0.52380, atol = 1e-4)
        @test isapprox(result_1[200, :petal_length_normalized], 0.66666, atol = 1e-4)
        @test isapprox(result_2[200, :petal_length_normalized], 0.66666, atol = 1e-4)
        @test first(result_1)[:species] == first(result_2)[:species] == "versicolor"
        @test last(result_1)[:species] == last(result_2)[:species] == "species_18"
        @test last(result_1)[:petal_length_normalized] ==
              last(result_1)[:petal_length_normalized] ==
              1.0
    end

    run_with_job("Multiple evaluations together (test 2) on small dataset") do job
        bucket = get_cluster_s3_bucket_name(get_cluster_name())
        # upload_iris_to_s3(bucket)
        iris = read_csv("s3://$(bucket)/iris_large.csv")

        # Inner join with itself on species to increase cardinality. Group by rounded sepal length
        # and compute number of rows. Sort by number of rows.
        iris = innerjoin(iris, iris[:, [:species]], on = :species)
        iris[:, :sepal_length_function] = map(sl -> round(sl), iris[:, :sepal_length])
        iris_sepal_length_groups = write_to_disk(
            sort(
                combine(
                    groupby(
                        rename(iris, :sepal_length_function => :sepal_length_rounded),
                        :sepal_length_rounded,
                    ),
                    nrow,
                ),
                :nrow,
            ),
        )
        res = collect(iris_sepal_length_groups)
        @test res[:, :nrow] == [1500, 1800, 7200, 14100, 20400]
        @test res[:, :sepal_length_rounded] == [4.0, 8.0, 7.0, 5.0, 6.0]
    end

    run_with_job("Multiple evaluations apart on small dataset") do job
        bucket = get_cluster_s3_bucket_name(get_cluster_name())
        # upload_iris_to_s3(bucket)
        iris = read_csv("s3://$(bucket)/iris_large.csv")

        # Inner join twice with itself on species to increase cardinality. Group by species and rounded petal length.
        # Compute cardinality of each group.
        # Also compute deviation from mean of petal_width of each group, filter for those with a deviation less that +/- 0.2,
        # and group on species.
        iris_large = groupby(
            innerjoin(
                innerjoin(iris, iris[:, [:species]], on = :species),
                iris[:, [:species]],
                on = :species,
            ),
            :species,
        )
        iris_large_gdf = groupby(
            rename(
                select(
                    iris_large,
                    :,
                    [:species, :petal_length] =>
                        (s, pl) -> s .* "-pl=" .* string.(round.(pl)),
                ),
                :species_petal_length_function => :species_pl,
            ),
            :species_pl,
        )
        @test length(iris_large_gdf) == 54

        lengths = collect(sort(combine(iris_large_gdf, nrow)))
        @test lengths[10, "species_pl"] == "species_12-pl=6.0"
        @test lengths[:, :nrow] == [
            57500,
            67500,
            57500,
            67500,
            7500,
            82500,
            35000,
            2500,
            52500,
            60000,
            10000,
            57500,
            67500,
            7500,
            82500,
            35000,
            2500,
            52500,
            60000,
            10000,
            57500,
            67500,
            7500,
            82500,
            35000,
            2500,
            52500,
            60000,
            10000,
            57500,
            67500,
            7500,
            82500,
            35000,
            2500,
            52500,
            60000,
            10000,
            57500,
            67500,
            7500,
            82500,
            35000,
            2500,
            52500,
            60000,
            10000,
            7500,
            82500,
            35000,
            2500,
            52500,
            60000,
            10000,
        ]

        iris_large_dev = sort(
            transform(iris_large_gdf, :petal_width => w -> (mean(w) .- w)),
            :petal_width_function,
        )
        lengths_dev = collect(
            filter(
                row -> (
                    row.petal_width_function >= -0.2 && row.petal_width_function <= 0.2
                ),
                iris_large_dev,
            ),
        )
        @test nrow(lengths_dev) == 1650000
    end
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
###DEPRECATED###WILL BE REMOVE SOON###

"""
    find_metabolite_by_formula(query::String, database)

Searching function for every metabolite that matches the chemical formula. Eases process of finding the desired metabolites faster in larger data arrays.
The `query` string should preferably be formatted according to the Hill notation. The `database` needs to be of one of the model types.
Please note that the pattern matching is exact and thus is prone to typing errors or alternate formulation (e.g. "NaCl" instead of "ClNa").
"""
function find_metabolite_by_formula(database, query::String)
    res_ls = String[]
    for met in metabolites(database)
        refstr = ""
        for (elem, elnum) in metabolite_formula(database, met)
            if elnum != 1
                refstr = string(refstr, elem, elnum)
            else
                refstr = string(refstr, elem)
            end
        end
        if refstr == query
            push!(res_ls, met)
        end
    end
    if length(res_ls) == 0
        println("No metabolite with matching formula found.")
    else
        return res_ls
    end
end

function print_formula(stoichiometry)
    leftstr = string()
    rightstr = string()
    for (k, v) in stoichiometry
        if v < 0
            leftstr = string(leftstr, string(abs(v)), k, " + ")
        elseif v > 0
            rightstr = string(rightstr, string(abs(v)), k, " + ")
        end
    end
    println(string(leftstr[1:end-3], " => ", rightstr[1:end-3]))
end

max_stars_repo_name: VictorVeraFrazao/OptStoic

max_issues_repo_name: VictorVeraFrazao/OptStoic

max_forks_repo_name: VictorVeraFrazao/OptStoic


text:
@testset "ElementVariableArray2MILP: $(fct_type), dimension $(dim), $(T)" for fct_type in ["vector of variables", "vector affine function"], dim in [2, 3], T in [Int, Float64]
    mock = MOIU.MockOptimizer(MILPModel{T}())
    model = COIB.ElementVariableArray2MILP{T}(mock)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{T},
        MOI.EqualTo{T},
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.ElementVariableArray,
    )

    x_index, _ = MOI.add_constrained_variable(model, MOI.Integer())
    if T == Int
        x_value, _ = MOI.add_constrained_variable(model, MOI.Integer())
        x_array, _ = MOI.add_constrained_variables(model, [MOI.Integer() for i in 1:dim])
    elseif T == Float64
        x_value = MOI.add_variable(model)
        x_array = MOI.add_variables(model, dim)
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables([x_value, x_index, x_array...])
    elseif fct_type == "vector affine function"
        MOIU.vectorize([x_value, x_index, x_array...])
    else
        @assert false
    end

    @test_throws AssertionError MOI.add_constraint(model, fct, CP.ElementVariableArray(dim))

    for i in 1:dim
        MOI.add_constraint(model, x_array[i], MOI.Interval(zero(T), one(T)))
    end

    c = MOI.add_constraint(model, fct, CP.ElementVariableArray(dim))

    @test MOI.is_valid(model, x_value)
    @test MOI.is_valid(model, x_index)
    for i in 1:dim
        @test MOI.is_valid(model, x_array[i])
    end
    @test MOI.is_valid(model, c)

    bridge = MOIBC.bridges(model)[MOI.ConstraintIndex{MOI.VectorOfVariables, CP.ElementVariableArray}(-1)]

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VectorOfVariables, CP.ElementVariableArray) == typeof(bridge)
        if T == Int
            @test MOIB.added_constrained_variable_types(typeof(bridge)) == [(MOI.ZeroOne,), (MOI.Integer,)]
        elseif T == Float64
            @test MOIB.added_constrained_variable_types(typeof(bridge)) == [(MOI.ZeroOne,)]
        end
        @test MOIB.added_constraint_types(typeof(bridge)) == [
            (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
            (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}),
            (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}),
        ]

        @test MOI.get(bridge, MOI.NumberOfVariables()) == 2 * dim
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VariableIndex, MOI.ZeroOne}()) == dim
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VariableIndex, MOI.Integer}()) == ((T == Int) ? dim : 0)
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}()) == 3
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}()) == dim
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}()) == dim

        @test MOI.get(bridge, MOI.ListOfVariableIndices()) == vcat(bridge.vars_unary, bridge.vars_product)
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VariableIndex, MOI.ZeroOne}()) == bridge.vars_unary_bin
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VariableIndex, MOI.Integer}()) == bridge.vars_product_int
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}()) == [bridge.con_unary, bridge.con_choose_one, bridge.con_value]
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}()) == bridge.con_product_lt
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}()) == bridge.con_product_gt
    end

    @testset "Set of binary variables" begin
        @test length(bridge.vars_unary) == dim
        @test length(bridge.vars_unary_bin) == dim

        for i in 1:dim
            @test MOI.is_valid(model, bridge.vars_unary[i])
            @test MOI.is_valid(model, bridge.vars_unary_bin[i])
            @test MOI.get(model, MOI.ConstraintFunction(), bridge.vars_unary_bin[i]) == bridge.vars_unary[i]
            @test MOI.get(model, MOI.ConstraintSet(), bridge.vars_unary_bin[i]) == MOI.ZeroOne()
        end
    end

    @testset "Set of product variables" begin
        @test length(bridge.vars_product) == dim
        @test length(bridge.vars_product_int) == ((T == Int) ? dim : 0)

        for i in 1:dim
            @test MOI.is_valid(model, bridge.vars_product[i])

            if T == Int
                @test MOI.is_valid(model, bridge.vars_product_int[i])
                @test MOI.get(model, MOI.ConstraintFunction(), bridge.vars_product_int[i]) == bridge.vars_product[i]
                @test MOI.get(model, MOI.ConstraintSet(), bridge.vars_product_int[i]) == MOI.Integer()
            end
        end
    end

    @testset "Unary decomposition of the index" begin
        @test MOI.is_valid(model, bridge.con_unary)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_unary)
        @test length(f.terms) == dim + 1
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con_unary) == MOI.EqualTo(zero(T))

        t = f.terms[1]
        @test t.coefficient === -one(T)
        @test t.variable == x_index

        for i in 1:dim
            t = f.terms[i + 1]
            @test t.coefficient === T(i)
            @test t.variable == bridge.vars_unary[i]
        end
    end

    @testset "Unicity of the unary decomposition" begin
        @test MOI.is_valid(model, bridge.con_choose_one)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_choose_one)
        @test length(f.terms) == dim
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con_choose_one) == MOI.EqualTo(one(T))

        for i in 1:dim
            t = f.terms[i]
            @test t.coefficient === one(T)
            @test t.variable == bridge.vars_unary[i]
        end
    end

    @testset "Value at the index" begin
        @test MOI.is_valid(model, bridge.con_value)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_value)
        @test length(f.terms) == dim + 1
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con_value) == MOI.EqualTo(zero(T))

        t = f.terms[1]
        @test t.coefficient === -one(T)
        @test t.variable == x_value

        for i in 1:dim
            t = f.terms[i + 1]
            @test t.coefficient === one(T)
            @test t.variable == bridge.vars_product[i]
        end
    end

    @testset "Constraining the product variables" begin
        @test length(bridge.con_product_lt) == dim
        @test length(bridge.con_product_gt) == dim

        for i in 1:dim
            @test MOI.is_valid(model, bridge.con_product_lt[i])
            flt = MOI.get(model, MOI.ConstraintFunction(), bridge.con_product_lt[i])
            @test length(flt.terms) == 2
            @test MOI.get(model, MOI.ConstraintSet(), bridge.con_product_lt[i]) == MOI.LessThan(zero(T))

            t1 = flt.terms[1]
            @test t1.coefficient === -one(T)
            @test t1.variable == bridge.vars_unary[i]

            t2 = flt.terms[2]
            @test t2.coefficient === one(T)
            @test t2.variable == bridge.vars_product[i]

            @test MOI.is_valid(model, bridge.con_product_gt[i])
            fgt = MOI.get(model, MOI.ConstraintFunction(), bridge.con_product_gt[i])
            @test length(fgt.terms) == 3
            @test MOI.get(model, MOI.ConstraintSet(), bridge.con_product_gt[i]) == MOI.GreaterThan(one(T))

            t1 = fgt.terms[1]
            @test t1.coefficient === -one(T)
            @test t1.variable == x_array[i]

            t2 = fgt.terms[2]
            @test t2.coefficient === one(T)
            @test t2.variable == bridge.vars_unary[i]

            t3 = fgt.terms[3]
            @test t3.coefficient === one(T)
            @test t3.variable == bridge.vars_product[i]
        end
    end
end


text:
using UnicodeFun
using Test

# write your own tests here
@test UnicodeFun.to_superscript(-1234567890) == "⁻¹²³⁴⁵⁶⁷⁸⁹⁰"
@test UnicodeFun.to_subscript(-1234567890) == "₋₁₂₃₄₅₆₇₈₉₀"

latexstring = "\\alpha^2 \\cdot \\alpha^{2+3} \\equiv \\alpha^7"
@test to_latex(latexstring) == "α² ⋅ α²⁺³ ≡ α⁷"
latexstring = "\\itA \\in \\bbR^{nxn}, \\bfv \\in \\bbR^n, \\lambda_i \\in \\bbR: \\itA\\bfv = \\lambda_i\\bfv"
@test to_latex(latexstring) == "𝐴 ∈ ℝⁿˣⁿ, 𝐯 ∈ ℝⁿ, λᵢ ∈ ℝ: 𝐴𝐯 = λᵢ𝐯"
@test to_latex("\\sqrt\\cbrt") == "√∛"

latexstring = "\\bf{boldface} \\it{italic} \\bb{blackboard} \\cal{calligraphic} \\frak{fraktur} \\mono{monospace}"
@test to_latex(latexstring) == "𝐛𝐨𝐥𝐝𝐟𝐚𝐜𝐞 𝑖𝑡𝑎𝑙𝑖𝑐 𝕓𝕝𝕒𝕔𝕜𝕓𝕠𝕒𝕣𝕕 𝓬𝓪𝓵𝓵𝓲𝓰𝓻𝓪𝓹𝓱𝓲𝓬 𝔣𝔯𝔞𝔨𝔱𝔲𝔯 𝚖𝚘𝚗𝚘𝚜𝚙𝚊𝚌𝚎"

@test to_fraction("a-123", 392) == "ᵃ⁻¹²³⁄₃₉₂"

@test to_fraction_nl("α² ⋅ α²⁺³ ≡ α⁷", "ℝ: 𝐴𝐯 = λᵢ𝐯") == "α̲²̲ ̲⋅̲ ̲α̲²̲⁺̲³̲ ̲≡̲ ̲α̲⁷̲
ℝ: 𝐴𝐯 = λᵢ𝐯"
@test to_overline("abc") == "a̅b̅c̅"
@test to_underline("abc") == "a̲b̲c̲"

@test to_root("542") == "√5̅4̅2̅"
@test to_root(3,"542") == "∛5̅4̅2̅"
@test to_root(4,"542") == "∜5̅4̅2̅"
@test to_root(17,"542") == "¹⁷√5̅4̅2̅"
@test to_root(1, "1") == "1"
@test to_root(1, "-1") == "-1"


text:
using Compat
################################################################################
# First try to detect and load existing libmxnet
################################################################################
libmxnet_detected = false
libmxnet_curr_ver = "master"

if haskey(ENV, "MXNET_HOME")
  info("MXNET_HOME environment detected: $(ENV["MXNET_HOME"])")
  info("Trying to load existing libmxnet...")
  lib = Libdl.find_library(["libmxnet", "libmxnet.so"], ["$(ENV["MXNET_HOME"])/lib"])
  if !isempty(lib)
    info("Existing libmxnet detected at $lib, skip building...")
    libmxnet_detected = true
  else
    info("Failed to load existing libmxnet, trying to build from source...")
  end
end

using BinDeps
@BinDeps.setup
if !libmxnet_detected
  ################################################################################
  # If not found, try to build automatically using BinDeps
  ################################################################################
  if is_windows()
    info("Please follow the libmxnet documentation on how to build manually")
    info("or to install pre-build packages:")
    info("http://mxnet.readthedocs.io/en/latest/how_to/build.html#building-on-windows")
    error("Automatic building libxmnet on Windows is currently not supported yet.")
  end

  blas_path = Libdl.dlpath(Libdl.dlopen(Base.libblas_name))

  if VERSION >= v"0.5.0-dev+4338"
    blas_vendor = Base.BLAS.vendor()
  else
    blas_vendor = Base.blas_vendor()
  end

  ilp64 = ""
  if blas_vendor == :openblas64
    ilp64 = "-DINTERFACE64"
  end

  if blas_vendor == :unknown
    info("Julia is built with an unkown blas library ($blas_path).")
    info("Attempting build without reusing the blas library")
    USE_JULIA_BLAS = false
  elseif !(blas_vendor in (:openblas, :openblas64))
    info("Unsure if we can build against $blas_vendor.")
    info("Attempting build anyway.")
    USE_JULIA_BLAS = true
  else
    USE_JULIA_BLAS = true
  end

  blas_name = blas_vendor == :openblas64 ? "openblas" : string(blas_vendor)
  MSHADOW_LDFLAGS = "MSHADOW_LDFLAGS=-lm $blas_path"

  #--------------------------------------------------------------------------------
  # Build libmxnet
  mxnet = library_dependency("mxnet", aliases=["mxnet", "libmxnet", "libmxnet.so"])

  _prefix = joinpath(BinDeps.depsdir(mxnet), "usr")
  _srcdir = joinpath(BinDeps.depsdir(mxnet), "src")
  _mxdir  = joinpath(_srcdir, "mxnet")
  _libdir = joinpath(_prefix, "lib")
  # We have do eagerly delete the installed libmxnet.so
  # Otherwise we won't rebuild on an update.
  run(`rm -f $_libdir/libmxnet.so`)
  provides(BuildProcess,
    (@build_steps begin
      CreateDirectory(_srcdir)
      CreateDirectory(_libdir)
      @build_steps begin
        BinDeps.DirectoryRule(_mxdir, @build_steps begin
          ChangeDirectory(_srcdir)
          `git clone --recursive https://github.com/dmlc/mxnet`
        end)
        @build_steps begin
          ChangeDirectory(_mxdir)
          # TODO(vchuravy). We have to reset mshadow/make/mshadow.mk
          `git -C mshadow checkout -- make/mshadow.mk`
          `git fetch`
          `git checkout $libmxnet_curr_ver`
          `git submodule update`
          `sed -i -s "s/MSHADOW_CFLAGS = \(.*\)/MSHADOW_CFLAGS = \1 $ilp64/" mshadow/make/mshadow.mk`
        end
        FileRule(joinpath(_mxdir, "config.mk"), @build_steps begin
          ChangeDirectory(_mxdir)
          if is_apple()
            `cp make/osx.mk config.mk`
          else
            `cp make/config.mk config.mk`
          end
          `sed -i -s 's/USE_OPENCV = 1/USE_OPENCV = 0/' config.mk`
        end)
        @build_steps begin
          ChangeDirectory(_mxdir)
          `cp ../../cblas.h include/cblas.h`
          if USE_JULIA_BLAS
            `make -j$(nprocs()) USE_BLAS=$blas_name $MSHADOW_LDFLAGS`
          else
            `make -j$(nprocs())`
          end
        end
        FileRule(joinpath(_libdir, "libmxnet.so"), @build_steps begin
          `cp $_mxdir/lib/libmxnet.so $_libdir/`
        end)
      end
    end), mxnet, installed_libpath=_libdir)

  @BinDeps.install Dict(:mxnet => :mxnet)
end


text:
import CUDAdrv: cuArray3DCreate, CUarray, CUarray_format, CUDA_ARRAY3D_DESCRIPTOR, cuArrayDestroy
import CUDAdrv: cuMemcpyDtoA, cuMemcpyHtoA, cuMemcpy2D, CUDA_MEMCPY2D, cuMemcpy3D, CUDA_MEMCPY3D
import CuArrays: CuArray

"""
Type to handle CUDA arrays which are opaque device memory buffers optimized for texture fetching.
The only way to initialize the content of this objects is by copying from host or device arrays using the constructor or `copyto!` calls.
"""
mutable struct CuTextureArray{T,N}
    handle::CUarray
    dims::Dims{N}
    # TODO: hold the CUDA context here so that it is not finalized before this object

    function CuTextureArray{T,N}(dims::Dims{N}) where {T,N}
        Ta = cuda_texture_alias_type(T)
        _assert_alias_size(T, Ta)
        nchan, format = _alias_type_to_nchan_and_format(Ta)

        if N == 2
            width, height = dims
            depth = 0
            @assert 1 <= width "CUDA 2D array (texture) width must be >= 1"
            # @assert witdh <= CU_DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE2D_WIDTH
            @assert 1 <= height "CUDA 2D array (texture) height must be >= 1"
            # @assert height <= CU_DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE2D_HEIGHT
        elseif N == 3
            width, height, depth = dims
            @assert 1 <= width "CUDA 3D array (texture) width must be >= 1"
            # @assert witdh <= CU_DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE3D_WIDTH
            @assert 1 <= height "CUDA 3D array (texture) height must be >= 1"
            # @assert height <= CU_DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE3D_HEIGHT
            @assert 1 <= depth "CUDA 3D array (texture) depth must be >= 1"
            # @assert depth <= CU_DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE3D_DEPTH
        elseif N == 1
            width = dims[1]
            height = depth = 0
            @assert 1 <= width "CUDA 1D array (texture) width must be >= 1"
            # @assert witdh <= CU_DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE1D_WIDTH
        else
            "CUDA arrays (texture memory) can only have 1, 2 or 3 dimensions"
        end

        allocateArray_ref = Ref(CUDA_ARRAY3D_DESCRIPTOR(width, # Width::Csize_t
                                                        height, # Height::Csize_t
                                                        depth, # Depth::Csize_t
                                                        format, # Format::CUarray_format
                                                        UInt32(nchan), # NumChannels::UInt32
                                                        0))

        handle_ref = Ref{CUarray}(C_NULL)
        cuArray3DCreate(handle_ref, allocateArray_ref)

        t = new{T,N}(handle_ref[], dims)
        finalizer(unsafe_free!, t)
        return t
    end
end

function unsafe_free!(t::CuTextureArray)
    if t.handle != C_NULL
        cuArrayDestroy(t.handle)
        t.handle = C_NULL
    end
    return nothing
end

CuTextureArray{T}(n::Int) where {T} = CuTextureArray{T,1}((n,))
CuTextureArray{T}(nx::Int, ny::Int) where {T} = CuTextureArray{T,2}((nx, ny))
CuTextureArray{T}(nx::Int, ny::Int, nz::Int) where {T} = CuTextureArray{T,3}((nx, ny, nz))

Base.eltype(tm::CuTextureArray{T,N}) where {T,N} = T
Base.size(tm::CuTextureArray) = tm.dims

### Memory transfer

function Base.copyto!(dst::CuTextureArray{T,1}, src::Union{Array{T,1}, CuArray{T,1}}) where {T}
    @assert dst.dims == size(src) "CuTextureArray and CuArray sizes must match"
    if isa(src, CuArray)
        cuMemcpyDtoA(dst.handle, 0, src.ptr, dst.dims[1] * sizeof(T))
    else
        cuMemcpyHtoA(dst.handle, 0, pointer(src), dst.dims[1] * sizeof(T))
    end
    return dst
end

function Base.copyto!(dst::CuTextureArray{T,2}, src::Union{Array{T,2}, CuArray{T,2}}) where {T}
    @assert dst.dims == size(src) "CuTextureArray and source array sizes must match"
    isdevmem = isa(src, CuArray)
    copy_ref = Ref(CUDA_MEMCPY2D(0, # srcXInBytes::Csize_t
        0, # srcY::Csize_t
        isdevmem ? CUDAdrv.CU_MEMORYTYPE_DEVICE : CUDAdrv.CU_MEMORYTYPE_HOST, # srcMemoryType::CUmemorytype
        isdevmem ? 0 : pointer(src), # srcHost::Ptr{Cvoid}
        isdevmem ? pointer(src) : 0, # srcDevice::CUdeviceptr
        0, # srcArray::CUarray
        0, # srcPitch::Csize_t ### TODO: check why this cannot be `size(src.dims, 1) * sizeof(T)` as it should
        0, # dstXInBytes::Csize_t
        0, # dstY::Csize_t
        CUDAdrv.CU_MEMORYTYPE_ARRAY, # dstMemoryType::CUmemorytype
        0, # dstHost::Ptr{Cvoid}
        0, # dstDevice::CUdeviceptr
        dst.handle, # dstArray::CUarray
        0, # dstPitch::Csize_t
        dst.dims[1] * sizeof(T), # WidthInBytes::Csize_t
        dst.dims[2], # Height::Csize_t
    ))
    cuMemcpy2D(copy_ref)
    return dst
end

function Base.copyto!(dst::CuTextureArray{T,3}, src::Union{Array{T,3}, CuArray{T,3}}) where {T}
    @assert dst.dims == size(src) "CuTextureArray and source array sizes must match"
    isdevmem = isa(src, CuArray)
    copy_ref = Ref(CUDA_MEMCPY3D(0, # srcXInBytes::Csize_t
        0, # srcY::Csize_t
        0, # srcZ::Csize_t
        0, # srcLOD::Csize_t
        isdevmem ? CUDAdrv.CU_MEMORYTYPE_DEVICE : CUDAdrv.CU_MEMORYTYPE_HOST, # srcMemoryType::CUmemorytype
        isdevmem ? 0 : pointer(src), # srcHost::Ptr{Cvoid}
        isdevmem ? pointer(src) : 0, # srcDevice::CUdeviceptr
        0, # srcArray::CUarray
        0, # reserved0::Ptr{Cvoid}
        size(src, 1) * sizeof(T), # srcPitch::Csize_t
        size(src, 2), # srcHeight::Csize_t
        0, # dstXInBytes::Csize_t
        0, # dstY::Csize_t
        0, # dstZ::Csize_t
        0, # dstLOD::Csize_t
        CUDAdrv.CU_MEMORYTYPE_ARRAY, # dstMemoryType::CUmemorytype
        0, # dstHost::Ptr{Cvoid}
        0, # dstDevice::CUdeviceptr
        dst.handle, # dstArray::CUarray
        0, # reserved1::Ptr{Cvoid}
        0, # dstPitch::Csize_t
        0, # dstHeight::Csize_t
        dst.dims[1] * sizeof(T), # WidthInBytes::Csize_t
        dst.dims[2], # Height::Csize_t
        dst.dims[3], # Depth::Csize_t
    ))
    cuMemcpy3D(copy_ref)
    return dst
end

function CuTextureArray(a::Union{Array{T,N}, CuArray{T,N}}) where {T,N}
    t = CuTextureArray{T}(size(a)...)
    copyto!(t, a)
    return t
end


text:
#Bayesian_Opt_RNN_XOR.jl

#Set the random seed:

srand(1234)

include("Reccurent_XOR.jl")
include("Kernals.jl")
include("gaussian_process.jl")

# Here we will create the reccurent neural network:

Layer_1=uniform(0,1,1,2)
Layer_2=uniform(0,1,2,1)
recurrent_layer=uniform(0,1,2,2)

Seq_Len=100
epochs=100 #1000 doesnt demonstrate much as run time is longer

function hyper_curry(h)
    return (x->sigmoid(x,h))
end

function hyper_curry_deriv(h)
    return (x->sigmoid_deriv(x,h))
end

#initialise learning rates and hyper-parameters

a=0.001  #Change to 0.001
b=1
c=0.001
d=1
N=10 #10 does well

#15 gives Minimum Average Square Error for Random Selection = 0.11429806093691719
#Minimum Average Square Error for Bayesian Optimization = 0.1226313941360362

Random_Learning_Rates=uniform(a,b,N,1)

Random_Hyperparameters=uniform(c,d,N,1)

Random_Mat=cat(2,Random_Learning_Rates,Random_Hyperparameters)
Random_Average_of_Temporal_Square_Errors=zeros(N)

tic()
for i=1:length(Random_Learning_Rates)
    node_function=hyper_curry(Random_Mat[i,2])
    node_deriv=hyper_curry_deriv(Random_Mat[i,2])
    learning_rate=Random_Mat[i,1]

    #Here P is the first epoch of the temporal and the last epoch of the temporal SE Values in a 1x2 matrix
    P=Train_Reccurent_Net_Loop(epochs,Layer_1,Layer_2,recurrent_layer,learning_rate,node_function,node_deriv,Seq_Len)

    Random_Average_of_Temporal_Square_Errors[i]=(1/length(P[2]))*sum(P[2])

    #Random_MSE[i]=Train_Neural_Net_Loop(epochs,Layer_1,Layer_2,learning_rate,node_function,node_deriv)[3]
    println("Epoch Complete")
end

println("Random Learning Rates Training Completed")
totaltimerand=toc()

# Plot the random hyper-parameter and Avergae SE values:

using PyPlot
# fig = figure("pyplot_subplot_mixed",figsize=(7,7))
ax=axes()

surf(reshape(Random_Learning_Rates,size(Random_Average_of_Temporal_Square_Errors)),reshape(Random_Hyperparameters,size(Random_Average_of_Temporal_Square_Errors)),Random_Average_of_Temporal_Square_Errors,alpha=0.65,color="#40d5bb")
title("Average Temporal SE Plot for varied Sigmoid Hyper-Parameters and Learning Rates")
xlabel("Learning Rates")
ylabel("Hyper-Parameters")
zlabel("Mean Square Error")
grid("on")
show()

#Bayesian Optimization Example===================================================

#Here are the points we can pick from in the Optimization
number_of_points=75 # at 75 was 0.12035 vs 0.11 #try

LR_Test=linspace(a,b,number_of_points)
HP_Test=linspace(c,d,number_of_points)

#Here is the carteisan product of these written as a vector
Test=gen_points([LR_Test,HP_Test])[1]

#We first have to pick a random point to begin bayesian optimization:

#currently starts with the midpoint, possibly randomise this:
Bayesian_Points=[Test[Int(round(length(Test)/2))]]

#Bayesian_Points is an vector of arrays where in each array first entry is LR second entry is Hyper-Parameters:

#Define hyperparemeter functions:
node_function=hyper_curry(Bayesian_Points[1][2])
node_deriv=hyper_curry_deriv(Bayesian_Points[1][2])

#Define Learning Rate:
learning_rate=Bayesian_Points[1][1]

#Run first train before Bayesian Optimization:
P=Train_Reccurent_Net_Loop(epochs,Layer_1,Layer_2,recurrent_layer,learning_rate,node_function,node_deriv,Seq_Len)
Bayesian_Temporal_Means=[(1/length(P[2]))*sum(P[2])]

mu=zeros(number_of_points)
sigma=zeros(number_of_points)
#Begin Bayesian Optimization:

# sigma_control=linspace(1,0.1,N) #Here we have added a parameter that fades away with time

tic()
for k=2:N
    D=[(Bayesian_Points[i],Bayesian_Temporal_Means[i]) for i=1:length(Bayesian_Points)]
    mu, sigma, D=gaussian_process_chol(std_exp_square_ker,D,1e-6,Test)  #This line here has been changed to chol version of gaussian_process
    println("Gaussian Process", k, " Complete","\r")
    mu=reshape(mu,length(mu));
    sigma=reshape(sigma,length(sigma))

    new_point=findmin(mu-sigma)[2]

    #Here we will need to change the number 2 to k
    Bayesian_Points=cat(1,Bayesian_Points,[Test[new_point]])

    learning_rate=Bayesian_Points[k][1]

    node_function=hyper_curry(Bayesian_Points[k][2])

    node_deriv=hyper_curry_deriv(Bayesian_Points[k][2])

    P=Train_Reccurent_Net_Loop(epochs,Layer_1,Layer_2,recurrent_layer,learning_rate,node_function,node_deriv,Seq_Len)
    value_to_be_appended=(1/length(P[2]))*sum(P[2])

    if value_to_be_appended !=Bayesian_Temporal_Means[k-1]
        Bayesian_Temporal_Means=cat(1,Bayesian_Temporal_Means,[value_to_be_appended])
        println("Epoch ", k, " Complete")

    else
        println("Found Optimum on the ", k-1, " iteration of ", N, " iterations")
        Bayesian_Points=Bayesian_Points[1:length(Bayesian_Points)-1]
        break
    end

end
totalbayestime=toc()

println("Bayesian_Learning_Rates Training Complete")
println("total time rand = ",totaltimerand)
println("total time Bayes = ",totalbayestime)

# Bayesian Plotting =========================================================

LR=[Bayesian_Points[i][1] for i=1:length(Bayesian_Points)]
HP=[Bayesian_Points[i][2] for i=1:length(Bayesian_Points)]

#Move this to the bottom
println("Minimum Average Square Error for Random Selection = ", minimum(Random_Average_of_Temporal_Square_Errors))
println("Minimum Average Square Error for Bayesian Optimization = ", minimum(Bayesian_Temporal_Means))
println("Maximum Average Square Error for Random Selection = ", maximum(Random_Average_of_Temporal_Square_Errors))
println("Maximum Average Square Error for Bayesian Optimization = ", maximum(Bayesian_Temporal_Means))

using PyPlot
# fig = figure("pyplot_subplot_mixed",figsize=(7,7))
# ax=axes()
scatter(LR,HP,Bayesian_Temporal_Means,alpha=0.65,color="#40d5bb")
surf(LR,HP,mu+2*sigma,alpha=0.3) #This should be just mu not 2*mu
surf(LR,HP,mu-2*sigma,alpha=0.3)

title("Average Temporal SE Plot for Optimized Sigmoid Hyper-Parameters and Learning Rates")
xlabel("Learning Rates")
ylabel("Hyper-Parameters")
zlabel("Average Square Error")

grid("off")
show()

"""
============================================== Results ==========================================

_________ Experiment 1_____________________

We will see how the error changes with various values First we will change
Seq_Len=100
epochs=100
N=10
number_of_points=125

Output---->
Minimum Average Square Error for Random Selection = 0.1203530371912469
Minimum Average Square Error for Bayesian Optimization = 0.11446107828879067
Maximum Average Square Error for Random Selection = 0.13370180297876086
Maximum Average Square Error for Bayesian Optimization = 0.12899905541674905

We now change the value N to see what happens:
==============================================
N=20

Minimum Average Square Error for Random Selection = 0.11933513840710926
Minimum Average Square Error for Bayesian Optimization = 0.11560440985781831
Maximum Average Square Error for Random Selection = 0.1303499358710453
Maximum Average Square Error for Bayesian Optimization = 0.131934570268544

as we can see this has far less impact

==============================================
N=5

Minimum Average Square Error for Random Selection = 0.12435012501543988
Minimum Average Square Error for Bayesian Optimization = 0.12383116207821351
Maximum Average Square Error for Random Selection = 0.12989541500123994
Maximum Average Square Error for Bayesian Optimization = 0.12737815821942364

More random ones

Seq_Len=100
epochs=100
N=100
number_of_points=75

see how this gaussian_process
output-->

Minimum Average Square Error for Random Selection = 0.11438941310870633
Minimum Average Square Error for Bayesian Optimization = 0.11159015815940883
Maximum Average Square Error for Random Selection = 0.13087656813355802
Maximum Average Square Error for Bayesian Optimization = 0.1300878843237718

==============================================

Lets try a large sequence length

Seq_Len=3000
epochs=100
N=15
number_of_points=75

Minimum Average Square Error for Random Selection = 0.12499468626485709
Minimum Average Square Error for Bayesian Optimization = 0.12493735720265979
Maximum Average Square Error for Random Selection = 0.12669453032837658
Maximum Average Square Error for Bayesian Optimization = 0.12662348856300

==============================================

Seq_Len=100
epochs=50
N=15
number_of_points=75

Minimum Average Square Error for Random Selection = 0.1179449010241154
Minimum Average Square Error for Bayesian Optimization = 0.12358375337646431
Maximum Average Square Error for Random Selection = 0.1296091692096744
Maximum Average Square Error for Bayesian Optimization = 0.13191330829847336

"""


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "libpng"
version = v"1.6.37"

# Collection of sources required to build libpng
sources = [
    ArchiveSource("https://sourceforge.net/projects/libpng/files/libpng16/$(version)/libpng-$(version).tar.gz",
                  "daeb2620d829575513e35fecc83f0d3791a620b9b93d800b763542ece9390fb4"),
]

version = v"1.6.38" # <--- This version number is a lie, we need to bump it to build for experimental platforms

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libpng-*/
mkdir build && cd build
if [[ "${target}" == aarch64-apple-darwin* ]]; then
    # Let CMake know this platform supports NEON extension
    FLAGS=(-DPNG_ARM_NEON=on)
fi
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" \
    -DCMAKE_BUILD_TYPE=Release \
    -DPNG_STATIC=OFF \
    "${FLAGS[@]}" \
    ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libpng16", :libpng)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Zlib_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
#
# Computes the density function of a list of values
#
struct Density
  x :: Vector{Float64}
  d :: Vector{Float64}
  step :: Float64
  norm :: Int64
end

"""
```
fitdensity(x; step, norm)
```

Obtains the density function given data sampled.

Use `step=(Float64)` to control the bin step. Use `norm=(0 or 1)` to set
if the number of data points or the probability of finding a data point
within ± step/2 will be output.

By default, `norm=1` (probability) and the step is `(xmax-xmin)/100`.

# Examples
```jldoctest
julia> x = randn(1000)

julia> d = fitdensity(x)

 ------------------- Density -------------

  `d` contains the number of data points within x ± 0.028325979340904375

 -----------------------------------------

```
"""
function fitdensity(v;nbins=nothing,
                      step=nothing, steptype="absolute",
                      vmin=nothing, vmax=nothing,
                      norm :: Int64 = 1)

  ndata = length(v)

  if vmin == nothing
    vmin = minimum(v)
  end
  if vmax == nothing
    vmax = maximum(v)
  end
  if nbins == nothing
    nbins = 100
  end
  if step == nothing
    step = (vmax - vmin)/nbins
  else
    # By default, the step size is absolute
    if steptype == "relative"
      step = step*(vmax-vmin)/nbins
    elseif steptype != "absolute"
      error(" steptype must be \"relative\" or \"absolute\"")
    end
  end

  x = Vector{Float64}(undef,nbins)
  df = Vector{Float64}(undef,nbins)

  binstep = (vmax - vmin)/nbins
  for i in 1:nbins
    x[i] = vmin + (i-1)*binstep + binstep/2
    nv = 0
    for j in 1:ndata
      if ( v[j] > x[i] - step/2 ) && ( v[j] <= x[i] + step/2 )
        nv = nv + 1
      end
    end
    binsize = min(vmax,x[i]+step/2) - max(vmin,x[i]-step/2)
    if norm == 0
      df[i] = nv
    elseif norm == 1
      df[i] = nv/(binsize*ndata)
    end
  end

  return Density(x,df,step,norm)
end
export fitdensity

function Base.show( io :: IO, d :: Density )
  if d.norm == 0
    s = "number of"
  end
  if d.norm == 1
    s = "probability of finding"
  end
  println("
 ------------------- Density -------------

  d contains the $s data points within x ± $(round(d.step/2,sigdigits=3))

 ----------------------------------------- ")
end


text:
using Plots

# Model
include_model("bimanipulation_block_wall")

# Horizon
T = 26

# Time step
tf = 2.5 # 2.5
h = tf / (T-1)

# Bounds
_uu = Inf * ones(model.m)
_uu[model.idx_u[1:4]] .= [Inf; 0.0; 0.0; 0.0]
# _uu[model.idx_u[1]] = 0.5
_ul = zeros(model.m)
_ul[model.idx_u[1:4]] .= [-Inf; -0.0; 0.0; 0.0]
# _ul[model.idx_u[1]] = 0.5
ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
q1 = [0.0, r, 0.0, -r, r + 0.05, -r, r + 0.05]
x1 = [q1; q1]
# qT = [0.0, r, 0.0, -2.0 * r, r, -2.0 * r, r]

# x1 = [q1; q1]
qT = [r, r, 0.0, -0, r, -0, r]
qT = [r, r, 0.5 * π, -0, r, -0, r]
ϕ_func(model, qT)
xT = [qT; qT]

# Trajectory initialization
q0 = [q1, linear_interpolation(q1, qT, T)...]
x0 = configuration_to_state(q0)

u0 = [[0.001; 0.001 * rand(model.m-1)] for t = 1:T-1] # random controls

xl, xu = state_bounds(model, T, x1 = x1, xT = [xT[1:3]; Inf * ones(4); xT[8:10]; Inf * ones(4)])
for t = 1:T
    xl[t][5] = 0.25 * r
    xl[t][12] = 0.25 * r
end
# Objective
include_objective(["velocity", "nonlinear_stage", "control_velocity"])
obj_velocity = velocity_objective(
    [t > T / 2 ? Diagonal(1.0 * ones(model.nq)) : Diagonal(1.0 * ones(model.nq)) for t = 1:T-1],
    model.nq,
    h = h,
    idx_angle = collect([3]))

obj_tracking = quadratic_tracking_objective(
    [Diagonal(1.0 * ones(model.n)) for t = 1:T],
    [Diagonal([1.0e-1 * ones(model.nu);
        zeros(model.nc);
        zeros(model.nb);
        zeros(model.m - model.nu - model.nc - model.nb)]) for t = 1:T-1],
    [xT for t = 1:T],
    [zeros(model.m) for t = 1:T])

obj_ctrl_vel = control_velocity_objective(Diagonal([1.0e-2 * ones(model.nu); 1.0 * ones(model.nc + model.nb); zeros(model.m - model.nu - model.nc - model.nb)]))

obj_penalty = PenaltyObjective(1.0e5, model.m)

obj = MultiObjective([obj_penalty, obj_velocity, obj_tracking, obj_ctrl_vel])

# Constraints
include_constraints(["contact"])
con_contact = contact_constraints(model, T)
con = multiple_constraints([con_contact])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
@time z̄, info = solve(prob, copy(z0), tol = 1.0e-2, c_tol = 1.0e-2, max_iter=2500)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)
q = state_to_configuration(x̄)
u = [u[model.idx_u] for u in ū]
γ = [u[model.idx_λ] for u in ū]
b = [u[model.idx_b] for u in ū]
h̄ = h

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
# open(vis)
visualize!(vis, model,
    # q0, u0,
	q, u,
    Δt = h)

plot(hcat(q...)[1:2, :]', width = 1.0)
plot(hcat([u..., u[end]]...)[model.idx_u[1:4], :]', linetype = :steppost, color = :black, width = 1.0, labels = "")

model.idx_u

max_issues_repo_name: thowell/DirectMotionPlanning


text:
# RK

"""
    RK(u,Δt,r₁;[rk::RKParams=RK31])

Construct an integrator to advance a system of the form

du/dt = r₁(u,t)

The resulting integrator will advance the state `u` by one time step, `Δt`.

# Arguments

- `u` : example of state vector data
- `Δt` : time-step size
- `r₁` : operator acting on type `u` and `t` and returning `u`
"""
struct RK{NS,FR1,TU}

  # time step size
  Δt :: Float64

  rk :: RKParams

  r₁ :: FR1  # function of u and t, returns TU

  # scratch space
  qᵢ :: TU
  w :: Vector{TU}

end

function (::Type{RK})(u::TU,Δt::Float64,rhs::FR1;
                          rk::RKParams{NS}=RK31) where {TU,FR1,NS}

    # check for methods for r₁ and r₂
    if hasmethod(rhs,(TU,Float64))
        r₁ = rhs
    else
        error("No valid operator for r₁ supplied")
    end

    # scratch space
    qᵢ = deepcopy(u)
    w = [deepcopy(u) for i = 1:NS] # one extra for last step in tuple form

    # fuse the time step size into the coefficients for some cost savings
    rkdt = deepcopy(rk)
    rkdt.a .*= Δt
    rkdt.c .*= Δt

    rksys = RK{NS,typeof(r₁),TU}(Δt,rkdt,r₁,qᵢ,w)

    #pre-compile
    rksys(0.0,u)

    return rksys
end

function Base.show(io::IO, scheme::RK{NS,FR1,TU}) where {NS,FR1,TU}
    println(io, "Order-$NS RK integator with")
    println(io, "   State of type $TU")
    println(io, "   Time step size $(scheme.Δt)")
end

# Advance the RK solution by one time step
# This form works when u is a tuple of state vectors
function (scheme::RK{NS,FR1,TU})(t::Float64,u::TU) where {NS,FR1,TU <: Tuple}
  @get scheme (Δt,rk,r₁,qᵢ,w)

  # Each of the coefficients includes the time step size

  i = 1
  tᵢ₊₁ = t
  for I in eachindex(u)
    qᵢ[I] .= u[I]
  end

  if NS > 1
    # first stage, i = 1

    w[i] = r₁(u,tᵢ₊₁)
    for I in eachindex(u)
      w[i][I] .*= rk.a[i,i]  # gᵢ
      u[I] .= qᵢ[I] .+ w[i][I]
    end
    tᵢ₊₁ = t + rk.c[i]

    # stages 2 through NS-1
    for i = 2:NS-1

      w[i] = r₁(u,tᵢ₊₁)
      for I in eachindex(u)
        w[i-1][I] ./= rk.a[i-1,i-1] # w(i,i-1)

        w[i][I] .*= rk.a[i,i] # gᵢ
        u[I] .= qᵢ[I] .+ w[i][I] # r₁
        for j = 1:i-1
          u[I] .+= rk.a[i,j]*w[j][I]
        end
      end
      tᵢ₊₁ = t + rk.c[i]

    end
    i = NS
    for I in eachindex(u)
      w[i-1][I] ./= rk.a[i-1,i-1] # w(i,i-1)
    end
  end

  # final stage (assembly)
  w[i] = r₁(u,tᵢ₊₁)
  for I in eachindex(u)
    w[i][I] .*= rk.a[i,i]
    u[I] .= qᵢ[I] .+ w[i][I] # r₁
    for j = 1:i-1
      u[I] .+= rk.a[i,j]*w[j][I] # r₁
    end
  end
  t = t + rk.c[i]

  return t, u

end

# Advance the RK solution by one time step
function (scheme::RK{NS,FR1,TU})(t::Float64,u::TU) where {NS,FR1,TU}
  @get scheme (Δt,rk,r₁,qᵢ,w)

  # Each of the coefficients includes the time step size

  i = 1
  tᵢ₊₁ = t
  qᵢ .= u

  if NS > 1
    # first stage, i = 1
    w[i] .= rk.a[i,i].*r₁(u,tᵢ₊₁) # gᵢ

    u .= qᵢ .+ w[i]
    tᵢ₊₁ = t + rk.c[i]

    # stages 2 through NS-1
    for i = 2:NS-1
      w[i-1] ./= rk.a[i-1,i-1] # w(i,i-1)
      w[i] .= rk.a[i,i].*r₁(u,tᵢ₊₁) # gᵢ

      u .= qᵢ .+ w[i] # r₁
      for j = 1:i-1
        u .+= rk.a[i,j]*w[j]
      end
      tᵢ₊₁ = t + rk.c[i]

    end
    i = NS
    w[i-1] ./= rk.a[i-1,i-1] # w(i,i-1)
  end

  # final stage (assembly)
  u .= qᵢ .+ rk.a[i,i].*r₁(u,tᵢ₊₁) # r₁
  for j = 1:i-1
    u .+= rk.a[i,j]*w[j] # r₁
  end
  t = t + rk.c[i]

  return t, u

end


text:
using ModelingToolkit
fn = "test/data/lorenz.xml"
xml = readxml(fn).root
eqs = parse_node(xml)
ODESystem(eqs)
@test !isnothing() # TODO


text:
module BipolarSphericalHarmonics

using SphericalHarmonicModes
using SphericalHarmonicModes: AbstractTriangleIterator
using SphericalHarmonicArrays
using SphericalHarmonics
using SphericalHarmonics: NorthPole, SouthPole, getY
import SphericalHarmonics: eltypeY
using VectorSphericalHarmonics
using SHTOOLS
using LinearAlgebra
using StaticArrays
using WignerSymbols
using Base: @propagate_inbounds

abstract type SHType end
struct SH <:SHType end
struct VSH{VSHType, Basis} <:SHType
    YT :: VSHType
    B :: Basis
end
struct GSH <:SHType end
Base.broadcastable(S::SHType) = Ref(S)

export monopolarharmonics
export monopolarharmonics!
export biposh
export biposh!
export biposh_flippoints
export SH
export VSH
export GSH

"""
    BSHCache

Wrapper around monopolar harmonics that are coupled to evaluate the bipolar harmonics. The cache may be constructed using
[`cache`](@ref).
"""
struct BSHCache{T, ST<:SHType, S, W}
    SHT :: ST
    S1 :: S
    S2 :: S
    W3j :: W
    C :: W
end

BSHCache{T}(SHT::SHType, S1, S2, W, C) where {T} = BSHCache{T, typeof(SHT), typeof(S1), typeof(W)}(SHT, S1, S2, W, C)
flip(B::BSHCache) = typeof(B)(B.SHT, B.S2, B.S1, B.W3j, B.C)

_cache(::SH, T, lmax) = SphericalHarmonics.cache(T, lmax)
_cache(V::VSH, T, lmax) = VectorSphericalHarmonics.VSHCache(T, V.YT, V.B, ML(ZeroTo(lmax)))
_cache(::GSH, T, lmax) = VectorSphericalHarmonics.VSHCache(T, ML(ZeroTo(lmax)))

_bshtype(SHT, C, Y1, Y2) = typeof(C * _kron(Y1, Y2))

"""
    cache(SHT, T::Type, j12max)

Allocate the arrays required to evaluate the monopolar harmonics for all modes `0 ≤ j1, j2 ≤ j12max`.
`SHT` may be one of `SH()`, `GSH()`, or `VSH(YT(), B())` where `YT` and `B` are vector harmonic and basis
types offered by [`VectorSphericalHarmonics.jl`](https://github.com/jishnub/VectorSphericalHarmonics.jl).
`T` is a real type that determines the precision of the monopolar harmonics.

!!! note
    The monopolar harmonics must be initialized by calling [`monopolarharmonics!`](@ref) before the cache is used.
"""
function cache(SHT::SHType, T::Type, lmax)
    S1 = _cache(SHT, T, lmax)
    S2 = _cache(SHT, T, lmax)
    W = zeros(Cdouble, 2lmax+1)
    C = copy(W)
    TB = _bshtype(SHT, first(C), _parent(first(getY(S1))), _parent(first(getY(S2))))
    BSHCache{TB}(SHT, S1, S2, C, W)
end
cache(SHT, args...) = cache(SHT, Float64, args...)

eltypeY(::BSHCache{T}) where {T} = T

SHType(B::BSHCache) = B.SHT

_parent(A::AbstractArray) = parent(A)
_parent(A::Number) = A

_kron(Y1, Y2) = kron(Y1, Y2)
if VERSION <= v"1.7.0-DEV.1046"
    # preserve sizes in kron
    function _kron(Y1::Diagonal{<:Any, <:SVector}, Y2::Diagonal{<:Any, <:SVector})
        YD1 = diag(Y1)
        YD2 = diag(Y2)
        Diagonal(kron(YD1, YD2))
    end
end

_kronindex(ind1, ind2) = 3(ind1 + 1) + ind2 + 2
"""
    kronindex(::GSH, ind1, ind2)

Given the indices of each `GSH` vector, return the corresponding index of the bipolar spherical harmonic.
The indices `ind1` and `ind2` must lie within `-1` and `1`.
"""
function kronindex(::GSH, ind1, ind2)
    @assert -1 <= ind1 <= 1 "ind1 must satisfy -1 ≤ ind1 ≤ 1"
    @assert -1 <= ind2 <= 1 "ind2 must satisfy -1 ≤ ind2 ≤ 1"
    _kronindex(ind1, ind2)
end

_basisindrange(::Union{Polar, Cartesian}) = 1:3
_basisindrange(::Union{HelicityCovariant, SphericalCovariant}) = -1:1
_basisindrange(::VSH{<:Any,B}) where {B} = _basisindrange(B())
_shiftbasisinds(ind, ::Union{HelicityCovariant, SphericalCovariant}) = ind
_shiftbasisinds(ind, ::Union{Polar, Cartesian}) = ind - 2 # move to the range -1:1
_shiftbasisinds(ind, ::VSH{<:Any,B}) where {B} =  _shiftbasisinds(ind, B())
"""
    kronindex(::VSH, ind11, ind12, ind21, ind22)

Given the indices of each `VSH` matrix, return the corresponding index of the bipolar spherical harmonic.
Note that the indices correspond to `ℓ - j` for the `Irreducible` vector spherical harmonic ``\\mathbf{Y}^{\\ell}_{jm}(\\hat{n})``.
"""
function kronindex(Y::VSH, ind11, ind12, ind21, ind22)
    @assert ind11 ∈ _basisindrange(Y) "ind11 must satisfy ind11 ∈ $(_basisindrange(Y))"
    @assert ind21 ∈ _basisindrange(Y) "ind21 must satisfy ind21 ∈ $(_basisindrange(Y))"
    @assert -1 <= ind12 <= 1 "ind12 must satisfy -1 ≤ ind12 ≤ 1"
    @assert -1 <= ind22 <= 1 "ind22 must satisfy -1 ≤ ind22 ≤ 1"
    ind11t = _shiftbasisinds(ind11, Y)
    ind21t = _shiftbasisinds(ind21, Y)
    CartesianIndex(_kronindex(ind11t, ind21t), _kronindex(ind12, ind22))
end

function Broadcast.broadcasted(::Broadcast.DefaultArrayStyle{1}, kronindex, ::Base.Ref{GSH}, i1::Integer, r2::AbstractUnitRange{<:Integer})
    ind1 = kronindex(GSH(), i1, first(r2))
    ind2 = kronindex(GSH(), i1, last(r2))
    ind1:ind2
end
function Broadcast.broadcasted(::Broadcast.DefaultArrayStyle{1}, kronindex, ::Base.Ref{GSH}, r1::AbstractUnitRange{<:Integer}, r2::Union{Integer, AbstractUnitRange{<:Integer}})
    ind1 = kronindex(GSH(), first(r1), first(r2))
    ind2 = kronindex(GSH(), last(r1), last(r2))
    s = length(r1) == 1 ? 1 : length(r2) == 1 ? 3 : 4
    ind1:s:ind2
end

_containertype(::Type{T}, ::Any) where {T} = Vector{T}
_containertype(::Type{T}, ::LM{SingleValuedRange, SingleValuedRange}) where {T} = MVector{1,T}
_container(el, modes) = [el for i = 1:length(modes)]
_container(el, modes::LM{SingleValuedRange, SingleValuedRange}) = MVector{1,typeof(el)}((el,))

_zerostype(::Type{T}, modes) where {T} = SHArray{T, 1, _containertype(T, modes), Tuple{typeof(modes)}}
_zeros(::Type{T}, modes) where {T} = SHArray(_container(zero(T), modes), modes)

_maybestripwrapper(T, ZT, ::Integer, ::Integer) = T
_maybestripwrapper(T, ZT, jm...) = ZT

_j1j2(::BSHCache{<:Any, SH}, j1, j2) = (j1,), (j2,)
_j1j2(::BSHCache, j1, j2) = (ML(ZeroTo(j1)),), (ML(ZeroTo(j2)),)
_j1j2(::BSHCache) = (), ()
_j1j2(B::BSHCache, j1j2modes) = _j1j2(B, _j12max(j1j2modes)...)

_j12max(j1j2modes) = maximum(first, j1j2modes), maximum(last, j1j2modes)

"""
    monopolarharmonics!(B::BSHCache, θ1, ϕ1, θ2, ϕ2, j1, j2)

Update `B` with the monopolar harmonics ``Y_{j_1,m_1}(\\theta_1, \\phi_1)`` and ``Y_{j_2,m_2}(\\theta_2, \\phi_2)``,
where ``Y`` may either be a scalar or a vector harmonic depending on `B`.
"""
function monopolarharmonics!(B::BSHCache{<:Any, SH}, θ1, ϕ1, θ2, ϕ2, j1j2...)
    j1tup, j2tup = _j1j2(B, j1j2...)
    computePlmcostheta!(B.S1, θ1, j1tup...)
    computeYlm!(B.S1, θ1, ϕ1, j1tup...)
    computePlmcostheta!(B.S2, θ2, j2tup...)
    computeYlm!(B.S2, θ2, ϕ2, j2tup...)
    return B
end

function monopolarharmonics!(B::BSHCache{<:Any, GSH}, θ1, ϕ1, θ2, ϕ2, j1j2...)
    j1tup, j2tup = _j1j2(B, j1j2...)
    VectorSphericalHarmonics.genspharm!(B.S1, j1tup..., θ1, ϕ1)
    VectorSphericalHarmonics.genspharm!(B.S2, j2tup..., θ2, ϕ2)
    return B
end

function monopolarharmonics!(B::BSHCache{<:Any, <:VSH}, θ1, ϕ1, θ2, ϕ2, j1j2...)
    j1tup, j2tup = _j1j2(B, j1j2...)
    VectorSphericalHarmonics.vshbasis!(B.S1, B.SHT.YT, B.SHT.B, j1tup..., θ1, ϕ1)
    VectorSphericalHarmonics.vshbasis!(B.S2, B.SHT.YT, B.SHT.B, j2tup..., θ2, ϕ2)
    return B
end

"""
    monopolarharmonics(SHT, θ1, ϕ1, θ2, ϕ2, j1, j2)

Evaluate the monopolar harmonics ``Y_{j_1,m_1}(\\theta_1, \\phi_1)`` and ``Y_{j_2,m_2}(\\theta_2, \\phi_2)``,
where ``Y`` may either be a scalar or a vector harmonic depending on `SHT`, which may be one of `SH()`, `GSH()`,
or `VSH(YT(), B())` where `YT` and `B` are vector harmonic and basis types offered by
[`VectorSphericalHarmonics.jl`](https://github.com/jishnub/VectorSphericalHarmonics.jl).
"""
function monopolarharmonics(SHT::SHType, θ1, ϕ1, θ2, ϕ2, j1, j2)
    T = float(mapreduce(typeof, promote_type, (θ1, ϕ1, θ2, ϕ2)))
    B = cache(SHT, T, max(j1, j2))
    monopolarharmonics!(B, θ1, ϕ1, θ2, ϕ2, j1, j2)
    return B
end

function monopolarharmonics(SHT::SHType, θ1, ϕ1, θ2, ϕ2, j1j2modes...)
    monopolarharmonics(SHT, θ1, ϕ1, θ2, ϕ2, _j12max(j1j2modes...)...)
end

monopolarharmonics(B::BSHCache) = getY(B.S1), getY(B.S2)

function degreerange(j1, j2, m::Integer)
    jmin = max(abs(j1 - j2), abs(m))
    jmax = j1 + j2
    jmin:jmax
end
degreerange(j1, j2, m::AbstractUnitRange) = degreerange(j1, j2, minimum(abs, m))
degreerange(j1, j2) = degreerange(j1, j2, 0)

neg1pow(x::Integer) = isodd(x) ? -1 : 1

_maybewrapj1j2(M, modes::AbstractTriangleIterator) = SHArray(M, modes)
_maybewrapj1j2(M, modes) = M

_modes(::Colon, ::Colon, j1, j2) = LM(degreerange(j1, j2))
_modes(::Colon, m::Union{Integer, AbstractUnitRange{<:Integer}}, j1, j2) = LM(degreerange(j1, j2, m), m)
_modes(j::Union{Integer, AbstractUnitRange{<:Integer}}, ::Colon, j1, j2) = LM(j)
_modes(j::Union{Integer, AbstractUnitRange{<:Integer}}, m::Union{Integer, AbstractUnitRange{<:Integer}}, j1, j2) = LM(j, m)
function _modes(j, m, j1j2modes)
    lmax = max(maximum(first, j1j2modes), maximum(last, j1j2modes))
    _modes(j, m, lmax, lmax)
end
_modes(jm::LM, j1j2modes) = jm

_nonnegative(r::AbstractUnitRange) = max(minimum(r), 0):maximum(r)
_negative(r::AbstractUnitRange) = minimum(r):min(-1, maximum(r))

_maybeonly(M, jm...) = M
_maybeonly(M::SHArray{<:Any,1,<:Any,Tuple{LM{SingleValuedRange, SingleValuedRange}}}, j::Integer, m::Integer) = M[1]

_modes_intersect(jm, j1, j2) = intersect(jm, LM(degreerange(j1, j2)))
_modes_intersect(j, m, j1, j2) = intersect(_modes(j, m, j1, j2), LM(degreerange(j1, j2)))

const perminds = SVector{9}((1,4,7,2,5,8,3,6,9))
const GSHPerm = CartesianIndex.(perminds)
const VSHPerm = CartesianIndex.(perminds, perminds')

_permutation(::VSH) = VSHPerm
_permutation(::GSH) = GSHPerm

function _permuteterms!(Y21, indj1j2, ::Nothing, SHT, j1, j2, jm)
    @inbounds Y21[indj1j2] = nothing
    return nothing
end
function _permuteterms!(Y21, indj1j2, Y12_j2j1::SHArray, SHT, j1, j2, jm)
    _Y21 = similar(Y12_j2j1)
    @inbounds begin
        for (indjm, (j,m)) in enumerate(first(SphericalHarmonicArrays.modes(Y12_j2j1)))
            phase = (-1)^(j1 + j2 - j)
            _Y21[indjm] = _permuteterms(Y12_j2j1[indjm], SHT) * phase
        end
        Y21[indj1j2] = _Y21
    end
    return nothing
end
function _permuteterms!(Y21, indj1j2, Y12_j2j1::Union{SArray,Diagonal{<:Any,<:SVector},Number}, SHT, j1, j2, jm::LM{SingleValuedRange, SingleValuedRange})
    j, _ = first(jm)
    phase = (-1)^(j1 + j2 - j)
    @inbounds Y21[indj1j2] = _permuteterms(Y12_j2j1, SHT) * phase
    return nothing
end
_permuteterms(Y12_j2j1_jm, SHT) = @inbounds Y12_j2j1_jm[_permutation(SHT)]
function _permuteterms(Y12_j2j1_jm::Diagonal{<:Any, <:SVector}, SHT)
    d = diag(Y12_j2j1_jm)
    v = typeof(d)(d[GSHPerm])
    Diagonal(v)
end
_permuteterms(Y12_j2j1_jm::Number, SHT) = Y12_j2j1_jm

"""
    biposh_flippoints(SHT, θ1, ϕ1, θ2, ϕ2, j, m, j1, j2)

Evaluate the bipolar harmonics ``Y_{jm}^{j_1, j_2}((\\theta_1, \\phi_1), (\\theta_2, \\phi_2))`` and
``Y_{jm}^{j_1, j_2}((\\theta_2, \\phi_2), (\\theta_1, \\phi_1))``
in one pass, utilizing symmetries of the Clebsch-Gordean coefficients.
"""
biposh_flippoints

function _biposh_flippoints(SHT::SHType, θ1, ϕ1, θ2, ϕ2, j1j2modes, jm...)
    B12 = monopolarharmonics(SHT, θ1, ϕ1, θ2, ϕ2, j1j2modes...);
    B21 = monopolarharmonics(SHT, θ2, ϕ2, θ1, ϕ1, j1j2modes...);
    biposh_flippoints(B12, B21, θ1, ϕ1, θ2, ϕ2, jm..., j1j2modes...)
end
biposh_flippoints(SHT::SHType, θ1, ϕ1, θ2, ϕ2, j, m, j1j2modes...) =
    _biposh_flippoints(SHT, θ1, ϕ1, θ2, ϕ2, j1j2modes, j, m)
biposh_flippoints(SHT::SHType, θ1, ϕ1, θ2, ϕ2, jm::LM, j1j2modes...) =
    _biposh_flippoints(SHT, θ1, ϕ1, θ2, ϕ2, j1j2modes, jm)

function _biposh_flippoints_j1j2(B12, B21, θ1, ϕ1, θ2, ϕ2, j1, j2, jm...)
    Y12 = biposh(B12, θ1, ϕ1, θ2, ϕ2, jm..., j1, j2)
    Y21 = biposh(B21, θ2, ϕ2, θ1, ϕ1, jm..., j1, j2)
    (Y12 === nothing || Y21 === nothing) && return nothing
    return Y12, Y21
end
biposh_flippoints(B12::BSHCache, B21::BSHCache, θ1, ϕ1, θ2, ϕ2, j, m, j1, j2) =
    _biposh_flippoints_j1j2(B12, B21, θ1, ϕ1, θ2, ϕ2, j1, j2, j, m)
biposh_flippoints(B12::BSHCache, B21::BSHCache, θ1, ϕ1, θ2, ϕ2, jm::LM, j1, j2) =
    _biposh_flippoints_j1j2(B12, B21, θ1, ϕ1, θ2, ϕ2, j1, j2, jm)

@inline function _biposh_flippoints_j1j2modes(B12::BSHCache{T}, B21::BSHCache{T}, θ1, ϕ1, θ2, ϕ2, j1j2modes, jm...) where {T}
    jm_modes = _modes(jm..., j1j2modes)
    TM = _maybestripwrapper(T, _zerostype(T, jm_modes), jm...)
    Y12 = Vector{Union{TM, Nothing}}(undef, length(j1j2modes))
    Y21 = similar(Y12)
    for (indj1j2, (j1, j2)) in enumerate(j1j2modes)
        if (j2,j1) in j1j2modes && modeindex(j1j2modes, j2, j1) < indj1j2
            indj2j1 = modeindex(j1j2modes, j2, j1)
            # in this case the bipolar harmonics for (j2,j1) have been computed
            # we may use these to compute the ones for (j1,j2)
            _permuteterms!(Y21, indj1j2, Y12[indj2j1], SHType(B12), j1, j2, jm_modes)
            _permuteterms!(Y12, indj1j2, Y21[indj2j1], SHType(B21), j1, j2, jm_modes)
        else
            Y12[indj1j2] = biposh(B12, θ1, ϕ1, θ2, ϕ2, jm..., j1, j2)
            Y21[indj1j2] = biposh(B21, θ2, ϕ2, θ1, ϕ1, jm..., j1, j2)
        end
    end
    _maybewrapj1j2(Y12, j1j2modes), _maybewrapj1j2(Y21, j1j2modes)
end
biposh_flippoints(B12::BSHCache, B21::BSHCache, θ1, ϕ1, θ2, ϕ2, j, m, j1j2modes) =
    _biposh_flippoints_j1j2modes(B12, B21, θ1, ϕ1, θ2, ϕ2, j1j2modes, j, m)
biposh_flippoints(B12::BSHCache, B21::BSHCache, θ1, ϕ1, θ2, ϕ2, jm::LM, j1j2modes) =
    _biposh_flippoints_j1j2modes(B12, B21, θ1, ϕ1, θ2, ϕ2, j1j2modes, jm)

"""
    biposh(SHT, θ1, ϕ1, θ2, ϕ2, j, m, j1, j2)

Evaluate the bipolar harmonic ``Y_{jm}^{j_1, j_2}((\\theta_1, \\phi_1), (\\theta_2, \\phi_2))``. `SHT` may be one of `SH()`,
`GSH()` and `VSH(YT(), B())` where `YT` and `B` are vector harmonic and basis types offered by
[`VectorSphericalHarmonics.jl`](https://github.com/jishnub/VectorSphericalHarmonics.jl).

    biposh(SHT, θ1, ϕ1, θ2, ϕ2, j, m, j1j2modes)

Evaluate the bipolar harmonic ``Y_{jm}^{j_1, j_2}((\\theta_1, \\phi_1), (\\theta_2, \\phi_2))`` for all `(j1,j2)` in
`j1j2modes`.

    biposh(SHT, θ1, ϕ1, θ2, ϕ2, jmcoll::SphericalHarmonicModes.LM, j1j2...)

Evaluate the bipolar harmonics ``Y_{jm}^{j_1, j_2}((\\theta_1, \\phi_1), (\\theta_2, \\phi_2))`` for all `(j,m)` in `jmcoll`,
where `j1j2` may be either a 2-Tuple `(j1,j2)` or a collection of 2-Tuples.
"""
biposh

function _biposh(SHT::SHType, θ1, ϕ1, θ2, ϕ2, j1j2modes, jm...)
    B = monopolarharmonics(SHT, θ1, ϕ1, θ2, ϕ2, j1j2modes...)
    biposh(B, θ1, ϕ1, θ2, ϕ2, jm..., j1j2modes...)
end
biposh(SHT::SHType, θ1, ϕ1, θ2, ϕ2, j, m, j1j2modes...) = _biposh(SHT, θ1, ϕ1, θ2, ϕ2, j1j2modes, j, m)
biposh(SHT::SHType, θ1, ϕ1, θ2, ϕ2, jm::LM, j1j2modes...) = _biposh(SHT, θ1, ϕ1, θ2, ϕ2, j1j2modes, jm)

@inline function _biposh_j1j2modes(B::BSHCache{T}, θ1, ϕ1, θ2, ϕ2, j1j2modes, jm...) where {T}
    jm_modes = _modes(jm..., j1j2modes)
    TM = _maybestripwrapper(T, _zerostype(T, jm_modes), jm...)
    M = Vector{Union{TM, Nothing}}(undef, length(j1j2modes))
    for (indj1j2, (j1, j2)) in enumerate(j1j2modes)
        M[indj1j2] = biposh(B, θ1, ϕ1, θ2, ϕ2, jm..., j1, j2)
    end
    _maybewrapj1j2(M, j1j2modes)
end
biposh(B::BSHCache, θ1, ϕ1, θ2, ϕ2, j, m, j1j2modes) = _biposh_j1j2modes(B, θ1, ϕ1, θ2, ϕ2, j1j2modes, j, m)
biposh(B::BSHCache, θ1, ϕ1, θ2, ϕ2, jm::LM, j1j2modes) = _biposh_j1j2modes(B, θ1, ϕ1, θ2, ϕ2, j1j2modes, jm)

@inline function _biposh_j1j2(B, θ1, ϕ1, θ2, ϕ2, j1, j2, jm...)
    jm_filt = _modes_intersect(jm..., j1, j2)
    jm_filt === nothing && return nothing
    Y = _zeros(eltypeY(B), jm_filt)
    biposh!(parent(Y), B, θ1, ϕ1, θ2, ϕ2, jm..., j1, j2)
    return _maybeonly(Y, jm...)
end
biposh(B::BSHCache, θ1, ϕ1, θ2, ϕ2, j, m, j1, j2) = _biposh_j1j2(B, θ1, ϕ1, θ2, ϕ2, j1, j2, j, m)
biposh(B::BSHCache, θ1, ϕ1, θ2, ϕ2, jm::LM, j1, j2) = _biposh_j1j2(B, θ1, ϕ1, θ2, ϕ2, j1, j2, jm)

const BasisPhase = SVector{9,Int}([(-1)^(α1 + α2) for α1 in -1:1, α2 in -1:1])
const VSHIndexPhase = BasisPhase'
const BasisPhase_VSHIndexPhase = BasisPhase*VSHIndexPhase

_reverse(Y::SArray; dims = :) = _reverse(Y, dims)
# special case this to use the reverse(::StaticArray) method defined by StaticArrays
_reverse(Y::SArray, ::Colon) = reverse(Y)
_reverse(Y::SArray, dims) = reverse(Y, dims = dims)
_reverse(Y::Diagonal{<:Any,<:SVector}) = Diagonal(reverse(diag(Y)))

_conjphase(Y, phase, ::SH) = conj(Y) .* phase
_conjphase(Y, phase, ::GSH) = conj(_reverse(Y)) .* BasisPhase .* phase
_conjphase(Y, phase, ::VSH{PB, SphericalCovariant}) = conj(_reverse(Y)) .* BasisPhase .* phase
function _conjphase(Y, phase, ::VSH{PB, HelicityCovariant})
    D = conj(_reverse(Y))
    Diagonal(diag(D) .* BasisPhase .* phase)
end
_conjphase(Y, phase, ::VSH{PB, <:Union{Cartesian, Polar}}) = (A = conj!(_reverse(Y, dims = 2)); A .*= phase; A)
function _conjphase(Y, phase, ::VSH{<:Any, <:Union{SphericalCovariant, HelicityCovariant}})
    A = conj!(_reverse(Y, dims = 1))
    A .*= BasisPhase_VSHIndexPhase * phase
    A
end
_conjphase(Y, phase, ::VSH{<:Any, <:Union{Cartesian, Polar}}) = conj(Y) .* VSHIndexPhase .* phase

"""
    biposh!(Y12::AbstractVetor, B::BSHCache, θ1, ϕ1, θ2, ϕ2, j, m, j1, j2)

Evaluate the bipolar harmonic ``Y_{jm}^{j_1, j_2}((\\theta_1, \\phi_1), (\\theta_2, \\phi_2))`` and store the result in `Y12`.
The cache `B` determines the type of harmonic evaluated.
"""
biposh!

biposh!(Y12::AbstractVector{T}, B::BSHCache{T}, θ1, ϕ1, θ2, ϕ2, j, m, j1::Integer, j2::Integer) where {T} =
    biposh!(Y12, B, θ1, ϕ1, θ2, ϕ2, _modes(j, m, j1, j2), j1, j2)

function biposh!(Y12::AbstractVector{T}, B::BSHCache{T}, θ1, ϕ1, θ2, ϕ2, jm::LM, j1::Integer, j2::Integer) where {T}
    jm_filt = _modes_intersect(jm, j1, j2)
    (jm_filt === nothing || isempty(jm_filt)) && return Y12
    length(Y12) >= length(jm_filt) || throw(ArgumentError("Y must have a a length >= length($jm_filt)"))
    mrange = m_range(jm_filt)

    Y1, Y2 = monopolarharmonics(B)

    j1 in l_range(first(SphericalHarmonicArrays.modes(Y1))) || throw("please compute monopolar harmonics for for j1 = $j1 first")
    j2 in l_range(first(SphericalHarmonicArrays.modes(Y2))) || throw("please compute monopolar harmonics for for j2 = $j2 first")

    Y = @view Y12[begin .+ (0:length(jm_filt)-1)]
    Y .= zero.(Y)

    for m in _nonnegative(mrange)
        jrange_m = l_range(jm_filt, m)
        modes_m = LM(jrange_m, m)
        modeindrange = modeindex(jm_filt, modes_m)
        Y12_m = @view Y12[modeindrange]
        biposh!(Y12_m, B, θ1, ϕ1, θ2, ϕ2, modes_m, j1, j2)
    end

    YS = SHArray(Y, jm_filt)
    for m in _negative(mrange)
        if -m in mrange
            jrange_m = l_range(jm_filt, m)
            jmin = minimum(jrange_m)
            # overall phase of (-1)^(j1 + j2 - j + m) appears in all terms
            overallphase = neg1pow(j1 + j2 + m - jmin)

            for j in jrange_m
                Yjminm = YS[(j,-m)]
                YS[(j,m)] = _conjphase(Yjminm, overallphase, SHType(B))
                overallphase *= -1
            end
        else
            jrange_m = l_range(jm_filt, m)
            modes_m = LM(jrange_m, m)
            modeindrange = modeindex(jm_filt, modes_m)
            Y12_m = @view Y12[modeindrange]
            biposh!(Y12_m, B, θ1, ϕ1, θ2, ϕ2, modes_m, j1, j2)
        end
    end

    return Y12
end

@inline function C_view_j1j2m_j(C, j1, j2, m, jrange)
    jmin_valid = minimum(degreerange(j1, j2, m))
    jmin = minimum(jrange)
    jmin_offset = jmin - jmin_valid
    Cv = @view C[(begin + jmin_offset) .+ (0:length(jrange)-1)]
end

function biposh!(Y12::AbstractVector{T}, B::BSHCache{T}, θ1, ϕ1, θ2, ϕ2, jm::LM{<:Any,SingleValuedRange}, j1::Integer, j2::Integer) where {T}
    jm_filt = _modes_intersect(jm, j1, j2)
    isempty(jm_filt) && return Y12
    jrange = l_range(jm_filt)
    length(Y12) >= length(jrange) || throw(ArgumentError("vector must have at least $(length(jrange)) elements, received $(length(Y12))"))
    Y = @view Y12[begin .+ (0:length(jrange)-1)]
    Y .= zero.(Y)

    Y1, Y2 = monopolarharmonics(B)

    j1 in l_range(first(SphericalHarmonicArrays.modes(Y1))) || throw("please compute monopolar harmonics for for j1 = $j1 first")
    j2 in l_range(first(SphericalHarmonicArrays.modes(Y2))) || throw("please compute monopolar harmonics for for j2 = $j2 first")

    m = only(m_range(jm_filt))
    Cv = C_view_j1j2m_j(B.C, j1, j2, m, jrange)

    for m1 in -j1:j1
        m2 = m - m1
        abs(m2) > j2 && continue
        clebschgordan!(B.C, B.W3j, j1, m1, j2, m)
        Y1_j1m1 = _parent(Y1[(j1,m1)])
        Y2_j2m2 = _parent(Y2[(j2,m2)])
        Y1Y2 = _kron(Y1_j1m1, Y2_j2m2)
        for ind in eachindex(Y, Cv)
            Y[ind] += Cv[ind] * Y1Y2
        end
    end

    return Y12
end

# All j for one m for (j1,j2) at the NorthPole
function biposh!(Y12::AbstractVector{T}, B::BSHCache{T,SH}, θ1::NorthPole, ϕ1, θ2, ϕ2, jm::LM{<:Any,SingleValuedRange}, j1::Integer, j2::Integer) where {T}
    jm_filt = _modes_intersect(jm, j1, j2)
    isempty(jm_filt) && return Y12
    jrange = l_range(jm_filt)
    length(Y12) >= length(jrange) || throw(ArgumentError("vector must have at least $(length(jrange)) elements, received $(length(Y12))"))
    Y12_section = @view Y12[begin .+ (0:length(jrange)-1)]

    m = only(m_range(jm_filt))
    if abs(m) > j2
        # in this case the value is known to be zero, as the spherical harmonic Ylm(0,ϕ1) == 0 for m != 0
        fill!(Y12_section, zero(eltype(Y12_section)))
        return Y12
    end
    _, Y2 = monopolarharmonics(B)
    clebschgordan!(B.C, B.W3j, j1, 0, j2, m)
    Cv = C_view_j1j2m_j(B.C, j1, j2, m, jrange)
    Y2_j2m = Y2[(j2,m)]
    Y1_l0_NP = √((2j1+1)/4pi)
    @. Y12_section = Cv * kron(Y1_l0_NP, Y2_j2m)
    return Y12
end

function biposh!(Y12::AbstractVector{T}, B::BSHCache{T,SH}, θ1::NorthPole, ϕ1, θ2::NorthPole, ϕ2, jm::LM{<:Any,SingleValuedRange}, j1::Integer, j2::Integer) where {T}
    jm_filt = _modes_intersect(jm, j1, j2)
    isempty(jm_filt) && return Y12
    jrange = l_range(jm_filt)
    length(Y12) >= length(jrange) || throw(ArgumentError("vector must have at least $(length(jrange)) elements, received $(length(Y12))"))
    Y12_section = @view Y12[begin .+ (0:length(jrange)-1)]

    m = only(m_range(jm_filt))
    if m != 0
        # in this case the value is known to be zero, as the spherical harmonic Ylm(0,ϕ1) == 0 for m != 0
        fill!(Y12_section, zero(eltype(Y12_section)))
        return Y12
    end
    clebschgordan!(B.C, B.W3j, j1, 0, j2, 0)
    Cv = C_view_j1j2m_j(B.C, j1, j2, m, jrange)
    Y1_j10_NP = √((2j1+1)/4pi)
    Y2_j20_NP = √((2j2+1)/4pi)
    @. Y12_section = Cv * kron(Y1_j10_NP, Y2_j20_NP)
    return Y12
end

function biposh!(Y12::AbstractVector{T}, B::BSHCache{T}, θ1, ϕ1, θ2::NorthPole, ϕ2, jm::LM{<:Any,SingleValuedRange}, j1::Integer, j2::Integer) where {T}
    biposh!(Y12, flip(B), θ2, ϕ2, θ1, ϕ1, jm, j2, j1)
    jm_filt = _modes_intersect(jm, j1, j2)
    Y12_section = @view Y12[begin .+ (0:length(jm_filt)-1)]
    Y = SHArray(Y12_section, jm_filt)

    @inbounds for m in m_range(jm_filt)
        jrange = l_range(jm_filt, m)
        jmin = minimum(jrange)
        phase = neg1pow(j1 + j2 - jmin)
        for j in jrange
            Y[(j,m)] *= phase
            phase *= -1
        end
    end
    return Y12
end

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

@inline function clebschgordan!(C::AbstractVector, w3j::AbstractVector, j1::Integer, m1::Integer, j2::Integer, m::Integer;
        primefactorization_cutoff = 160)

    m2 = m - m1
    jrange = max(abs(j1-j2), abs(m)):j1+j2
    @assert length(C) >= length(jrange) "Clebsch-Gordan vector does not contain enough elements"

    if abs(j1) <= primefactorization_cutoff && abs(j2) <= primefactorization_cutoff
        # Evaluate the Wigner 3j symbol
        #   /  j  j1  j2 \
        #   \ -m  m1  m2 /
        SHTOOLS.Wigner3j!(w3j, j1, j2, -m, m1, m2)

        phase = neg1pow(j1 - j2 + m)
        jrange_SHTOOLS = intersect(jrange, 0:primefactorization_cutoff)
        jrange_WigSym = primefactorization_cutoff+1:maximum(jrange)
        @inbounds for (ind, j) in enumerate(jrange_SHTOOLS)
            C[begin - 1 + ind] = w3j[begin - 1 + ind] * √(2j+1) * phase
        end
        if !isempty(jrange_WigSym)
            indshift = length(jrange_SHTOOLS)
            @inbounds for (indWS, j) in enumerate(jrange_WigSym)
                ind = indWS + indshift
                C[begin - 1 + ind] = clebschgordan(eltype(C), j1, m1, j2, m2, j, m)
            end
        end
    else
        @inbounds for (ind, j) in enumerate(jrange)
            C[begin - 1 + ind] = clebschgordan(eltype(C), j1, m1, j2, m2, j, m)
        end
    end

    return C
end

include("precompile.jl")

end


text:
type DirtyCall <: CallType end
type CleanCall <: CallType end

type Price{T <: CallType}
  amount::Float64
  callType::T
end

type Callability{T <: OptionType}
  price::Price
  optionType::T
  date::Date
end

typealias CallabilitySchedule Vector{Callability}

function has_occurred(callability::Callability, ref_date::Date, include_ref_date::Bool = true)
  # will need to expand this
  if ref_date < callability.date || (ref_date == callability.date && include_ref_date)
    return false
  else
    return true
  end
end


text:
@testset "shapes" begin
    inputsize = 20
    poollength = 10
    datalength = 100
    x = rand(Float32, inputsize, poollength, 1, datalength)
    m = DSANet(inputsize, poollength, 3, 3, 4, 1, 3, 2)
    @test size(m(x)) == (1, datalength)
    if Flux.CUDA.functional()
        @test size(gpu(m)(gpu(x))) == (1, datalength)
    end
end

@testset "misc" begin
    @test repr(DSANet(10, 10, 3, 3, 4, 1, 3, 2)) == "DSANet(10, 10, 3, 3, 4, 1, 3, 2)"
    @test repr(FluxArchitectures.Global_SelfAttn(10, 10, 3, 3, 4, 1, 3, 2)) == "Global_SelfAttn(10, 10, 3, 3, 4, 1, 3, 2, 0.1)"
    @test repr(FluxArchitectures.Local_SelfAttn(10, 10, 3, 3, 4, 4, 3, 3, 2)) == "Local_SelfAttn(10, 10, 3, 3, 4, 4, 3, 3, 2, 0.1)"
    @test_throws ArgumentError FluxArchitectures.Local_SelfAttn(10, 10, 3, 3, 4, 1, 3, 3, 2)
    @test_throws ArgumentError DSANet(10, 10, 3, 3, 4, 1, 3, 2, 20)
end

max_stars_repo_name: sdobber/FluxArchitectures

max_issues_repo_name: sdobber/FluxArchitectures

max_forks_repo_name: sdobber/FluxArchitectures