text:
using Turing, Random
using Turing: Selector, reconstruct, invlink, CACHERESET,
    SampleFromPrior, Sampler, runmodel!, SampleFromUniform
using Turing.RandomVariables
using Turing.RandomVariables: uid, _getidcs,
    set_retained_vns_del_by_spl!, is_flagged,
    set_flag!, unset_flag!, VarInfo, TypedVarInfo,
    getlogp, setlogp!, resetlogp!, acclogp!
using Distributions
using ForwardDiff: Dual
using Test

i, j, k = 1, 2, 3

dir = splitdir(splitdir(pathof(Turing))[1])[1]
include(dir*"/test/test_utils/AllUtils.jl")

@testset "RandomVariables.jl" begin
    @turing_testset "TypedVarInfo" begin
        @model gdemo(x, y) = begin
            s ~ InverseGamma(2,3)
            m ~ TruncatedNormal(0.0,sqrt(s),0.0,2.0)
            x ~ Normal(m, sqrt(s))
            y ~ Normal(m, sqrt(s))
        end
        model = gdemo(1.0, 2.0)

        vi = VarInfo()
        model(vi, SampleFromUniform())
        tvi = TypedVarInfo(vi)

        meta = vi.metadata
        for f in fieldnames(typeof(tvi.metadata))
            fmeta = getfield(tvi.metadata, f)
            for vn in fmeta.vns
                @test tvi[vn] == vi[vn]
                ind = meta.idcs[vn]
                tind = fmeta.idcs[vn]
                @test meta.dists[ind] == fmeta.dists[tind]
                @test meta.orders[ind] == fmeta.orders[tind]
                @test meta.gids[ind] == fmeta.gids[tind]
                for flag in keys(meta.flags)
                    @test meta.flags[flag][ind] == fmeta.flags[flag][tind]
                end
                range = meta.ranges[ind]
                trange = fmeta.ranges[tind]
                @test all(meta.vals[range] .== fmeta.vals[trange])
            end
        end
    end
    @turing_testset "Base" begin
        # Test Base functions:
        #   string, Symbol, ==, hash, in, keys, haskey, isempty, push!, empty!,
        #   getindex, setindex!, getproperty, setproperty!
        csym = gensym()
        vn1 = VarName(csym, :x, "[1][2]", 1)
        @test string(vn1) == "{$csym,x[1][2]}:1"
        @test string(vn1, all=false) == "x[1][2]"
        @test Symbol(vn1) == Symbol("x[1][2]")

        vn2 = VarName(csym, :x, "[1][2]", 1)
        @test vn2 == vn1
        @test hash(vn2) == hash(vn1)
        @test in(vn1, (:x,))

        function test_base!(vi)
            empty!(vi)
            @test vi.logp == 0
            @test vi.num_produce == 0

            vn = VarName(gensym(), :x, "", 1)
            dist = Normal(0, 1)
            r = rand(dist)
            gid = Selector()

            @test isempty(vi)
            @test ~haskey(vi, vn)
            push!(vi, vn, r, dist, gid)
            @test ~isempty(vi)
            @test haskey(vi, vn)

            @test length(vi[vn]) == 1
            @test length(vi[SampleFromPrior()]) == 1

            @test vi[vn] == r
            @test vi[SampleFromPrior()][1] == r
            vi[vn] = [2*r]
            @test vi[vn] == 2*r
            @test vi[SampleFromPrior()][1] == 2*r
            vi[SampleFromPrior()] = [3*r]
            @test vi[vn] == 3*r
            @test vi[SampleFromPrior()][1] == 3*r

            empty!(vi)
            @test isempty(vi)
            push!(vi, vn, r, dist, gid)

            function test_in()
                space = (:x, :y, :(z[1]))
                vn1 = genvn(:x)
                vn2 = genvn(:y)
                vn3 = genvn(:(x[1]))
                vn4 = genvn(:(z[1][1]))
                vn5 = genvn(:(z[2]))
                vn6 = genvn(:z)

                @test in(vn1, space)
                @test in(vn2, space)
                @test in(vn3, space)
                @test in(vn4, space)
                @test ~in(vn5, space)
                @test ~in(vn6, space)
            end
            test_in()
        end
        vi = VarInfo()
        test_base!(vi)
        test_base!(empty!(TypedVarInfo(vi)))
    end
    @turing_testset "runmodel!" begin
        # Test that eval_num is incremented when calling runmodel!
        @model testmodel() = begin
            x ~ Normal()
        end
        alg = HMC(0.1, 5)
        spl = Sampler(alg, testmodel())
        vi = VarInfo()
        m = testmodel()
        m(vi)
        runmodel!(m, vi, spl)
        @test spl.state.eval_num == 1
        runmodel!(m, vi, spl)
        @test spl.state.eval_num == 2
    end
    @turing_testset "flags" begin
        # Test flag setting:
        #    is_flagged, set_flag!, unset_flag!
        function test_varinfo!(vi)
            vn_x = VarName(gensym(), :x, "", 1)
            dist = Normal(0, 1)
            r = rand(dist)
            gid = Selector()

            push!(vi, vn_x, r, dist, gid)

            # del is set by default
            @test is_flagged(vi, vn_x, "del") == false

            set_flag!(vi, vn_x, "del")
            @test is_flagged(vi, vn_x, "del") == true

            unset_flag!(vi, vn_x, "del")
            @test is_flagged(vi, vn_x, "del") == false
        end
        vi = VarInfo()
        test_varinfo!(vi)
        test_varinfo!(empty!(TypedVarInfo(vi)))
    end
    @turing_testset "link!" begin
        # Test linking spl and vi:
        #    link!, invlink!, istrans
        @model gdemo(x, y) = begin
            s ~ InverseGamma(2,3)
            m ~ Uniform(0, 2)
            x ~ Normal(m, sqrt(s))
            y ~ Normal(m, sqrt(s))
        end
        model = gdemo(1.0, 2.0)

        vi = VarInfo()
        model(vi, SampleFromUniform())

        @test all(x -> ~istrans(vi, x), vi.vns)
        alg = HMC(0.1, 5)
        spl = Sampler(alg, model)
        v = copy(vi.vals)
        link!(vi, spl)
        @test all(x -> istrans(vi, x), vi.vns)
        invlink!(vi, spl)
        @test all(x -> ~istrans(vi, x), vi.vns)
        @test vi.vals == v

        vi = TypedVarInfo(vi)
        alg = HMC(0.1, 5)
        spl = Sampler(alg, model)
        @test all(x -> ~istrans(vi, x), vi.metadata.s.vns)
        @test all(x -> ~istrans(vi, x), vi.metadata.m.vns)
        v_s = copy(vi.metadata.s.vals)
        v_m = copy(vi.metadata.m.vals)
        link!(vi, spl)
        @test all(x -> istrans(vi, x), vi.metadata.s.vns)
        @test all(x -> istrans(vi, x), vi.metadata.m.vns)
        invlink!(vi, spl)
        @test all(x -> ~istrans(vi, x), vi.metadata.s.vns)
        @test all(x -> ~istrans(vi, x), vi.metadata.m.vns)
        @test vi.metadata.s.vals == v_s
        @test vi.metadata.m.vals == v_m
    end
    @turing_testset "setgid!" begin
        vi = VarInfo()
        vn = VarName(gensym(), :x, "", 1)
        dist = Normal(0, 1)
        r = rand(dist)
        gid1 = Selector()
        gid2 = Selector(2, :HMC)

        push!(vi, vn, r, dist, gid1)
        @test vi.gids[vi.idcs[vn]] == Set([gid1])
        setgid!(vi, gid2, vn)
        @test vi.gids[vi.idcs[vn]] == Set([gid1, gid2])

        vi = empty!(TypedVarInfo(vi))
        push!(vi, vn, r, dist, gid1)
        @test vi.metadata.x.gids[vi.metadata.x.idcs[vn]] == Set([gid1])
        setgid!(vi, gid2, vn)
        @test vi.metadata.x.gids[vi.metadata.x.idcs[vn]] == Set([gid1, gid2])
    end
    @testset "orders" begin
        function randr(vi::VarInfo, vn::VarName, dist::Distribution, spl::Turing.Sampler)
            if ~haskey(vi, vn)
                r = rand(dist)
                Turing.push!(vi, vn, r, dist, spl)
                r
            elseif is_flagged(vi, vn, "del")
                unset_flag!(vi, vn, "del")
                r = rand(dist)
                vi[vn] = Turing.vectorize(dist, r)
                Turing.setorder!(vi, vn, vi.num_produce)
                r
            else
                Turing.updategid!(vi, vn, spl)
                vi[vn]
            end
        end

        csym = gensym() # unique per model
        vn_z1 = VarName(csym, :z, "[1]", 1)
        vn_z2 = VarName(csym, :z, "[2]", 1)
        vn_z3 = VarName(csym, :z, "[3]", 1)
        vn_z4 = VarName(csym, :z, "[4]", 1)
        vn_a1 = VarName(csym, :a, "[1]", 1)
        vn_a2 = VarName(csym, :a, "[2]", 1)
        vn_b = VarName(csym, :b, "", 1)

        vi = VarInfo()
        dists = [Categorical([0.7, 0.3]), Normal()]

        spl1 = Turing.Sampler(PG(5), empty_model())
        spl2 = Turing.Sampler(PG(5), empty_model())

        # First iteration, variables are added to vi
        # variables samples in order: z1,a1,z2,a2,z3
        vi.num_produce += 1
        randr(vi, vn_z1, dists[1], spl1)
        randr(vi, vn_a1, dists[2], spl1)
        vi.num_produce += 1
        randr(vi, vn_b, dists[2], spl2)
        randr(vi, vn_z2, dists[1], spl1)
        randr(vi, vn_a2, dists[2], spl1)
        vi.num_produce += 1
        randr(vi, vn_z3, dists[1], spl1)
        @test vi.orders == [1, 1, 2, 2, 2, 3]
        @test vi.num_produce == 3

        vi.num_produce = 0
        set_retained_vns_del_by_spl!(vi, spl1)
        @test is_flagged(vi, vn_z1, "del")
        @test is_flagged(vi, vn_a1, "del")
        @test is_flagged(vi, vn_z2, "del")
        @test is_flagged(vi, vn_a2, "del")
        @test is_flagged(vi, vn_z3, "del")

        vi.num_produce += 1
        randr(vi, vn_z1, dists[1], spl1)
        randr(vi, vn_a1, dists[2], spl1)
        vi.num_produce += 1
        randr(vi, vn_z2, dists[1], spl1)
        vi.num_produce += 1
        randr(vi, vn_z3, dists[1], spl1)
        randr(vi, vn_a2, dists[2], spl1)
        @test vi.orders == [1, 1, 2, 2, 3, 3]
        @test vi.num_produce == 3

        vi = empty!(TypedVarInfo(vi))
        # First iteration, variables are added to vi
        # variables samples in order: z1,a1,z2,a2,z3
        vi.num_produce += 1
        randr(vi, vn_z1, dists[1], spl1)
        randr(vi, vn_a1, dists[2], spl1)
        vi.num_produce += 1
        randr(vi, vn_b, dists[2], spl2)
        randr(vi, vn_z2, dists[1], spl1)
        randr(vi, vn_a2, dists[2], spl1)
        vi.num_produce += 1
        randr(vi, vn_z3, dists[1], spl1)
        @test vi.metadata.z.orders == [1, 2, 3]
        @test vi.metadata.a.orders == [1, 2]
        @test vi.metadata.b.orders == [2]
        @test vi.num_produce == 3

        vi.num_produce = 0
        set_retained_vns_del_by_spl!(vi, spl1)
        @test is_flagged(vi, vn_z1, "del")
        @test is_flagged(vi, vn_a1, "del")
        @test is_flagged(vi, vn_z2, "del")
        @test is_flagged(vi, vn_a2, "del")
        @test is_flagged(vi, vn_z3, "del")

        vi.num_produce += 1
        randr(vi, vn_z1, dists[1], spl1)
        randr(vi, vn_a1, dists[2], spl1)
        vi.num_produce += 1
        randr(vi, vn_z2, dists[1], spl1)
        vi.num_produce += 1
        randr(vi, vn_z3, dists[1], spl1)
        randr(vi, vn_a2, dists[2], spl1)
        @test vi.metadata.z.orders == [1, 2, 3]
        @test vi.metadata.a.orders == [1, 3]
        @test vi.metadata.b.orders == [2]
        @test vi.num_produce == 3
    end
    @turing_testset "replay" begin
        # Generate synthesised data
        xs = rand(Normal(0.5, 1), 100)

        # Define model
        @model priorsinarray(xs, ::Type{T}=Float64) where {T} = begin
            priors = Vector{T}(undef, 2)
            priors[1] ~ InverseGamma(2, 3)
            priors[2] ~ Normal(0, sqrt(priors[1]))
            for i = 1:length(xs)
                xs[i] ~ Normal(priors[2], sqrt(priors[1]))
            end
            priors
        end

        # Sampling
        chain = sample(priorsinarray(xs), HMC(0.01, 10), 10)
    end
    @turing_testset "varname" begin
        csym = gensym()
        vn1 = VarName(csym, :x, "[1]", 1)
        @test vn1 == VarName{:x}(csym, "[1]", 1)
        @test vn1 == VarName{:x}(csym, "[1]")
        @test vn1 == VarName([csym, :x], "[1]")
        vn2 = VarName(csym, :x, "[2]", 1)
        @test vn2 == VarName(vn1, "[2]")

        # Symbol
        v_sym = string(:x)
        @test v_sym == "x"

        # Array
        v_arr = eval(varname(:(x[i]))[1])
        @test v_arr == "[1]"

        # Matrix
        v_mat = eval(varname(:(x[i,j]))[1])
        @test v_mat== "[1,2]"

        v_mat = eval(varname(:(x[i,j,k]))[1])
        @test v_mat== "[1,2,3]"

        v_mat = eval(varname(:((x[1,2][1+5][45][3][i])))[1])
        @test v_mat == "[1,2][6][45][3][1]"

        @model mat_name_test() = begin
            p = Array{Any}(undef, 2, 2)
            for i in 1:2, j in 1:2
                p[i,j] ~ Normal(0, 1)
            end
            p
        end
        chain = sample(mat_name_test(), HMC(0.2, 4), 1000)
        check_numerical(chain, ["p[1, 1]"], [0], atol = 0.25)

        # Multi array
        v_arrarr = eval(varname(:(x[i][j]))[1])
        @test v_arrarr == "[1][2]"

        @model marr_name_test() = begin
            p = Array{Array{Any}}(undef, 2)
            p[1] = Array{Any}(undef, 2)
            p[2] = Array{Any}(undef, 2)
            for i in 1:2, j in 1:2
                p[i][j] ~ Normal(0, 1)
            end
            p
        end

        chain = sample(marr_name_test(), HMC(0.2, 4), 1000)
        check_numerical(chain, ["p[1][1]"], [0], atol = 0.25)
    end
    @turing_testset "varinfo" begin
        dists = [Normal(0, 1), MvNormal([0; 0], [1.0 0; 0 1.0]), Wishart(7, [1 0.5; 0.5 1])]
        function test_varinfo!(vi)
            @test getlogp(vi) == 0
            setlogp!(vi, 1)
            @test getlogp(vi) == 1
            acclogp!(vi, 1)
            @test getlogp(vi) == 2
            resetlogp!(vi)
            @test getlogp(vi) == 0

            spl2 = Turing.Sampler(PG(5, :w, :u), empty_model())
            vn_w = VarName(gensym(), :w, "", 1)
            randr(vi, vn_w, dists[1], spl2, true)

            vn_x = VarName(gensym(), :x, "", 1)
            vn_y = VarName(gensym(), :y, "", 1)
            vn_z = VarName(gensym(), :z, "", 1)
            vns = [vn_x, vn_y, vn_z]

            spl1 = Turing.Sampler(PG(5, :x, :y, :z), empty_model())
            for i = 1:3
                r = randr(vi, vns[i], dists[i], spl1, false)
                val = vi[vns[i]]
                @test sum(val - r) <= 1e-9
            end

            idcs = _getidcs(vi, spl1)
            if idcs isa NamedTuple
                @test sum(length(getfield(idcs, f)) for f in fieldnames(typeof(idcs))) == 3
            else
                @test length(idcs) == 3
            end
            @test length(vi[spl1]) == 7

            idcs = _getidcs(vi, spl2)
            if idcs isa NamedTuple
                @test sum(length(getfield(idcs, f)) for f in fieldnames(typeof(idcs))) == 1
            else
                @test length(idcs) == 1
            end
            @test length(vi[spl2]) == 1

            vn_u = VarName(gensym(), :u, "", 1)
            randr(vi, vn_u, dists[1], spl2, true)

            idcs = _getidcs(vi, spl2)
            if idcs isa NamedTuple
                @test sum(length(getfield(idcs, f)) for f in fieldnames(typeof(idcs))) == 2
            else
                @test length(idcs) == 2
            end
            @test length(vi[spl2]) == 2
        end
        vi = VarInfo()
        test_varinfo!(vi)
        test_varinfo!(empty!(TypedVarInfo(vi)))

        @model igtest() = begin
            x ~ InverseGamma(2,3)
            y ~ InverseGamma(2,3)
            z ~ InverseGamma(2,3)
            w ~ InverseGamma(2,3)
            u ~ InverseGamma(2,3)
        end

        # Test the update of group IDs
        g_demo_f = igtest()

        # This test section no longer seems as applicable, considering the
        # user will never end up using an UntypedVarInfo. The `VarInfo`
        # Varible is also not passed around in the same way as it used to be.
        g = Turing.Sampler(Gibbs(PG(10, :x, :y, :z), HMC(0.4, 8, :w, :u)), g_demo_f)
        pg, hmc = g.state.samplers
        vi = VarInfo()
        g_demo_f(vi, SampleFromPrior())
        Turing.Inference.step!(Random.GLOBAL_RNG, g_demo_f, pg, 1)
        vi1 = pg.state.vi
        @test mapreduce(x -> x.gids, vcat, vi1.metadata) ==
            [Set([pg.selector]), Set([pg.selector]), Set([pg.selector]), Set{Selector}(), Set{Selector}()]

        g_demo_f(vi1, hmc)
        @test mapreduce(x -> x.gids, vcat, vi1.metadata) ==
            [Set([pg.selector]), Set([pg.selector]), Set([pg.selector]), Set([hmc.selector]), Set([hmc.selector])]

        g = Turing.Sampler(Gibbs(PG(10, :x, :y, :z), HMC(0.4, 8, :w, :u)), g_demo_f)
        pg, hmc = g.state.samplers
        vi = empty!(TypedVarInfo(vi))
        g_demo_f(vi, SampleFromPrior())
        pg.state.vi = vi
        Turing.Inference.step!(Random.GLOBAL_RNG, g_demo_f, pg, 1)
        vi = pg.state.vi
        g_demo_f(vi, hmc)
        @test vi.metadata.x.gids[1] == Set([pg.selector])
        @test vi.metadata.y.gids[1] == Set([pg.selector])
        @test vi.metadata.z.gids[1] == Set([pg.selector])
        @test vi.metadata.w.gids[1] == Set([hmc.selector])
        @test vi.metadata.u.gids[1] == Set([hmc.selector])
    end
end


text:
@testset "Test utility functions" begin
    concrete_types = IS.get_all_concrete_subtypes(IS.InfrastructureSystemsComponent)
    @test length([x for x in concrete_types if isconcretetype(x)]) == length(concrete_types)
end

@testset "Test strip_module_name" begin
    @test IS.strip_module_name("PowerSystems.HydroDispatch") == "HydroDispatch"

    @test IS.strip_module_name(
        "InfrastructureSystems.SingleTimeSeries{PowerSystems.HydroDispatch}",
    ) == "SingleTimeSeries{PowerSystems.HydroDispatch}"

    @test IS.strip_module_name("SingleTimeSeries{PowerSystems.HydroDispatch}") ==
          "SingleTimeSeries{PowerSystems.HydroDispatch}"
end

@testset "Test strip_parametric_type" begin
    @test IS.strip_parametric_type("SingleTimeSeries{PowerSystems.HydroDispatch}") ==
          "SingleTimeSeries"

    @test IS.strip_parametric_type(
        "InfrastructureSystems.SingleTimeSeries{PowerSystems.HydroDispatch}",
    ) == "InfrastructureSystems.SingleTimeSeries"
end

@testset "Test exported names" begin
    @test IS.validate_exported_names(IS)
end

IS.@scoped_enum Fruit APPLE = 1 ORANGE = 2

@testset "Test scoped_enum" begin
    @test Fruit.APPLE isa Fruit
    @test Fruit.ORANGE isa Fruit
    @test sort([Fruit.ORANGE, Fruit.APPLE]) == [Fruit.APPLE, Fruit.ORANGE]

    struct Foo
        fruit::Fruit
    end
    @test Foo(1) == Foo(Fruit.APPLE)
end

@testset "Test undef component prints" begin
    v = Vector{IS.InfrastructureSystemsComponent}(undef, 3)
    @test sprint(show, v) ==
          "InfrastructureSystems.InfrastructureSystemsComponent[#undef, #undef, #undef]"
end

struct FakeTimeSeries <: InfrastructureSystems.TimeSeriesData end
Base.length(::FakeTimeSeries) = 42

@testset "Test TimeSeriesData printing" begin
    @test sprint(show, MIME("text/plain"), FakeTimeSeries()) ==
          "FakeTimeSeries time_series (42):"
end

ext: jl


text:
#=
To provide some helping functions which are more MATLAB like style
=#

function modified_reshape(matrix; row_number=0, column_number=0)
  #=
  This function allows you not to specify all the dimension of a reshape matrix
  Simply write the dimension that you want to change and let the remaining part
  is done by function itself
  =#
  row, column = size(matrix)
  number_of_entries = row * column

  if row_number == 0 && column_number == 0
    return matrix

  elseif row_number != 0 && column_number == 0
    if number_of_entries % row_number != 0
      throw(DimensionMismatch("Dimensions must be consistent with input matrix"))
    end
    return reshape(matrix, row_number, convert(Int64, number_of_entries / row_number))

  elseif row_number == 0 && column_number != 0
    if number_of_entries % column_number != 0
      throw(DimensionMismatch("Dimensions must be consistent with input matrix"))
    end
    return reshape(matrix, convert(Int64, number_of_entries / column_number), column_number)

  else
    if row_number * column_number != number_of_entries
      throw(DimensionMismatch("Dimensions must be consistent with input matrix"))
    end
      return reshape(matrix, row_number, column_number)
  end
end

function network_initializer(init_wt, numhid1, numhid2, numwords, vocab_size)
  #=
  This function is designed to create initial weights and delta matrices
  for backpropagation algorithm

  =#
  word_embedding_weights = init_wt * randn(vocab_size, numhid1)
  embed_to_hid_weights = init_wt * randn(numwords * numhid1, numhid2)
  hid_to_output_weigths = init_wt * randn(numhid2, vocab_size)
  hid_bias = zeros(numhid2, 1)
  output_bias = zeros(vocab_size, 1)

  word_embedding_weights_delta = zeros(vocab_size, numhid1)
  word_embedding_weights_gradient = zeros(vocab_size, numhid1)

  embed_to_hid_weights_delta = zeros(numwords * numhid1, numhid2)

  hid_to_output_weights_delta = zeros(numhid2, vocab_size)

  hid_bias_delta = zeros(numhid2, 1)

  output_bias_delta = zeros(vocab_size, 1)

return (word_embedding_weights,
        embed_to_hid_weights,
        hid_to_output_weigths,
        hid_bias,
        output_bias,
        word_embedding_weights_delta,
        word_embedding_weights_gradient,
        embed_to_hid_weights_delta,
        hid_to_output_weights_delta,
        hid_bias_delta,
        output_bias_delta)

end

function findmaxcols(matrix)
  #=
  Input: matrix by size of M X N
  returns : a row vector containing maximum of each column in rows.

  =#
  result = zeros(1, size(matrix,2))
  for i=1:size(matrix, 2)
    result[1,i] = maximum(matrix[:,i])
  end
  return result
end

function sum_columns(matrix)
  #=
  MATLAB like summer to sum up the columns of a matrix

  =#
  result = zeros(1, size(matrix,2))
  for i=1:size(matrix, 2)
    result[1, i] = sum(matrix[:, i])
  end
  return result
end

function sum_rows(matrix)
  #=
  MATLAB like summer to sum up the rows of a matrix
  returns a column vector in where the each column is corresponding sum.
  =#
  result = zeros(size(matrix, 1), 1)
  for i=1:size(matrix, 1)
    result[i, 1] = sum(matrix[i, :])
  end
  return result
end

function output_mapper(expansion_matrix, target_batch)
#=
Since Julia suggests devoctarizes implementations,
this function creates a mapping matrix which maps the target values a vector of
1 and 0s
=#
  batchsize = size(target_batch)[2]
  result_matrix = zeros(size(expansion_matrix)[1], batchsize)
  for i=1:batchsize
    result_matrix[:, i] = expansion_matrix[:, target_batch[i]]
  end
  return result_matrix
end

max_stars_repo_name: kirnap/Machine-Learning-in-Julia

max_issues_repo_name: kirnap/Machine-Learning-in-Julia

max_forks_repo_name: kirnap/Machine-Learning-in-Julia


text:
module JLD2Tests

using Test
using Gridap.Helpers
import Gridap.Io: to_jld2_file
import Gridap.Io: from_jld2_file

struct Foo <: GridapType
  bar::Int
end

bar = 3
foo = Foo(bar)

d = mktempdir()
f = joinpath(d,"foo.jld2")

to_jld2_file(foo,f)
@test foo == from_jld2_file(typeof(foo),f)

f = joinpath(d,"dict.jld2")
foo = Dict("a"=>Int32(1),2=>Int(3),4.0=>Float32(5),"six"=>Float64(7),:s=>"Symbol")

to_jld2_file(foo,f)
@test foo == from_jld2_file(typeof(foo),f)

rm(d,recursive=true)

end # module


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    PointwiseLearn(model)

A learning solver that converts spatial data to a tabular format
with features (and possibly labels) for each point, and then solves
the problem with statistical learning `model`.

## Parameters

* `model` - Learning model (e.g. SVM, Logistic Regression, K-means)

## Notes

Any model implementing the `MLJModelInterface.jl` can be used in
pointwise learning. Please refer to the `MLJ` documentation for
a list of available models.
"""
struct PointwiseLearn{M} <: AbstractLearningSolver
  model::M
end

function solve(problem::LearningProblem, solver::PointwiseLearn)
  ptask = task(problem)
  model = solver.model

  # assert model is compatible with task
  @assert iscompatible(model, ptask) "$model is not compatible with $ptask"

  # learn model on source data
  lmodel = learn(ptask, sourcedata(problem), model)

  # apply model to target data
  result = perform(ptask, targetdata(problem), lmodel)

  LearningSolution(domain(targetdata(problem)), result)
end

# ------------
# IO methods
# ------------
function Base.show(io::IO, solver::PointwiseLearn)
  print(io, "PointwiseLearn")
end

function Base.show(io::IO, ::MIME"text/plain", solver::PointwiseLearn)
  println(io, solver)
  print(io, "  └─model ⇨ ")
  show(IOContext(io, :compact => true), solver.model)
  println(io, "")
end


text:
using Base.Test
using RiemannTheta
reload("RiemannTheta")

g, num_vectors = 5, 10
z = [ rand(ComplexF64, g) for i in 1:num_vectors]
mode, ϵ = 0, 1e-6
accuracy_radius = 5.
tttt = rand(5,5)
Y = tttt * tttt'
Ω = Complex.(rand(g, g), Y)

derivs = Vector{ComplexF64}[]
oscillatory_part(z, Ω, 0, ϵ, derivs, accuracy_radius)
@btime oscillatory_part($z, $Ω, 0, $ϵ, $derivs, $accuracy_radius) # 380-390ms / 430 Mb

derivs = [ rand(ComplexF64, 5) for i in 1:5 ]
@btime oscillatory_part($z, $Ω, 0, $ϵ, $derivs, $accuracy_radius) # 6.9s / 5.6Gb

@testset "oscillatory_part derivs are correct" begin
    z₀ = rand(ComplexF64, g) - Complex(0.5, 0.5)
    δ = [Complex(1e-8,1e-8); zeros(ComplexF64, g-1)]
    derivs = Vector{ComplexF64}[]
    ϵ = 1e-8
    # calculate function at slightly shifted z₀
    z = [[z₀] ; [z₀ + circshift(δ, i) for i in 0:g-1]]
    res = oscillatory_part(z, Ω, 0, ϵ, derivs, accuracy_radius)
    dres = [ res[i] - res[1] for i in 2:g+1 ]
    dres

    # and compare to calculated derivates
    for i in 0:g-1
        derivs = [ circshift(δ, i) ]
        res2 = oscillatory_part([z₀], Ω, 0, ϵ, derivs, accuracy_radius)
        @test isapprox(res2[1], dres[i+1], rtol=1e-4)
    end
end


text:
__precompile__()
module InvSmooth2d_mod
using LinearAlgebra
export Smooth2d!

include("lsq_cgsolver.jl")

using  .LSQ_CGSolver

mutable struct CGsmooth2d
#
# CG Data Structure:
n1::Int64
n2::Int64
nd::Int64
dat::Array{Float64,1}
scl::Float64
# cg precond
nlag::Int64
pcd::Float64
#
# Constructor:
function CGsmooth2d(img::Array{Float64,2})
n1::Int64=size(img,1)
n2::Int64=size(img,2)
nd::Int64 = n1*n2
dat::Array{Float64,1}=zeros(Float64,nd)
#
id::Int64=0
for i2=1:n2
    for i1=1:n1
        id+=1
        dat[id]=img[i1,i2]
    end
end
#
scl::Float64=sqrt(Float64(nd))
pnlag::Int64=5
pcd::Float64=0.0
#
@inbounds for i2::Int64 in -pnlag:pnlag
@inbounds for i1::Int64 in -pnlag:pnlag
    pcd+=1.0
end
end
pcd=1.0/pcd
#
new(n1,n2,nd,dat,scl,pnlag,pcd)
#
end
end
#
#traits
has_dataweight(cg::CGsmooth2d)=true
has_operator(cg::CGsmooth2d)=true
has_preconditioning(cg::CGsmooth2d)=true
has_modelweight(cg::CGsmooth2d)=true
has_cg_operators(cg::CGsmooth2d)=true
#
# Identity operator
function oper!(adj::Int64,
               add::Int64,
               nm::Int64,
               mdl::Array{Float64,1},
               nd::Int64,
               dat::Array{Float64,1},
               cg::CGsmooth2d)
iexit::Int64=0
nd::Int64
if adj == 0
   if add == 0
      @inbounds dat[1:nd].=0.0
   end
@simd   for id in 1:nd
       @inbounds dat[id] = mdl[id]
   end
else
   if add == 0
      @inbounds mdl[1:nd].=0.0
   end
@simd   for id  in 1:nd
       @inbounds mdl[id] = dat[id]
   end
end
return iexit
end
#
# Box preconditioning
function pcnd!(adj::Int64,
               add::Int64,
               nm::Int64,
               mdl::Array{Float64,1},
               nd::Int64,
               dat::Array{Float64,1},
               cg::CGsmooth2d)

ic::Int64=0
iexit::Int64=0
if adj == 0
   #
   if add == 0
      @inbounds dat[1:nd] .= 0.0
   end
   #
   #
   for i2=cg.nlag+1:cg.n2-cg.nlag
   for i1=cg.nlag+1:cg.n1-cg.nlag
       #
       ic=i1+(i2-1)*cg.n1
       #
       for l2::Int64 in -cg.nlag:cg.nlag
           j2::Int64=i2+l2
           for l1::Int64 in -cg.nlag:cg.nlag
               j1::Int64=i1+l1
               if 0 < j2 && j2 <= cg.n2
               if 0 < j1 && j1 <= cg.n1
               idx::Int64=j1+(j2-1)*cg.n1
               @inbounds dat[ic] += cg.pcd*mdl[idx]
               end
               end
           end
       end
   end
   end
#
else
#
   if add == 0
      @inbounds mdl[1:nm].=0.0
   end
   #
   for i2=cg.nlag+1:cg.n2-cg.nlag
   for i1=cg.nlag+1:cg.n1-cg.nlag
       #
       ic=i1+(i2-1)*cg.n1
       #
       for l2::Int64 in -cg.nlag:cg.nlag
           j2::Int64=i2+l2
           for l1::Int64 in -cg.nlag:cg.nlag
               j1::Int64=i1+l1
               if 0 < j2 && j2 <= cg.n2
               if 0 < j1 && j1 <= cg.n1
                  idx::Int64=j1+(j2-1)*cg.n1
                  @inbounds mdl[idx] += cg.pcd*dat[ic]
               end
               end
           end
       end
       #
   end
   end
   #
end
return iexit
end

# data weight
wdat!(adj::Int64,
           nd::Int64,
           dat::Array{Float64,1},
           cg::CGsmooth2d)=0

# model selection weight
wmdl!(adj::Int64,
      nc::Int64,
      mdl::Array{Float64,1},
      cg::CGsmooth2d)=0

function Smooth2d!(img::Array{Float64,2},
                   imgs::Array{Float64,2},
                   lmbd::Float64)
#
n1::Int64=size(img,1)
n2::Int64=size(img,2)
#
@assert size(imgs,1)==n1 && size(imgs,2)==n2
#
cg::CGsmooth2d=CGsmooth2d(img)
mdl::Array{Float64,1}=zeros(Float64,n1*n2)
niter::Int64=n1*n2
eps::Float64=0.000000000001
#
iexit::Int64=lsq_cgd(cg,
                     InvSmooth2d_mod.wdat!,
                     InvSmooth2d_mod.oper!,
                     InvSmooth2d_mod.wmdl!,
                     InvSmooth2d_mod.pcnd!,
                     cg.nd,mdl,
                     cg.nd,cg.dat,
                     niter,
                     lmbd*cg.scl,
                     eps)
#

nb::Int64=3*cg.nlag+1
id::Int64=0
for i2=nb+1:n2-nb
    for i1=nb+1:n1-nb
        id=i1+(i2-1)*n1
        imgs[i1,i2]=mdl[id]
    end
end
#=
 Boundary conditions:
=#

for i2=1:nb
   for i1=nb+1:n1-nb
       imgs[i1,i2]=imgs[i1,nb+1]
   end
end
for i2=n2-nb+1:n2
   for i1=nb+1:n1-nb
       imgs[i1,i2]=imgs[i1,n2-nb]
   end
end

for i2=1:n2
    for i1=1:nb
        imgs[i1,i2]      =imgs[nb+1,i2]
        imgs[n1-nb+i1,i2]=imgs[n1-nb,i2]
    end
end

scl::Float64=0.0
for i2=1:n2
    for i1=1:n1
        scl+= imgs[i1,i2]^2
    end
end

if scl > 0.0
   scl=sqrt(dot(cg.dat,cg.dat))/sqrt(scl)
end

for i2=1:n2
    for i1=1:n1
        imgs[i1,i2]=scl*imgs[i1,i2]
    end
end
return nothing

end

end #module


text:
struct LazDataset
    filename::String
    filehandle::Ptr{Cvoid}
    header::LazHeader  # this enables iterating without unsafe_load everytime
    point::Ptr{LazPoint}
end

function Base.show(io::IO, ds::LazDataset)
    n = max(ds.header.number_of_point_records, ds.header.extended_number_of_point_records)
    println(io, "LazDataset of $(ds.filename) with $n points.")
end

function bounds(ds::LazDataset)
    bounds(ds.header)
end

function open(f::AbstractString)
    # Setup laszip reader
    laszip_reader = Ref{Ptr{Cvoid}}()
    @check laszip_reader[] laszip_create(laszip_reader)

    # Open lasfile
    @check laszip_reader[] laszip_open_reader(laszip_reader[], f, Ref{Cint}(0))

    # Get a pointer to the header
    header_ptr = Ref{Ptr{LazHeader}}()
    @check laszip_reader[] laszip_get_header_pointer(laszip_reader[], header_ptr)
    header = unsafe_load(header_ptr[])

    # Get a pointer to the points that will be read
    point_ptr = Ref{Ptr{LazPoint}}()
    @check laszip_reader[] laszip_get_point_pointer(laszip_reader[], point_ptr)

    LazDataset(f, laszip_reader[], header, point_ptr[])
end

"""Iteration of LAZ file."""
function Base.iterate(ds::LazDataset, state::Int)
    if state >= length(ds)
        return nothing
    else
        laszip_read_point(ds.filehandle)
        return unsafe_load(ds.point), state + 1
    end
end

function Base.iterate(ds::LazDataset)
    # Be sure to seek to 0 at begin of iteration
    if length(ds) == 0
        nothing
    else
        laszip_seek_point(ds.filehandle, 0)
        laszip_read_point(ds.filehandle)
        unsafe_load(ds.point), 1
    end
end

Base.eltype(::LazDataset) = LazPoint
Base.length(ds::LazDataset) = Int(max(ds.header.number_of_point_records, ds.header.extended_number_of_point_records))

function Base.close(ds::LazDataset)
    reader = ds.filehandle
    LazIO.@check reader LazIO.laszip_close_reader(reader)
    LazIO.@check reader LazIO.laszip_destroy(reader)
end

function LasIO.boundingbox(h::LazHeader)
    (xmin = h.min_x, ymin = h.min_y, zmin = h.min_z, xmax = h.max_x, ymax = h.max_y, zmax = h.max_z)
end
LasIO.boundingbox(ds::LazDataset) = boundingbox(ds.header)

max_stars_repo_name: evetion/LazIO

max_issues_repo_name: evetion/LazIO

max_forks_repo_name: evetion/LazIO


text:
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{VAF},
    ::Type{DisjunctiveConstraints.DisjunctiveSet{Float64}},
)
    return true
end

function MOI.is_valid(
    opt::Optimizer,
    c::CI{VAF,DisjunctiveConstraints.DisjunctiveSet{Float64}},
)
    return 1 <= c.value <= num_disjunctive_constraints(opt.form)
end

function MOI.add_constraint(
    opt::Optimizer,
    _f::VAF,
    s::DisjunctiveConstraints.DisjunctiveSet{Float64},
)
    f = convert.(CSAF, MOIU.scalarize(_f))
    disj = Disjunction(f, s)
    # TODO: Make this configurable
    df = DisjunctiveFormulater(
        disj,
        DisjunctiveConstraints.NaiveBigM(opt.config.activity_method),
    )
    attach_formulater!(opt.form, df)
    return CI{VAF,DisjunctiveConstraints.DisjunctiveSet{Float64}}(
        num_disjunctive_constraints(opt.form),
    )
end

function MOI.get(
    opt::Optimizer,
    ::MOI.NumberOfConstraints{
        VAF,
        DisjunctiveConstraints.DisjunctiveSet{Float64},
    },
)
    return num_disjunctive_constraints(opt.form)
end

function MOI.get(
    opt::Optimizer,
    ::MOI.ListOfConstraintIndices{
        VAF,
        DisjunctiveConstraints.DisjunctiveSet{Float64},
    },
)
    return [
        CI{VAF,DisjunctiveConstraints.DisjunctiveSet{Float64}}(i) for
        i in 1:num_disjunctive_constraints(opt.form)
    ]
end


text:
include(Pkg.dir("FinancialAnalytics/src/helpers.jl"))
using MarketData

r = percentchange(cl)

facts("keyword checking is correct") do

  context("simple returns selected") do
      @fact sum(keyword_check(r,false,false).values) => -0.8392600357038683
  end

  context("lsum(og returns selected") do
      @fact sum(keyword_check(percentchange(cl,method="log"),false,true).values) => -0.8392600357038683
  end

  context("ssum(imple prices selected") do
      @fact sum(keyword_check(cl,true,false).values) => -0.8392600357038683
  end

  context("lsum(og prices selected") do
      @fact sum(keyword_check(basecall(cl,log),true,true).values) => -0.8392600357038683
  end
end


text:
import PredictMD
import Test

Test.@test_throws ErrorException PredictMD.trapz([1, 2, 3], [])
Test.@test_throws ErrorException PredictMD.trapz([], [])
Test.@test_throws ErrorException PredictMD.trapz([1, 2, 1], [4, 5, 4])

x_1 = collect(0:0.00001:1)
f_1(t) = t^2
y_1 = f_1.(x_1)
I_1 = PredictMD.trapz(x_1, y_1)
Test.@test isapprox(I_1, 1/3; atol=0.00000001)

x_2 = collect(0:0.00001:1)
f_2(t) = (t^2) * (tanh(t))
y_2 = f_2.(x_2)
I_2 = PredictMD.trapz(x_2, y_2)
Test.@test isapprox(I_2, 0.207068855098706896; atol=0.00000001)


text:
include("grammar.jl")

import Base

"""
This type represents an [`Automaton`](@ref) transition. It has a `from` starting
*state* and a `to` destination *state* and a `label`

    Transition(from::Union{AbstractString, Iterable}, label::AbstractString, to::Union{AbstractString, Iterable})
Build a transition based on the given states.
"""
struct Transition
    "The starting state(s) of the transition"
    from::Union{AbstractString, Set}
    "The label of the transition"
    label::AbstractString
    "The destination starte(s) of the transition"
    to::Union{AbstractString, Set}
    function Transition(from::Union{AbstractString, Iterable}, label::AbstractString, to::Union{AbstractString, Iterable})
        check(s::AbstractString)::Bool = ~isempty(s)
        check(s::AbstractProduction)::Bool = true
        check(s)::Bool = false
        check(s::Iterable) = all(x -> check(x), Set(collect(s))) && ~isempty(s)

        f = nothing
        l = nothing
        t = nothing

        if check(from)
            f = isa(from, Iterable) ? Set(collect(from)) : from
        else
            throw(ArgumentError("The from state is not a nonempty string, or a nonempty set of nonempty strings/items"))
        end

        if check(label)
            l = label
        else
            throw(ArgumentError("The label is not a nonempty string"))
        end

        if check(to)
            t = isa(to, Iterable) ? Set(collect(to)) : to
        else
            throw(ArgumentError("The to state is not a nonempty string, or a nonempty set of nonempty strings/items"))
        end

        return new(f, l, t)
    end
end

### Functions ###

"""
    parsetransitions(t::AbstractString)::Array{Transition}
Builds an array of [`Transition`](@ref) obtained from the given string.
"""
function parsetransitions(t::AbstractString)::Array{Transition}
    res = Transition[]
    for t ∈ split(t, '\n')
        if isempty(strip(t))
            continue
        end
        from, label, to = split(t, ',')
        push!(res, Transition(strip(from), strip(label), strip(to)))
    end
    return res
end

### Operators ###

Base.:<(t1::Transition, t2::Transition) = (t1.from, t1.label, t1.to) < (t2.from, t2.label, t2.to)
Base.:(==)(t1::Transition, t2::Transition) = (t1.from, t1.label, t1.to) == (t2.from, t2.label, t2.to)
Base.:>(t1::Transition, t2::Transition) = !(t1 == t2 || t1 < t2)
Base.hash(t::Transition) = Base.hash((t.from, t.label, t.to))
Base.show(io::IO, e::Transition) = Base.show(io, string(e.from, "-", e.label, "->", e.to))
Base.iterate(x::Transition, i...) = Base.iterate((x.from, x.label, x.to))


text:
# Autogenerated wrapper script for LibYAML_jll for x86_64-unknown-freebsd
export libyaml

JLLWrappers.@generate_wrapper_header("LibYAML")
JLLWrappers.@declare_library_product(libyaml, "libyaml-0.so.2")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libyaml,
        "lib/libyaml-0.so.2",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
baremodule HistoricAnalysis

function simplified end
function unstack_tasks end

module Internal

using DataFrames
using Tables

import Historic
using Historic.Internal: EventRecord

using ..HistoricAnalysis: HistoricAnalysis

include("utils.jl")
include("simplified.jl")
include("unstack_tasks.jl")

end  # module Internal

end  # baremodule HistoricAnalysis


text:
##

# Convert Julia types to LLVM types
const LLVMTYPE = Dict{DataType,String}(
    Bool => "i8",   # Julia represents Tuple{Bool} as [1 x i8]
    Int8 => "i8",
    Int16 => "i16",
    Int32 => "i32",
    Int64 => "i64",
    Int128 => "i128",
    UInt8 => "i8",
    UInt16 => "i16",
    UInt32 => "i32",
    UInt64 => "i64",
    UInt128 => "i128",
    Float16 => "half",
    Float32 => "float",
    Float64 => "double",
    Nothing => "void"
)
llvmtype(x)::String = LLVMTYPE[x]
const JuliaPointerType = LLVMTYPE[Int]

const NativeTypes = Union{Bool, Base.HWReal}
# The `@eval` loop is to avoid type ambiguities on Julia 1.0. Later versions did not have this problem

# function vload_quote(::Type{T})
# end

const SCOPE_METADATA = """
!1 = !{!\"noaliasdomain\"}
!2 = !{!\"noaliasscope\", !1}
!3 = !{!2}
"""
const LOAD_SCOPE_TBAA = SCOPE_METADATA * """
!4 = !{!"jtbaa"}
!5 = !{!6, !6, i64 0, i64 0}
!6 = !{!"jtbaa_arraybuf", !4, i64 0}
"""
const STORE_TBAA = """
!4 = !{!"jtbaa", !5, i64 0}
!5 = !{!"jtbaa"}
!6 = !{!"jtbaa_data", !4, i64 0}
!7 = !{!8, !8, i64 0}
!8 = !{!"jtbaa_arraybuf", !6, i64 0}
"""

function vload_quote(::Type{T}) where {T <: NativeTypes}
    ptyp = JuliaPointerType
    typ = llvmtype(T)
    instrs = String[]
    decl = LOAD_SCOPE_TBAA
    alignment = Base.datatype_alignment(T)
    # push!(flags, "!noalias !0")
    push!(instrs, "%ptr = inttoptr $ptyp %0 to $typ*")
    push!(instrs, "%res = load $typ, $typ* %ptr, align $alignment, !alias.scope !3, !tbaa !5")
    push!(instrs, "ret $typ %res")
    :(llvmcall($((decl,join(instrs, "\n"))), $T, Tuple{Ptr{$T}}, ptr))
end
function vload_quote(::Type{T}, ::Type{I}) where {T <: NativeTypes, I <: Integer}
    ityp = 'i' * string(8sizeof(I))
    ptyp = JuliaPointerType
    typ = llvmtype(T)
    instrs = String[]
    decl = LOAD_SCOPE_TBAA
    alignment = Base.datatype_alignment(T)
    # push!(flags, "!noalias !0")
    push!(instrs, "%typptr = inttoptr $ptyp %0 to i8*")
    push!(instrs, "%iptr = getelementptr inbounds i8, i8* %typptr, $ityp %1")
    push!(instrs, "%ptr = bitcast i8* %iptr to $typ*")
    push!(instrs, "%res = load $typ, $typ* %ptr, align $alignment, !alias.scope !3, !tbaa !5")
    push!(instrs, "ret $typ %res")
    :(llvmcall($((decl, join(instrs, "\n"))), $T, Tuple{Ptr{$T}, $I}, ptr, i))
end
function vstore_quote(::Type{T}, alias) where {T <: NativeTypes}
    ptyp = JuliaPointerType
    typ = llvmtype(T)
    instrs = String[]
    alignment = Base.datatype_alignment(T)
    # push!(flags, "!noalias !0")
    decl = alias ? STORE_TBAA : SCOPE_METADATA * STORE_TBAA
    push!(instrs, "%ptr = inttoptr $ptyp %0 to $typ*")
    aliasstoreinstr = "store $typ %1, $typ* %ptr, align $alignment, !tbaa !7"
    noaliasstoreinstr = "store $typ %1, $typ* %ptr, align $alignment, !noalias !3, !tbaa !7"
    push!(instrs, alias ? aliasstoreinstr : noaliasstoreinstr)
    push!(instrs, "ret void")
    :(llvmcall($((decl, join(instrs, "\n"))), Cvoid, Tuple{Ptr{$T}, $T}, ptr, v))
end
function vstore_quote(::Type{T}, ::Type{I}, alias) where {T <: NativeTypes, I <: Integer}
    ityp = 'i' * string(8sizeof(I))
    ptyp = JuliaPointerType
    typ = llvmtype(T)
    instrs = String[]
    decl = alias ? STORE_TBAA : SCOPE_METADATA * STORE_TBAA
    alignment = Base.datatype_alignment(T)
    # push!(flags, "!noalias !0")
    push!(instrs, "%typptr = inttoptr $ptyp %0 to i8*")
    push!(instrs, "%iptr = getelementptr inbounds i8, i8* %typptr, $ityp %2")
    push!(instrs, "%ptr = bitcast i8* %iptr to $typ*")
    aliasstoreinstr = "store $typ %1, $typ* %ptr, align $alignment, !tbaa !7"
    noaliasstoreinstr = "store $typ %1, $typ* %ptr, align $alignment, !noalias !3, !tbaa !7"
    push!(instrs, alias ? aliasstoreinstr : noaliasstoreinstr)
    push!(instrs, "ret void")
    :(llvmcall($((decl,join(instrs, "\n"))), Cvoid, Tuple{Ptr{$T}, $T, $I}, ptr, v, i))
end

function gepquote(::Type{T}, ::Type{I}, byte::Bool) where {T <: NativeTypes, I <: Integer}
    ptyp = JuliaPointerType
    ityp = llvmtype(I)
    typ = byte ? "i8" : llvmtype(T)
    instrs = String[]
    push!(instrs, "%ptr = inttoptr $ptyp %0 to $typ*")
    push!(instrs, "%offsetptr = getelementptr inbounds $typ, $typ* %ptr, $ityp %1")
    push!(instrs, "%iptr = ptrtoint $typ* %offsetptr to $ptyp")
    push!(instrs, "ret $ptyp %iptr")
    quote
        llvmcall(
            $(join(instrs, "\n")),
            Ptr{$T}, Tuple{Ptr{$T}, $I},
            ptr, i
        )
    end
end

@inline vload(x::Number, args...) = x
for T ∈ [Bool,Int8,Int16,Int32,Int64,UInt8,UInt16,UInt32,UInt64,Float32,Float64]
    @eval @inline vload(ptr::Ptr{$T}) = $(vload_quote(T))
    @eval @inline vstore!(ptr::Ptr{$T}, v::$T) = $(vstore_quote(T, true))
    @eval @inline vnoaliasstore!(ptr::Ptr{$T}, v::$T) = $(vstore_quote(T, false))
    for I ∈ [Int32,UInt32,Int64,UInt64]
        @eval @inline vload(ptr::Ptr{$T}, i::$I) = $(vload_quote(T, I))
        @eval @inline vstore!(ptr::Ptr{$T}, v::$T, i::$I) = $(vstore_quote(T, I, true))
        @eval @inline vnoaliasstore!(ptr::Ptr{$T}, v::$T, i::$I) = $(vstore_quote(T, I, false))
        # @eval @inline gep(ptr::Ptr{$T}, i::$I) = $(gepquote(T, I, false))
        @eval @inline Base.@pure gep(ptr::Ptr{$T}, i::$I) = $(gepquote(T, I, false))
        @eval @inline Base.@pure gepbyte(ptr::Ptr{$T}, i::$I) = $(gepquote(T, I, true))
        # @eval @inline gep(ptr::Ptr{$T}, i::$I) = $(gepquote(T, I, false))
        # @eval @inline gepbyte(ptr::Ptr{$T}, i::$I) = $(gepquote(T, I, true))
    end
end
@inline vload(ptr, ::Static{N}) where {N} = vload(ptr, N)
@inline vstore!(ptr, v, ::Static{N}) where {N} = vstore!(ptr, v, N)

# Fall back definitions
@inline vload(ptr::Ptr) = Base.unsafe_load(ptr)
# @inline vload(ptr::Ptr, i::Int) = vload(gep(ptr, i))
@inline vstore!(ptr::Ptr{T}, v::T) where {T} = Base.unsafe_store!(ptr, v)
@inline vload(::Type{T1}, ptr::Ptr{T2}) where {T1, T2} = vload(Base.unsafe_convert(Ptr{T1}, ptr))
@inline vstore!(ptr::Ptr{T1}, v::T2) where {T1,T2} = vstore!(ptr, convert(T1, v))
@inline vstore!(ptr::Ptr{T1}, v::T2) where {T1<:Integer,T2<:Integer} = vstore!(ptr, v % T1)
# @inline vstore!(ptr::Ptr{T1}, v::T2) where {T1,T2} = vstore!(ptr, convert(T1, v))

@inline gep(ptr::Ptr{Mask{8,UInt8}}, i) = Base.unsafe_convert(Ptr{Mask{8,UInt8}}, gep(Base.unsafe_convert(Ptr{UInt8}, ptr), i >> 3))
@inline gepbyte(ptr::Ptr{Mask{8,UInt8}}, i::Integer) = Base.unsafe_convert(Ptr{Mask{8,UInt8}}, gep(Base.unsafe_convert(Ptr{UInt8}, ptr), i))
@inline vload(ptr::Ptr{Mask{8,UInt8}}, i) = (vload(Base.unsafe_convert(Ptr{UInt8}, ptr), i >> 3) >> (i & 7)) % Bool
@inline vload(ptr::Ptr{Mask{8,UInt8}}, i::_MM{8}) = Mask{8}(vload(Base.unsafe_convert(Ptr{UInt8}, ptr), i.i >> 3))
@inline vload(ptr::Ptr{Mask{8,UInt8}}, i::_MM{16}) = Mask{16}(vload(Base.unsafe_convert(Ptr{UInt16}, gep(ptr, i.i))))
@inline vload(ptr::Ptr{Mask{8,UInt8}}, i::_MM{32}) = Mask{32}(vload(Base.unsafe_convert(Ptr{UInt32}, gep(ptr, i.i))))
@inline vload(ptr::Ptr{Mask{8,UInt8}}, i::_MM{64}) = Mask{64}(vload(Base.unsafe_convert(Ptr{UInt64}, gep(ptr, i.i))))
@inline function vstore!(ptr::Ptr{Mask{8,UInt8}}, v::Bool, i)
    p = gep(Base.unsafe_convert(Ptr{UInt8}, ptr), i >> 3)
    u = vload(p)
    m = 0x01 << (i & 7)
    u = ifelse(v, u | m, u ⊻ m)
    vstore!(p, u); v
end
@inline vstore!(ptr::Ptr{Mask{8,UInt8}}, v::UInt8, i::_MM{8}) = Mask{8}(vstore!(Base.unsafe_convert(Ptr{UInt8}, ptr), i.i >> 3))
@inline vstore!(ptr::Ptr{Mask{8,UInt8}}, v::UInt16, i::_MM{16}) = Mask{16}(vstore!(Base.unsafe_convert(Ptr{UInt16}, gep(ptr, i.i >> 3))))
@inline vstore!(ptr::Ptr{Mask{8,UInt8}}, v::UInt32, i::_MM{32}) = Mask{32}(vstore!(Base.unsafe_convert(Ptr{UInt32}, gep(ptr, i.i >> 3))))
@inline vstore!(ptr::Ptr{Mask{8,UInt8}}, v::UInt64, i::_MM{64}) = Mask{64}(vstore!(Base.unsafe_convert(Ptr{UInt64}, gep(ptr, i.i >> 3))))

@inline tdot(::Tuple{}, ::Tuple{}) = Zero()
@inline tdot(a::Tuple{I1}, b::Tuple{I2}) where {I1,I2} = @inbounds vmul(a[1], b[1])
@inline tdot(a::Tuple{I1,I3}, b::Tuple{I2,I4}) where {I1,I2,I3,I4} = @inbounds vadd(vmul(a[1],b[1]), vmul(a[2],b[2]))
@inline tdot(a::Tuple{I1,I3,I5}, b::Tuple{I2,I4,I6}) where {I1,I2,I3,I4,I5,I6} = @inbounds vadd(vadd(vmul(a[1],b[1]), vmul(a[2],b[2])), vmul(a[3],b[3]))
# @inline tdot(a::NTuple{N,Int}, b::NTuple{N,Int}) where {N} = @inbounds a[1]*b[1] + tdot(Base.tail(a), Base.tail(b))

@inline tdot(a::Tuple{I}, ::Tuple{}) where {I} = Zero() #@inbounds a[1]
@inline tdot(::Tuple{}, b::Tuple{I}) where {I} = Zero() #@inbounds b[1]
@inline tdot(a::Tuple{I1,Vararg}, b::Tuple{I2}) where {I1,I2} = @inbounds vmul(a[1],b[1])
@inline tdot(a::Tuple{I1}, b::Tuple{I2,Vararg}) where {I1,I2} = @inbounds vmul(a[1],b[1])
@inline tdot(a::Tuple{I1,Vararg}, b::Tuple{I2,Vararg}) where {I1,I2} = @inbounds vadd(vmul(a[1],b[1]), tdot(Base.tail(a), Base.tail(b)))

"""
A wrapper to the base pointer type, that supports pointer arithmetic.
Note that `VectorizationBase.vload` and `VectorizationBase.vstore!` are 0-indexed,
while `Base.unsafe_load` and `Base.unsafe_store!` are 1-indexed.
x = [1, 2, 3, 4, 5, 6, 7, 8];
ptrx = Pointer(x);
vload(ptrx)
# 1
vload(ptrx + 1)
# 2
ptrx[]
# 1
(ptrx+1)[]
# 2
ptrx[1]
# 1
ptrx[2]
# 2
"""
abstract type AbstractPointer{T} end

@generated function gepbyte(ptr::Ptr{T}, i::I) where {T, I <: Integer}
    q = gepquote(T, I, true)
    pushfirst!(q.args, Expr(:meta, :inline))
    q
end
@generated function gep(ptr::Ptr{T}, i::I) where {T <: NativeTypes, I <: Integer}
    q = gepquote(T, I, false)
    pushfirst!(q.args, Expr(:meta, :inline))
    q
end

@generated function gep(ptr::Ptr{T}, i::_Vec{_W,I}) where {_W, T <: NativeTypes, I <: Integer}
    W = _W + 1
    ptyp = JuliaPointerType
    typ = llvmtype(T)
    ityp = llvmtype(I)
    vityp = "<$W x $ityp>"
    vptyp = "<$W x $ptyp>"
    instrs = String[]
    push!(instrs, "%ptr = inttoptr $ptyp %0 to $typ*")
    push!(instrs, "%offsetptr = getelementptr inbounds $typ, $typ* %ptr, $vityp %1")
    push!(instrs, "%iptr = ptrtoint <$W x $typ*> %offsetptr to $vptyp")
    push!(instrs, "ret $vptyp %iptr")
    quote
        $(Expr(:meta, :inline))
        llvmcall(
            $(join(instrs, "\n")),
            NTuple{$W,Core.VecElement{Ptr{$T}}}, Tuple{Ptr{$T}, NTuple{$W,Core.VecElement{$I}}},
            ptr, i
        )
    end
end
@inline gep(ptr::Ptr{T}, i::_Vec{_W,I}, ::Val{0}) where {T <: NativeTypes, _W, I <: Integer} = ptr
@generated function gep(ptr::Ptr{T}, i::_Vec{_W,I}, ::Val{N}) where {T <: NativeTypes, _W, I <: Integer, N}
    if N ∉ [1,2,4,8]
        tz = trailing_zeros(N)
        if iszero(tz) || isone(count_ones(N))
            return :(Base.@_inline_meta; gepbyte(ptr, vmul(i,N)))
        else
            return :(Base.@_inline_meta; gep(ptr, vmul(i, $(N >> tz)), Val{$(1 << tz)}()))
        end
    end
    W = _W + 1
    ptyp = "i$(8 * sizeof(Int))"
    # We can remove a bitcast from the optimized IR
    styp = N == sizeof(T) ? llvmtype(T) : "i$(8 * N)"
    ityp = "i$(8 * sizeof(I))"
    instrs = String[]
    push!(instrs, "%ptr = inttoptr $ptyp %0 to $styp*")
    push!(instrs, "%offsetptr = getelementptr inbounds $styp, $styp* %ptr, <$W x $ityp> %1")
    push!(instrs, "%iptr = ptrtoint <$W x $styp*> %offsetptr to <$W x $ptyp>")
    push!(instrs, "ret <$W x $ptyp> %iptr")
    quote
        Base.@_inline_meta
        llvmcall(
            $(join(instrs, "\n")),
            NTuple{$W,Core.VecElement{Ptr{$T}}}, Tuple{Ptr{$T}, NTuple{$W,Core.VecElement{$I}}},
            ptr, i
        )
    end
end
@inline gep(ptr::Ptr{T}, i::I, ::Val{0}) where {T <: NativeTypes, I <:Integer} = ptr
@generated function gep(ptr::Ptr{T}, i::I, ::Val{N}) where {T <: NativeTypes, I <: Integer, N}
    if N ∉ [1,2,4,8]
        tz = trailing_zeros(N)
        if iszero(tz) || isone(count_ones(N))
            return :(Base.@_inline_meta; gepbyte(ptr, vmul(i,N)))
        else
            return :(Base.@_inline_meta; gep(ptr, vmul(i, $(N >> tz)), Val{$(1 << tz)}()))
        end
    end
    ptyp = "i$(8 * sizeof(Int))"
    # We can remove a bitcast from the optimized IR
    styp = N == sizeof(T) ? llvmtype(T) : "i$(8 * N)"
    ityp = "i$(8 * sizeof(I))"
    instrs = String[]
    push!(instrs, "%ptr = inttoptr $ptyp %0 to $styp*")
    push!(instrs, "%offsetptr = getelementptr inbounds $styp, $styp* %ptr, $ityp %1")
    push!(instrs, "%iptr = ptrtoint $styp* %offsetptr to $ptyp")
    push!(instrs, "ret $ptyp %iptr")
    quote
        Base.@_inline_meta
        llvmcall(
            $(join(instrs, "\n")),
            Ptr{$T}, Tuple{Ptr{$T}, $I},
            ptr, i
        )
    end
end

@inline gep(ptr::Ptr, ::Static{0}) = ptr
@inline gepbyte(ptr::Ptr, ::Static{0}) = ptr
# @inline gepbyte(ptr::AbstractPointer, ::Static{0}) = pointer(ptr)
@inline gep(ptr::Ptr, ::Static{N}) where {N} = gep(ptr, N)
@inline gepbyte(ptr::Ptr, ::Static{N}) where {N} = gepbyte(ptr, N)
@inline gepbyte(ptr::Ptr, i::LazyStaticMul{N}) where {N} = gep(ptr, i.data, Val{N}())
@inline gep(ptr::Ptr, v::SVec) = gep(ptr, extract_data(v))
@inline gep(ptr::Ptr{Cvoid}, i::Integer) = gepbyte(ptr, i)

# struct Reference{T} <: AbstractPointer{T}
    # ptr::Ptr{T}
    # @inline Reference(ptr::Ptr{T}) where {T} = new{T}(ptr)
# end
@inline Base.eltype(::AbstractPointer{T}) where {T} = T
# @inline gep(ptr::Pointer, i::Tuple{<:Integer}) = gep(ptr, first(i))
# @inline vstore!(ptr::AbstractPointer{T1}, v::T2, args...) where {T1,T2} = vstore!(ptr, convert(T1, v), args...)

@inline vload(::Type{SVec{W,T}}, ptr::AbstractPointer) where {W,T} = vload(SVec{W,T}, pointer(ptr))
@inline vloadstride(ptr::AbstractPointer{T}, str1, strd::Union{_MM,AbstractSIMDVector}) where {T} = vload(pointer(ptr), vadd(vmul(static_sizeof(T), str1), strd))
@inline vloadstride(ptr::AbstractPointer{T}, str1, strd::Union{_MM,AbstractSIMDVector}, u::Union{AbstractMask,Unsigned}) where {T} = vload(pointer(ptr), vadd(vmul(static_sizeof(T), str1), strd), u)

@inline vloadstride(ptr::AbstractPointer{T}, str1, strd) where {T} = vload(gepbyte(pointer(ptr), strd), vmulnp(static_sizeof(T), str1))
@inline vloadstride(ptr::AbstractPointer{T}, str1, strd, u::Union{AbstractMask,Unsigned}) where {T} = vload(gepbyte(pointer(ptr), strd), vmulnp(static_sizeof(T), str1), u)

# @inline vloadstride(ptr::AbstractPointer{T}, str1::Integer, strd) where {T} = vload(gep(gepbyte(pointer(ptr), strd), str1))
# @inline vloadstride(ptr::AbstractPointer{T}, str1::Integer, strd, u::Union{AbstractMask,Unsigned}) where {T} = vload(gep(gepbyte(pointer(ptr), strd), str1), u)
# @inline vloadstride(ptr::AbstractPointer{T}, str1::SVec{<:Any,<:Integer}, strd) where {T} = vload(gep(gepbyte(pointer(ptr), strd), str1))
# @inline vloadstride(ptr::AbstractPointer{T}, str1::SVec{<:Any,<:Integer}, strd, u::Union{AbstractMask,Unsigned}) where {T} = vload(gep(gepbyte(pointer(ptr), strd), str1), u)

@inline vload(ptr::AbstractPointer{T}, i::Tuple) where {T} = vloadstride(ptr, stride1offset(ptr, i), stridedoffset(ptr, i))
@inline vload(ptr::AbstractPointer{T}, i::Tuple, u::Union{AbstractMask,Unsigned}) where {T} = vloadstride(ptr, stride1offset(ptr, i), stridedoffset(ptr, i), u)

for (fs, f) ∈ [(:vstorestride!, :vstore!), (:vnoaliasstorestride!, :vnoaliasstore!)]
    @eval begin
        @inline function $fs(ptr::AbstractPointer{T}, v, str1, strd) where {T}
            $f(gepbyte(pointer(ptr), strd), v, vmulnp(static_sizeof(T), str1))
        end
        @inline function $fs(ptr::AbstractPointer{T}, v, str1, strd::Union{_MM,AbstractSIMDVector}) where {T}
            $f(pointer(ptr), v, vadd(vmul(static_sizeof(T), str1), strd))
        end
        @inline function $fs(ptr::AbstractPointer{T}, v, str1, strd, u::Union{AbstractMask,Unsigned}) where {T}
            $f(gepbyte(pointer(ptr), strd), v, vmulnp(static_sizeof(T), str1), u)
        end
        @inline function $fs(ptr::AbstractPointer{T}, v, str1, strd::Union{_MM,AbstractSIMDVector}, u::Union{AbstractMask,Unsigned}) where {T}
            $f(pointer(ptr), v, vadd(vmul(static_sizeof(T), str1), strd), u)
        end
    end
end

@inline vstore!(ptr::AbstractPointer, v, i::Tuple) = vstorestride!(ptr, v, stride1offset(ptr, i), stridedoffset(ptr, i))
@inline vstore!(ptr::AbstractPointer, v, i::Tuple, u::Union{AbstractMask,Unsigned}) = vstorestride!(ptr, v, stride1offset(ptr, i), stridedoffset(ptr, i), u)
@inline vnoaliasstore!(ptr::AbstractPointer, v, i::Tuple) = vnoaliasstorestride!(ptr, v, stride1offset(ptr, i), stridedoffset(ptr, i))
@inline vnoaliasstore!(ptr::AbstractPointer, v, i::Tuple, u::Union{AbstractMask,Unsigned}) = vnoaliasstorestride!(ptr, v, stride1offset(ptr, i), stridedoffset(ptr, i), u)

@inline vnoaliasstore!(args...) = vstore!(args...) # generic fallback

@inline vstore!(ptr::Ptr{T}, v::Number, i::Integer) where {T <: Number} = vstore!(ptr, convert(T, v), i)
@inline vstore!(ptr::Ptr{T}, v::Integer, i::Integer) where {T <: Integer} = vstore!(ptr, v % T, i)
# @inline vstore!(ptr::AbstractPointer{T}, v::T, i...) where {T<:Integer} = vstore!(pointer(ptr), v, offset(ptr, i...))
# @inline vstore!(ptr::AbstractPointer{T1}, v::T2, i...) where {T1<:Integer,T2<:Integer} = vstore!(pointer(ptr), v % T1, offset(ptr, i...))
# @inline vstore!(ptr::AbstractPointer{T1}, v::T2, i...) where {T1<:Number,T2<:Number} = vstore!(pointer(ptr), convert(T1,v), offset(ptr, i...))
# @inline vstore!(ptr::AbstractPointer{T}, v::T, i...) where {T<:Number} = vstore!(pointer(ptr), v, offset(ptr, i...))
# @inline vstore!(ptr::AbstractPointer{T}, v::T, i...) where {T<:Integer} = vstore!(pointer(ptr), v, offset(ptr, i...))
# @inline vstore!(ptr::AbstractPointer{T1}, v::T2, i...) where {T1<:Integer,T2<:Integer} = vstore!(pointer(ptr), v % T1, offset(ptr, i...))
# @inline vstore!(ptr::AbstractPointer{T1}, v::T2, i...) where {T1<:Number,T2<:Number} = vstore!(pointer(ptr), convert(T1,v), offset(ptr, i...))

abstract type AbstractStridedPointer{T} <: AbstractPointer{T} end
abstract type AbstractColumnMajorStridedPointer{T,N} <: AbstractStridedPointer{T} end
struct PackedStridedPointer{T,N} <: AbstractColumnMajorStridedPointer{T,N}
    ptr::Ptr{T}
    strides::NTuple{N,Int}
end
struct PackedStridedBitPointer{Nm1,N} <: AbstractColumnMajorStridedPointer{Bool,Nm1}
    ptr::Ptr{UInt64}
    strides::NTuple{Nm1,Int}
    offsets::NTuple{N,Int}
end
@inline function gesp(ptr::PackedStridedBitPointer, i::Tuple)
    PackedStridedBitPointer(ptr.ptr, ptr.strides, vadd(i, ptr.offsets))
end
struct PaddedStridedBitPointer{N} <: AbstractColumnMajorStridedPointer{Bool,N}
    ptr::Ptr{Mask{8,UInt8}}#use Mask for dispatch
    strides::NTuple{N,Int}
end
PaddedStridedBitPointer(p::Ptr{<:Unsigned}, s) = PaddedStridedBitPointer(Base.unsafe_convert(Ptr{Mask{8,UInt8}}, p), s)
PaddedStridedBitPointer(A::StridedArray{T}) where {T <: Unsigned} = PaddedStridedBitPointer(pointer(A), staticmul(T,tailstrides(A)))

abstract type AbstractRowMajorStridedPointer{T,N} <: AbstractStridedPointer{T} end
struct RowMajorStridedPointer{T,N} <: AbstractRowMajorStridedPointer{T,N}
    ptr::Ptr{T}
    strides::NTuple{N,Int}
end
struct RowMajorStridedBitPointer{N} <: AbstractRowMajorStridedPointer{Bool,N}
    ptr::Ptr{UInt64}
    strides::NTuple{N,Int}
end

struct PermutedDimsStridedPointer{S1,S2,T,P<:AbstractStridedPointer{T}} <: AbstractStridedPointer{T}
    ptr::P
end
@inline PermutedDimsStridedPointer{S1,S2}(ptr::P) where {S1,S2, T, P <: AbstractStridedPointer{T}} = PermutedDimsStridedPointer{S1,S2,T,P}(ptr)
@inline function resort_tuple(i::Tuple{Vararg{<:Any,N}}, ::Val{S}) where {N,S}
    ntuple(Val{N}()) do n
        i[S[n]]
    end
end
@inline Base.pointer(ptr::PermutedDimsStridedPointer) = pointer(ptr.ptr)
@inline function stridedpointer(A::PermutedDimsArray{T,N,S1,S2}) where {T,N,S1,S2}
    PermutedDimsStridedPointer{S1,S2}(stridedpointer(parent(A)))
end
@inline pointerforcomparison(ptr::AbstractStridedPointer, i::Tuple) = gep(ptr, i)
@inline pointerforcomparison(ptr::AbstractStridedPointer) = pointer(ptr)

@generated function subsetview(ptr::PermutedDimsStridedPointer{S1,S2}, ::Val{I}, i) where {S1,S2,I}
    s1new = Expr(:tuple)
    s2new = Expr(:tuple)
    resize!(s2new.args, length(S1)-1)
    Iinner = S1[I]
    i = 0
    for s ∈ S1
        if s < Iinner
            i += 1
            push!(s1new.args, s)
            s2new.args[s] = i
        elseif s > Iinner
            i += 1
            push!(s1new.args, s - 1)
            s2new.args[s-1] = i
        end
    end
    sv = :(subsetview(ptr.ptr, Val{$Iinner}(), i))
    quote
        $(Expr(:meta,:inline))
        PermutedDimsStridedPointer{$s1new,$s2new}($sv)
    end
end

# @inline function stridedpointer(A::PermutedDimsArray{T,2,(2,1),(2,1)}) where {T}
#     RowMajorStridedPointer(stridedp
# end

abstract type AbstractSparseStridedPointer{T,N} <: AbstractStridedPointer{T} end
struct SparseStridedPointer{T,N} <: AbstractSparseStridedPointer{T,N}
    ptr::Ptr{T}
    strides::NTuple{N,Int}
end
# struct SparseStridedBitPointer{N} <: AbstractSparseStridedPointer{Bool,N}
#     ptr::Ptr{UInt64}
#     strides::NTuple{N,Int}
# end

abstract type AbstractStaticStridedPointer{T,X} <: AbstractStridedPointer{T} end
struct StaticStridedPointer{T,X} <: AbstractStaticStridedPointer{T,X}
    ptr::Ptr{T}
end
struct StaticStridedBitPointer{X} <: AbstractStaticStridedPointer{Bool,X}
    ptr::Ptr{UInt64}
end
const AbstractBitPointer = Union{PackedStridedBitPointer, RowMajorStridedBitPointer, StaticStridedBitPointer}#, SparseStridedBitPointer
@inline pointerforcomparison(ptr::AbstractBitPointer, i::Tuple) = gesp(ptr, i)
@inline pointerforcomparison(ptr::AbstractBitPointer) = ptr

@inline offset(ptr::AbstractStridedPointer{T}, i::Tuple) where {T} = vmuladdnp(static_sizeof(T), stride1offset(ptr, i), stridedoffset(ptr, i))

@inline stride1offset(::AbstractColumnMajorStridedPointer, i::Tuple{}) = Zero()
@inline stride1offset(::AbstractColumnMajorStridedPointer, i::Tuple) = first(i)
@inline stridedoffset(ptr::AbstractColumnMajorStridedPointer, i::Tuple{I}) where {I} = Zero()
@inline stridedoffset(ptr::AbstractColumnMajorStridedPointer, i::Tuple{I1,I2,Vararg}) where {I1,I2} = tdot(Base.tail(i), ptr.strides)

# @inline offset(ptr::AbstractSparseStridedPointer, i::Integer) = i * @inbounds ptr.strides[1]
# @inline offset(ptr::AbstractStaticStridedPointer{<:Any,<:Tuple{1,Vararg}}, i::Integer) = i
# @inline offset(ptr::AbstractStaticStridedPointer{<:Any,<:Tuple{M,Vararg}}, i::Integer) where {M} = M*i
@inline gep(ptr::AbstractStridedPointer, i::Tuple) = gep(gepbyte(pointer(ptr), stridedoffset(ptr, i)), extract_data(stride1offset(ptr, i)))
# @inline gep(ptr::AbstractStridedPointer, i::Tuple{I}) where {I} = gepbyte(ptr.ptr, first(offset(ptr, i)))
@inline gepbyte(ptr::AbstractStridedPointer, i::Tuple) = gep(gepbyte(pointer(ptr), stridedoffset(ptr, i)), extract_data(stride1offset(ptr, i)))
# @inline gepbyte(ptr::AbstractStridedPointer, i::Tuple{I}) where {I} = gepbyte(ptr.ptr, first(offset(ptr, i)))
@inline gesp(ptr::AbstractStridedPointer, i) = similar(ptr, gep(ptr, i))

@inline Base.similar(p::PackedStridedPointer, ptr::Ptr) = PackedStridedPointer(ptr, p.strides)
@inline Base.similar(p::PaddedStridedBitPointer, ptr::Ptr) = PaddedStridedBitPointer(ptr, p.strides)
@inline Base.similar(p::PackedStridedBitPointer{Nm1,N}, ptr::Ptr) where {Nm1,N} = PackedStridedBitPointer(ptr, p.strides, ntuple(_ -> 0, Val{N}()))
@inline Base.similar(p::RowMajorStridedPointer, ptr::Ptr) = RowMajorStridedPointer(ptr, p.strides)
@inline Base.similar(p::RowMajorStridedBitPointer, ptr::Ptr) = RowMajorStridedBitPointer(ptr, p.strides)
@inline Base.similar(p::SparseStridedPointer, ptr::Ptr) = SparseStridedPointer(ptr, p.strides)
# @inline Base.similar(p::SparseStridedBitPointer, ptr::Ptr) = SparseStridedBitPointer(ptr, p.strides)
@inline Base.similar(::StaticStridedPointer{T,X}, ptr::Ptr) where {T,X} = StaticStridedPointer{T,X}(ptr)
@inline Base.similar(::StaticStridedBitPointer{X}, ptr::Ptr) where {X} = StaticStridedBitPointer{X}(ptr)
@inline Base.similar(p::PermutedDimsStridedPointer{S1,S2}, ptr) where {S1,S2} = PermutedDimsStridedPointer{S1,S2}(similar(p.ptr, ptr))

@inline LinearAlgebra.transpose(ptr::RowMajorStridedPointer) = PackedStridedPointer(ptr.ptr, ptr.strides)
@inline stride1offset(::AbstractRowMajorStridedPointer{T,0}, i::Tuple{I}) where {T,I} = first(i)
@inline stride1offset(::AbstractRowMajorStridedPointer{T,1}, i::Tuple{I1,I2}) where {T,I1,I2} = last(i)
@inline stride1offset(::AbstractRowMajorStridedPointer{T,1}, i::Tuple{I1}) where {T,I1} = Zero()
@inline stride1offset(::AbstractRowMajorStridedPointer{T,N}, i::Tuple{I1,Vararg{<:Any,N}}) where {T,I1,N} = last(i)
# @inline stride1offset(::AbstractRowMajorStridedPointer{T,0}, i::Tuple) where {T} = last(i)

@inline stridedoffset(ptr::AbstractRowMajorStridedPointer{T,0}, i::Tuple) where {T} = Zero()
@inline stridedoffset(ptr::AbstractRowMajorStridedPointer{T,1}, i::Tuple{I1,I2}) where {T,I1,I2} = @inbounds vmul(i[1],ptr.strides[1])
@inline stridedoffset(ptr::AbstractRowMajorStridedPointer{T,2}, i::Tuple{I1,I2,I3}) where {T,I1,I2,I3} = @inbounds vadd(vmul(i[1],ptr.strides[2]), vmul(i[2],ptr.strides[1]))
@inline stridedoffset(ptr::AbstractRowMajorStridedPointer, i::Tuple) = tdot(reverse(ptr.strides), i)

@inline stride1offset(::AbstractSparseStridedPointer, ::Tuple) = Zero()
@inline stridedoffset(ptr::AbstractSparseStridedPointer, i::Tuple) = tdot(i, ptr.strides)

@generated function LinearAlgebra.transpose(ptr::StaticStridedPointer{T,X}) where {T,X}
    tup = Expr(:curly, :Tuple)
    N = length(X.parameters)
    for n ∈ N:-1:1
        push!(tup.args, (X.parameters[n])::Int)
    end
    Expr(:block, Expr(:meta, :inline), Expr(:call, Expr(:curly, :StaticStridedPointer, T, tup), Expr(:(.), :ptr, QuoteNode(:ptr))))
end

# @inline offset(ptr::AbstractStaticStridedPointer{T,<:Tuple{1,Vararg}}, i::Integer) where {T} = vmulnp(static_sizeof(T), i)
# @inline offset(ptr::AbstractStaticStridedPointer{T,<:Tuple{N,Vararg}}, i::Integer) where {N,T} = vmul(N * static_sizeof(T), i)
@inline offset(ptr::AbstractStaticStridedPointer{T,<:Tuple{1,Vararg}}, i::Tuple{I}) where {T,I<:Integer} = vmulnp(static_sizeof(T), first(i))
@inline offset(ptr::AbstractStaticStridedPointer{T,<:Tuple{N,Vararg}}, i::Tuple{I}) where {N,T,I<:Integer} = vmul(vmul(Static(N), static_sizeof(T)), first(i))

# @inline offset(ptr::AbstractStaticStridedPointer{T,Tuple{1,1}}, i::Tuple{<:_MM,<:Any}) where {T} = @inbounds vmuladdnp(static_sizeof(T), i[1], vmul(static_sizeof(T), i[2]))
# @inline offset(ptr::AbstractStaticStridedPointer{T,Tuple{1,1}}, i::Tuple{<:Any,<:_MM}) where {T} = @inbounds vmuladdnp(static_sizeof(T), i[2], vmul(static_sizeof(T), i[1]))
# @inline offset(ptr::AbstractStaticStridedPointer{T,Tuple{1,1}}, i::Tuple{<:Any,<:Any}) where {T} = @inbounds vadd(vmul(static_sizeof(T), i[1]), vmul(static_sizeof(T), i[2]))

function indprod(X, i, st)
    Xᵢ = (X[i])::Int
    iᵢ = Expr(:ref, :i, i)
    Expr(:call, :vmul, st * Xᵢ, iᵢ)
end
@generated function stride1offset(ptr::AbstractStaticStridedPointer{T,X}, i::I) where {T,X,I<:Tuple}
    N = length(I.parameters)
    Xv = X.parameters
    M = min(N, length(X.parameters))
    Xvv = Int[(Xv[m])::Int for m in 1:M]
    mstart = findfirst(isequal(1), Xvv)
    isnothing(mstart) && return Expr(:block, Expr(:meta, :inline), Zero())
    ret = Expr(:ref, :i, mstart)
    for m ∈ 1+mstart:M
        if isone(Xvv[m])
            ret = Expr(:call, :vadd, ret, Expr(:ref, :i, m))
        end
    end
    Expr(
        :block,
        Expr(:meta,:inline),
        Expr(
            :macrocall,
            Symbol("@inbounds"),
            LineNumberNode(@__LINE__, Symbol(@__FILE__)), ret
         )
    )
end
@generated function stridedoffset(ptr::AbstractStaticStridedPointer{T,X}, i::I) where {T,X,I<:Tuple}
    N = length(I.parameters)
    Xv = X.parameters
    M = min(N, length(X.parameters))
    Xvv = Int[(Xv[m])::Int for m in 1:M]
    mstart = findfirst(!isequal(1), Xvv)
    isnothing(mstart) && return Expr(:block, Expr(:meta, :inline), Zero())
    st = sizeof(T)
    ind = indprod(Xvv, mstart, st)
    for m ∈ 1+mstart:M
        if !isone(Xvv[m])
            ind = Expr(:call, :vadd, indprod(Xv, m, st), ind)
        end
    end
    Expr(
        :block,
        Expr(:meta,:inline),
        Expr(
            :macrocall,
            Symbol("@inbounds"),
            LineNumberNode(@__LINE__, Symbol(@__FILE__)), ind
         )
    )
end

struct StaticStridedStruct{T,X,S} <: AbstractStaticStridedPointer{T,X}
    ptr::S
    offset::Int # keeps track of offset, incase of nested offset calls
end

const AbstractInitializedStridedPointer{T} = Union{
    PackedStridedPointer{T},
    RowMajorStridedPointer{T},
    SparseStridedPointer{T},
    StaticStridedPointer{T},
    StaticStridedStruct{T}
}
const AbstractInitializedPointer{T} = Union{
    # Pointer{T},
    PackedStridedPointer{T},
    RowMajorStridedPointer{T},
    SparseStridedPointer{T},
    StaticStridedPointer{T},
    StaticStridedStruct{T}
}

# @inline Base.stride(ptr::AbstractColumnMajorStridedPointer, i) = isone(i) ? 1 : @inbounds ptr.strides[i-1]
# @inline Base.stride(ptr::AbstractSparseStridedPointer, i) = @inbounds ptr.strides[i]
# @generated function Base.stride(::AbstractStaticStridedPointer{T,X}, i) where {T,X}
#     Expr(:block, Expr(:meta, :inline), Expr(:getindex, Expr(:tuple, X.parameters...), :i))
# end
# @inline LinearAlgebra.stride1(ptr::AbstractColumnMajorStridedPointer) = 1
# @inline LinearAlgebra.stride1(ptr::AbstractSparseStridedPointer) = @inbounds first(ptr.strides)
# @inline LinearAlgebra.stride1(::AbstractStaticStridedPointer{T,<:Tuple{X,Vararg}}) where {T,X} = X

@inline offset(ptr::AbstractPointer, i::CartesianIndex) = offset(ptr, i.I)

@inline Base.:+(ptr::AbstractPointer{T}, i) where {T} = similar(ptr, gep(pointer(ptr), i))
@inline Base.:+(i, ptr::AbstractPointer{T}) where {T} = similar(ptr, gep(pointer(ptr), i))
@inline Base.:-(ptr::AbstractPointer{T}, i) where {T} = similar(ptr, gep(pointer(ptr), - i))

@inline vload(ptr::AbstractPointer) = vload(pointer(ptr))
@inline Base.unsafe_load(p::AbstractPointer) = vload(pointer(p))
@inline Base.getindex(p::AbstractPointer) = vload(pointer(p))

@inline vload(ptr::AbstractPointer{T}, ::Tuple{}) where {T} = vload(pointer(ptr))
# @inline vload(ptr::AbstractPointer, i::Tuple) = vload(pointer(ptr), offset(ptr, i))
# @inline vload(ptr::AbstractPointer, i::Tuple, u::Unsigned) = vload(pointer(ptr), offset(ptr, i), u)
@inline Base.unsafe_load(ptr::AbstractPointer, i) = vload(pointer(ptr), offset(ptr, vsub(i, 1)))
@inline Base.getindex(ptr::AbstractPointer, i) = vload(ptr, (i, ))
@inline Base.getindex(ptr::AbstractPointer, i, j) = vload(ptr, (i, j))
@inline Base.getindex(ptr::AbstractPointer, i, j, k) = vload(ptr, (i, j, k))
@inline Base.getindex(ptr::AbstractPointer, i, j, k, rests...) = vload(ptr, (i, j, k, rests...))

@inline vstore!(ptr::AbstractPointer{T}, v::T) where {T} = vstore!(pointer(ptr), v)
@inline Base.unsafe_store!(ptr::AbstractPointer{T}, v::T) where {T} = vstore!(pointer(ptr), v)
@inline Base.setindex!(ptr::AbstractPointer{T}, v::T) where {T} = vstore!(pointer(ptr), v)

# @inline vstore!(ptr::AbstractPointer{T}, v, i::Tuple) where {T} = vstore!(pointer(ptr), v, offset(ptr, i))
# @inline vstore!(ptr::AbstractPointer{T}, v, i::Tuple, u::Unsigned) where {T} = vstore!(pointer(ptr), v, offset(ptr, i), u)
@inline Base.unsafe_store!(ptr::AbstractPointer{T}, v::T, i) where {T} = vstore!(pointer(ptr), v, offset(ptr, vsub(i, 1)))
@inline Base.setindex!(ptr::AbstractPointer{T}, v::T, i) where {T} = vstore!(pointer(ptr), v, offset(ptr, i))

# @inline Pointer(A) = Pointer(pointer(A))
@inline Base.pointer(ptr::AbstractPointer) = ptr.ptr
@inline Base.unsafe_convert(::Type{Ptr{T}}, ptr::AbstractPointer{T}) where {T} = pointer(ptr)

@inline stridedpointer(x) = x#Pointer(x)

"""
    stridedpointer(x, i)

Gets the stridedpointer pointing to i, using 1-based indexes.
"""
@inline function stridedpointer(x, i)
    gesp(stridedpointer(x), (staticm1(i),))
end
@inline function stridedpointer(x, i1, i2, I...)
    gesp(stridedpointer(x), staticm1((i1, i2, I...)))
end
@inline stridedpointer(x::Ptr) = PackedStridedPointer(x, tuple())
# @inline stridedpointer(x::Union{LowerTriangular,UpperTriangular}) = stridedpointer(parent(x))
# @inline stridedpointer(x::AbstractArray) = stridedpointer(parent(x))
@inline tailstrides(A::AbstractArray) = Base.tail(strides(A))
@inline tailstrides(A::AbstractVector) = tuple()
@inline tailstrides(A::AbstractMatrix) = (stride(A,2),)
@inline tailstrides(A::BitArray{1}) = tuple()
@inline tailstrides(A::BitArray{2}) = (size(A,1),)
@inline tailstrides(A::BitArray{3}) = (size(A,1),size(A,1)*size(A,2))
@generated function tailstrides(A::BitArray{N}) where {N}
    quote
        (Base.Cartesian.@ntuple $(N-1) s) = size(A)
        Base.Cartesian.@nexprs $(N-2) n -> s_{n+1} = vmul(s_n, s_{n+1})
        (Base.Cartesian.@ntuple $(N-1) s)
    end
end
@inline stridedpointer(A::AbstractArray{T}) where {T} = PackedStridedPointer(pointer(A), staticmul(T,tailstrides(A)))
@inline stridedpointer(A::BitArray{N}) where {N} = PackedStridedBitPointer(pointer(A.chunks), tailstrides(A), ntuple(_ -> 0, Val{N}()))
@inline stridedpointer(A::AbstractArray{T,0}) where {T} = pointer(A)
@inline stridedpointer(A::SubArray{T,0,P,S}) where {T,P,S <: Tuple{Int,Vararg}} = pointer(A)
@inline stridedpointer(A::SubArray{T,N,P,S}) where {T,N,P,S <: Tuple{<:StepRange,Vararg}} = SparseStridedPointer(pointer(A), staticmul(T, strides(A)))
@inline stridedpointer(A::SubArray{T,N,P,S}) where {T,N,P,S <: Tuple{Int,Vararg}} = SparseStridedPointer(pointer(A), staticmul(T, strides(A)))
@inline stridedpointer(A::SubArray{T,N,P,S}) where {T,N,P,S} = PackedStridedPointer(pointer(A), staticmul(T, tailstrides(A)))
# Slow fallback
@inline stridedpointer(A::SubArray{T,N,P,S}) where {T,N,P<:PermutedDimsArray,S} = SparseStridedPointer(pointer(A), staticmul(T, strides(A)))
@inline function stridedpointer(A::SubArray{T,N,P,S}) where {S1,S2,T,N,P<:PermutedDimsArray{<:StridedArray,N,S1,S2},S<:Tuple{Vararg{AbstractUnitRange}}}
    ptr = similar(stridedpointer(parent(parent(A))), pointer(A))
    PermutedDimsStridedPointer{S1,S2}(ptr)
end
@inline stridedpointer(B::Union{Adjoint{T,A},Transpose{T,A}}) where {T, A <: DenseVector{T}} = StaticStridedPointer{T,Tuple{1,1}}(pointer(parent(B)))
@inline stridedpointer(B::Union{Adjoint{T,A},Transpose{T,A}}) where {T, A <: StridedVector{T}} = adjoint_vector_strided_pointer(stridedpointer(parent(B)))
@inline function adjoint_vector_strided_pointer(ptr::SparseStridedPointer)
    x = ptr.strides[1]
    SparseStridedPointer(pointer(ptr), (x,x))
end
@inline adjoint_vector_strided_pointer(ptr::PackedStridedPointer{T}) where {T} = StaticStridedPointer{T,Tuple{1,1}}(pointer(ptr))
@inline adjoint_vector_strided_pointer(ptr::StaticStridedPointer{T,Tuple{X}}) where {T,X} = StaticStridedPointer{T,Tuple{X,X}}(pointer(ptr))

@inline function stridedpointer(B::Union{Adjoint{T,A},Transpose{T,A}}) where {T,N,A <: AbstractArray{T,N}}
    pB = parent(B)
    RowMajorStridedPointer(pointer(pB), staticmul(T, tailstrides(pB)))
end
@inline function stridedpointer(B::Union{Adjoint{Bool,A},Transpose{Bool,A}}) where {N,A <: BitArray{N}}
    pB = parent(B)
    RowMajorStridedBitPointer(pointer(pB.chunks), tailstrides(pB))
end
@inline function stridedpointer(C::Union{Adjoint{T,A},Transpose{T,A}}) where {T, P, B, A <: SubArray{T,2,P,Tuple{Int,Vararg},B}}
    pC = parent(C)
    SparseStridedPointer(pointer(pC), staticmul(T, reverse(strides(pC))))
end

@inline stridedpointer(x::Number) = x
@inline stridedpointer(x::AbstractRange) = x
# @inline stridedpointer(ptr::Pointer) = PackedStridedPointer(pointer(ptr), tuple())
@inline stridedpointer(ptr::AbstractPointer) = ptr

@generated function noalias!(ptr::Ptr{T}) where {T <: NativeTypes}
    ptyp = JuliaPointerType
    typ = llvmtype(T)
    if Base.libllvm_version < v"10"
        funcname = "noalias" * typ
        decls = "define noalias $typ* @$(funcname)($typ *%a) willreturn noinline { ret $typ* %a }"
        instrs = """
            %ptr = inttoptr $ptyp %0 to $typ*
            %naptr = call $typ* @$(funcname)($typ* %ptr)
            %jptr = ptrtoint $typ* %naptr to $ptyp
            ret $ptyp %jptr
        """
    else
        decls = "declare void @llvm.assume(i1)"
        instrs = """
            %ptr = inttoptr $ptyp %0 to $typ*
            call void @llvm.assume(i1 true) ["noalias"($typ* %ptr)]
            %int = ptrtoint $typ* %ptr to $ptyp
            ret $ptyp %int
        """
    end
    quote
        $(Expr(:meta,:inline))
        llvmcall(
            $((decls, instrs)),
            Ptr{$T}, Tuple{Ptr{$T}}, ptr
        )
    end
end
@inline noalias!(ptr::AbstractPointer) = similar(ptr, noalias!(pointer(ptr)))
@inline noalias!(x::Any) = x
@inline noaliasstridedpointer(x) = stridedpointer(x)
@inline noaliasstridedpointer(A::AbstractRange) = stridedpointer(x)
@inline noaliasstridedpointer(A::AbstractArray) = noalias!(stridedpointer(A))
# @inline StaticStridedStruct{T,X}(s::S) where {T,X,S} = StaticStridedStruct{T,X,S}(s, 0)
# @inline StaticStridedStruct{T,X}(s::S, i::Int) where {T,X,S} = StaticStridedStruct{T,X,S}(s, i)
# @inline offset(ptr::StaticStridedStruct{T,X,S}, i::Integer) where {T,X,S} = StaticStridedStruct{T,X,S}(pointer(ptr), ptr.offset + i)
# # Trying to avoid generated functions
# @inline offset(ptr::StaticStridedStruct{T,X,S}, i::Tuple{<:Integer}) where {T,X,S} = @inbounds StaticStridedStruct{T,X,S}(pointer(ptr), ptr.offset + i[1])
# @inline function offset(ptr::StaticStridedStruct{T,Tuple{A},S}, i::Tuple{<:Integer,<:Integer}) where {T,A,S}
#     @inbounds StaticStridedStruct{T,Tuple{A},S}(pointer(ptr), ptr.offset + i[1] + A * i[2])
# end
# @inline function offset(ptr::StaticStridedStruct{T,Tuple{A,B},S}, i::Tuple{<:Integer,<:Integer,<:Integer}) where {T,A,B,S}
#     @inbounds StaticStridedStruct{T,Tuple{A,B},S}(pointer(ptr), ptr.offset + i[1] + A*i[2] + B*i[3])
# end
# @inline function offset(ptr::StaticStridedStruct{T,Tuple{A,B,C},S}, i::Tuple{<:Integer,<:Integer,<:Integer,<:Integer}) where {T,A,B,C,S}
#     @inbounds StaticStridedStruct{T,Tuple{A,B,C},S}(pointer(ptr), ptr.offset + i[1] + A*i[2] + B*i[3] + C*i[4] )
# end

# @generated tupletype_to_tuple(::Type{T}) where {T<:Tuple} = Expr(:block, Expr(:meta,:inline), Expr(:tuple, T.parameters...))
# @inline function offset(ptr::StaticStridedStruct{T,X,S}, i::NTuple{N}) where {T,X,S,N}
#     strides = tupletype_to_tuple(X)
#     StaticStridedStruct{T,X,S}(pointer(ptr), ptr.offset + first(i) + tdot(strides, Base.tail(i)))
# end

struct RangeWrapper{R <: AbstractRange, I}
    r::R
    i::I
end
@inline gesp(r::AbstractRange, i::I) where {I} = RangeWrapper(r, i)
@inline gesp(r::RangeWrapper, i::I) where {I} = RangeWrapper(r.r, r.i .+ i)
@inline vload(rw::RangeWrapper, i, mask) = vload(rw, i)
@inline vload(rw::RangeWrapper, i::Tuple{I}) where {I}  = vload(rw.r, @inbounds (vadd(rw.i[1], i[1]),))
@inline vload(rw::RangeWrapper, i::Tuple{I1,I2}) where {I1,I2}  = vload(rw.r, @inbounds (vadd(rw.i[1], i[1]), vadd(rw.i[2], i[2])))
@inline vload(rw::RangeWrapper, i::Tuple{I1,I2,I3}) where {I1,I2,I3}  = vload(rw.r, @inbounds (vadd(rw.i[1], i[1]), vadd(rw.i[2], i[2]), vadd(rw.i[3], i[3])))
@inline vload(rw::RangeWrapper, i)  = vload(rw.r, vadd.(rw.i, i))

@inline vload(r::AbstractRange, i::Tuple{I}) where {I} = @inbounds r[vadd(i[1], one(I))]
@inline vload(r::LinearIndices, i::Tuple) = @inbounds r[staticp1(i)...]

@inline subsetview(ptr::PackedStridedPointer, ::Val{1}, i::Integer) = SparseStridedPointer(gep(ptr.ptr, i), ptr.strides)
@generated function subsetview(ptr::PackedStridedPointer{T, N}, ::Val{I}, i::Integer) where {I, T, N}
    I > N + 1 && return :ptr
    strides = Expr(:tuple, [Expr(:ref, :s, n) for n ∈ 1:N if n != I-1]...)
    offset = Expr(:call, :vmul, :i, Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), Expr(:ref, :s, I - 1)))
    Expr(
        :block,
        Expr(:meta, :inline),
        Expr(:(=), :p, Expr(:(.), :ptr, QuoteNode(:ptr))),
        Expr(:(=), :s, Expr(:(.), :ptr, QuoteNode(:strides))),
        Expr(:(=), :strides, Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), strides)),
        Expr(:(=), :gp, Expr(:call, :gepbyte, :p, offset)),
        Expr(:call, :PackedStridedPointer, :gp, :strides)
    )
end

@generated function subsetview(ptr::RowMajorStridedPointer{T, N}, ::Val{I}, i::Integer) where {I, T, N}
    if N + 1 == I
        return Expr(
            :block, Expr(:meta, :inline),
            Expr(:(=), :p, Expr(:call, :gep, Expr(:(.), :ptr, QuoteNode(:ptr)), :i)),
            Expr(:call, :SparseStridedPointer, :p, Expr(:call, :reverse, Expr(:(.), :ptr, QuoteNode(:strides))))
        )
    elseif N + 1 < I
        return :ptr
    end
    strideind = N + 1 - I
    strides = Expr(:tuple, [Expr(:ref, :s, n) for n ∈ 1:N if n != strideind]...)
    offset = Expr(:call, :vmul, :i, Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), Expr(:ref, :s, strideind)))
    Expr(
        :block,
        Expr(:meta, :inline),
        Expr(:(=), :p, Expr(:(.), :ptr, QuoteNode(:ptr))),
        Expr(:(=), :s, Expr(:(.), :ptr, QuoteNode(:strides))),
        Expr(:(=), :strides, Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), strides)),
        Expr(:(=), :gp, Expr(:call, :gepbyte, :p, offset)),
        Expr(:call, :RowMajorStridedPointer, :gp, :strides)
    )
end

@generated function subsetview(ptr::SparseStridedPointer{T, N}, ::Val{I}, i::Integer) where {I, T, N}
    I > N && return :ptr
    strides = Expr(:tuple, [Expr(:ref, :s, n) for n ∈ 1:N if n != I]...)
    offset = Expr(:call, :vmul, :i, Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), Expr(:ref, :s, I)))
    Expr(
        :block,
        Expr(:meta, :inline),
        Expr(:(=), :p, Expr(:(.), :ptr, QuoteNode(:ptr))),
        Expr(:(=), :s, Expr(:(.), :ptr, QuoteNode(:strides))),
        Expr(:(=), :strides, Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), strides)),
        Expr(:(=), :gp, Expr(:call, :gepbyte, :p, offset)),
        Expr(:call, :SparseStridedPointer, :gp, :strides)
    )
end

@generated function subsetview(ptr::StaticStridedPointer{T, X}, ::Val{I}, i::Integer) where {I, T, X}
    I > length(X.parameters) && return :ptr
    Xa = Expr(:curly, :Tuple)
    Xparam = X.parameters
    for n ∈ 1:length(Xparam)
        n == I && continue
        push!(Xa.args, Xparam[n])
    end
    Xᵢ = (Xparam[I])::Int
    offset = Expr(:call, isone(Xᵢ) ? :vmulnp : :vmul, :i, Xᵢ * sizeof(T))
    Expr(
        :block,
        Expr(:meta, :inline),
        Expr(:(=), :p, Expr(:(.), :ptr, QuoteNode(:ptr))),
        Expr(:(=), :gp, Expr(:call, :gepbyte, :p, offset)),
        Expr(:call, Expr(:curly, :StaticStridedPointer, T, Xa), :gp)
    )
end

@generated function subsetview(ptr::StaticStridedStruct{T, X}, ::Val{I}, i::Integer) where {I, T, X}
    I > length(X.parameters) && return :ptr
    Xa = Expr(:curly, :Tuple)
    Xparam = X.parameters
    for n ∈ 1:length(Xparam)
        n == I && continue
        push!(Xa.args, Xparam[n])
    end
    offset = Expr(:call, :vmul, :i, Xparam[I])
    Expr(
        :block,
        Expr(:meta, :inline),
        Expr(:(=), :p, Expr(:(.), :ptr, QuoteNode(:ptr))),
        Expr(:(=), :offset, Expr(:call, :vadd, Expr(:(.), :ptr, QuoteNote(:offset)), offset)),
        Expr(:call, Expr(:curly, :StaticStridedStruct, T, Xa), :p, :offset)
    )
end

@inline stridedpointer(A::AbstractArray, indices::Tuple) = stridedpointer(view(A, indices...))
@inline stridedpointer(A::AbstractArray, ::Type{Transpose}) = stridedpointer(transpose(A))
@inline stridedpointer(A::AbstractArray, ::Type{Adjoint}) = stridedpointer(adjoint(A))
@inline stridedpointer(A::AbstractArray, ::Nothing) = stridedpointer(A)

@inline filter_strides_by_dimequal1(sz::NTuple{N,Int}, st::NTuple{N,Int}) where {N} = @inbounds ntuple(n -> isone(sz[n]) ? 0 : st[n], Val{N}())

@inline stridedpointer_for_broadcast(A::AbstractRange) = A
@inline function stridedpointer_for_broadcast(A::AbstractArray{T,N}) where {T,N}
    PackedStridedPointer(pointer(A), staticmul(T, filter_strides_by_dimequal1(Base.tail(size(A)), tailstrides(A))))
end
@inline stridedpointer_for_broadcast(B::Union{Adjoint{T,A},Transpose{T,A}}) where {T,A <: AbstractVector{T}} = stridedpointer_for_broadcast(parent(B))
@inline stridedpointer_for_broadcast(A::SubArray{T,0,P,S}) where {T,P,S <: Tuple{Int,Vararg}} = pointer(A)
@inline function stridedpointer_for_broadcast(A::SubArray{T,N,P,S}) where {T,N,P,S <: Tuple{Int,Vararg}}
    SparseStridedPointer(pointer(A), staticmul(T, filter_strides_by_dimequal1(size(A), strides(A))))
end
@inline function stridedpointer_for_broadcast(A::SubArray{T,N,P,S}) where {T,N,P,S}
    PackedStridedPointer(pointer(A), staticmul(T, filter_strides_by_dimequal1(Base.tail(size(A)), tailstrides(A))))
end
@inline stridedpointer_for_broadcast(A::BitArray) = stridedpointer(A)

struct MappedStridedPointer{F, T, P <: AbstractPointer{T}}
    f::F
    ptr::P
end
@inline vload(ptr::MappedStridedPointer) = ptr.f(vload(ptr.ptr))

"""
These use 0 based indexing for convenience with respect to the packed and row major strided pointers.
Preserves the first index, dropping the second.
"""
@inline function double_index(ptr::PackedStridedPointer{T}, ::Val{0}, ::Val{M}) where {M,T}
    x = ptr.strides
    SparseStridedPointer(pointer(ptr), @inbounds Base.setindex(x, vadd(x[M], sizeof(T)), M))
end
@generated function double_index(ptr::PackedStridedPointer{T,K}, ::Val{N}, ::Val{M}) where {K,N,M,T}
    tup = Expr(:tuple)
    for k ∈ 1:K
        k == N && continue
        if k == M
            push!(tup.args, Expr(:call, :vadd, Expr(:ref, :x, N), Expr(:ref, :x, M)))
        else
            push!(tup.args, Expr(:ref, :x, k))
        end
    end
    quote
        $(Expr(:meta,:inline))
        x = ptr.strides
        xabridged = @inbounds $tup
        PackedStridedPointer(pointer(ptr), xabridged)
    end
end
@inline function double_index(ptr::RowMajorStridedPointer{T,1}, ::Val{0}, ::Val{1}) where {T}
    x = ptr.strides
    @inbounds SparseStridedPointer(pointer(ptr), (vadd(x[1], sizeof(T)),))
end
@generated function double_index(ptr::SparseStridedPointer{T,K}, ::Val{Nm1}, ::Val{Mm1}) where {K,Nm1,Mm1, T}
    tup = Expr(:tuple)
    N = Nm1 + 1; M = Mm1 + 1
    for k ∈ 1:K
        k == N && continue
        if k == M
            push!(tup.args, Expr(:call, :vadd, Expr(:ref, :x, N), Expr(:ref, :x, M)))
        else
            push!(tup.args, Expr(:ref, :x, k))
        end
    end
    quote
        $(Expr(:meta,:inline))
        x = ptr.strides
        xabridged = @inbounds $tup
        SparseStridedPointer(pointer(ptr), xabridged)
    end
end
@generated function double_index(ptr::StaticStridedPointer{T,X}, ::Val{Nm1}, ::Val{Mm1}) where {X,Nm1,Mm1,T}
    tup = Expr(:curly, :Tuple)
    Xparam = X.parameters
    N = Nm1 + 1; M = Mm1 + 1
    for k ∈ eachindex(Xparam)
        k == N && continue
        if k == M
            push!(tup.args, (Xparam[N])::Int + (Xparam[M])::Int)
        else
            push!(tup.args, (Xparam[k])::Int)
        end
    end
    Expr(:block, Expr(:meta,:inline), :(StaticStridedPointer{$T,$tup}(pointer(ptr))))
end

struct ZeroInitializedStridedPointer{T, P <: AbstractStridedPointer{T}} <: AbstractStridedPointer{T}; ptr::P; end
@inline ZeroInitializedStridedPointer(A) = ZeroInitializedStridedPointer(stridedpointer(A))
@inline vload(::ZeroInitializedStridedPointer{T}) where {T} = zero(T)

@inline vec_length() = Val{1}()
@inline vec_length(::Integer) = Val{1}()
@inline vec_length(::Static) = Val{1}()
@inline vec_length(::_MM{W}) where {W} = Val{W}()
@inline vec_length(::Vec{W}) where {W} = Val{W}()
@inline vec_length(::SVec{W}) where {W} = Val{W}()
@inline promote_val_width(::Val{1}, ::Val{1}) = Val{1}()
@inline promote_val_width(::Val{W}, ::Val{1}) where {W} = Val{W}()
@inline promote_val_width(::Val{1}, ::Val{W}) where {W} = Val{W}()
@inline promote_val_width(::Val{W}, ::Val{W}) where {W} = Val{W}()

@inline vec_length(i::Tuple{}) = Val{1}()
@inline vec_length(i::Tuple{I1}) where {I1} = vec_length(first(i))
@inline vec_length(i::Tuple{I1,I2,Vararg}) where {I1,I2} = promote_val_width(vec_length(first(i)), vec_length(Base.tail(i)))

@inline vload(::ZeroInitializedStridedPointer{T}, ::Val{1}) where {T} = zero(T)
@inline vload(::ZeroInitializedStridedPointer{T}, ::Val{W}) where {T,W} = vzero(SVec{W,T})

@inline vload(::Val{W}, ::ZeroInitializedStridedPointer{T}, ::Any) where {W, T} = vzero(SVec{W,T})
@inline vload(::Val{W}, ::ZeroInitializedStridedPointer{T}, ::Any, ::Mask) where {W, T} = vzero(SVec{W,T})

@inline vload(ptr::ZeroInitializedStridedPointer, i::Tuple) = vload(ptr, vec_length(i))
@inline vload(ptr::ZeroInitializedStridedPointer, i::Tuple, ::Mask) = vload(ptr, vec_length(i))

@inline vstore!(ptr::ZeroInitializedStridedPointer, v::Number, i::Tuple) = vstore!(ptr.ptr, v, i)
@inline vstore!(ptr::ZeroInitializedStridedPointer, v::Number, i::Tuple, m::Mask) = vstore!(ptr.ptr, v, i, m)
@inline vstore!(ptr::ZeroInitializedStridedPointer, v::AbstractStructVec, i::Tuple) = vstore!(ptr.ptr, v, i)
@inline vstore!(ptr::ZeroInitializedStridedPointer, v::AbstractStructVec, i::Tuple, m::Mask) = vstore!(ptr.ptr, v, i, m)

@inline vnoaliasstore!(ptr::ZeroInitializedStridedPointer, v::Number, i::Tuple) = vnoaliasstore!(ptr.ptr, v, i)
@inline vnoaliasstore!(ptr::ZeroInitializedStridedPointer, v::Number, i::Tuple, m::Mask) = vnoaliasstore!(ptr.ptr, v, i, m)
@inline vnoaliasstore!(ptr::ZeroInitializedStridedPointer, v::AbstractStructVec, i::Tuple) = vnoaliasstore!(ptr.ptr, v, i)
@inline vnoaliasstore!(ptr::ZeroInitializedStridedPointer, v::AbstractStructVec, i::Tuple, m::Mask) = vnoaliasstore!(ptr.ptr, v, i, m)

@inline stride1offset(ptr::ZeroInitializedStridedPointer, i) = stride1offset(ptr.ptr, i)
@inline stridedoffset(ptr::ZeroInitializedStridedPointer, i) = stridedoffset(ptr.ptr, i)
@inline gep(ptr::ZeroInitializedStridedPointer, i::Tuple) = gep(ptr.ptr, i)
@inline gesp(ptr::ZeroInitializedStridedPointer, i::Tuple) = ZeroInitializedStridedPointer(gesp(ptr.ptr, i))

@inline Base.pointer(ptr::ZeroInitializedStridedPointer) = pointer(ptr.ptr)


text:
using ConstraintDomains
using Dictionaries
using Intervals
using Test

@testset "ConstraintDomains.jl" begin
    @testset "EmptyDomain" begin
        ed = domain()
        @test domain_size(ed) == 0 == length(ed)
        @test isempty(ed)
        @test π ∉ ed
    end

    @testset "DiscreteDomain" begin
        d1 = domain([4,3,2,1])
        d2 = domain(1)
        foreach(i -> add!(d2, i), 2:4)
        domains = Dictionary(1:2, [d1, d2])

        for d in domains
        # constructors and ∈
            for x in [1,2,3,4]
                @test x ∈ d
            end
        # length
            @test length(d) == 4
        # draw and ∈
            @test rand(d) ∈ d
        # add!
            add!(d, 5)
            @test 5 ∈ d
        # delete!
            delete!(d, 5)
            @test 5 ∉ d
            @test domain_size(d) == 3
        end
    end

    @testset "ContinuousDomain" begin
        d1 = domain(1.0..3.15)
        d2 = domain(Interval{Open, Open}(-42.42, 5.0))
        # d3 = domain([d1, d2])
        domains = Dictionary(1:2, [d1, d2])

        for d in domains
            for x in [1, 2.3, π]
                @test x ∈ d
            end
            for x in [5.1, π^π, Inf]
                @test x ∉ d
            end
            @test rand(d) ∈ d
            @test rand(d, 1) ∈ d
            @test domain_size(d) > 0.0
        end

    end

    @testset "RangeDomain" begin
        d1 = domain(1:5)
        d2 = domain(1:.5:5)
        domains = Dictionary(1:2, [d1, d2])

        for d in domains
            for x in [1, 2, 3, 4, 5]
                @test x ∈ d
            end
            for x in [42]
                @test x ∉ d
            end
            @test rand(d) ∈ d
        end

    end
end