text:
using ParallelAccelerator

@acc function opt_update_w( k, m, W, Y, E1 )
	 for j = 1:k
            #w = view(W,:,j)
            #y = view(Y,:,j)
            e1 = E1[j]
            for i = 1:m
                #w[i] = y[i] - e1 * w[i]
		W[i,j] = Y[i,j] - e1 * W[i,j]
            end
        end
end

function main(m::Int, k::Int)
    Y  = Array(Float64, m, k)    # to store E{x g(w'x)} for components
    E1 = Array(Float64, k)       # store E{g'(w'x)} for components

    W  = Array(Float64, m, k)    # to store E{x g(w'x)} for components

    fill!(W, 3)
    fill!(Y, 9)
    fill!(E1, 2)
    opt_update_w(k,m,W,Y,E1)

    return W

end

main(100, 100)
println("Done testing test_update_w.")


text:
"""
Li, Yuhua, 2011.
Selecting training points for one-class support vector machines.
Pattern recognition letters 32.11: 1517-1522.

We dub the method boundary point selector (BPS).

Parameter guideline from the paper, given a data set of size N:
k = floor(10 ln N)
eps ∈ [0, 0.2]
"""
mutable struct BPS <: Sampler
    k::Union{Int, Nothing}
    eps::Real
    function BPS(k::Union{Int, Nothing}=nothing, eps::Real=0.2)
        isa(k, Int) && (k > 0 || throw(ArgumentError("Invalid number of neighbors to consider $(k).")))
        eps >= 0 || throw(ArgumentError("Invalid threshold $(eps)."))
        return new(k, eps)
    end
end

function sample(sampler::BPS, data::DataSet, labels::Vector{Symbol})::BitArray{1}
    if sampler.k === nothing
        sampler.k = round(Int, floor(10 * log(size(data, 2))))
    end
    sampler.k = min(size(data, 2) - 1, sampler.k)
    tree = KDTree(data)
    sample_mask = falses(size(data, 2))

    for i in 1:size(data, 2)
        idx, dist = knn(tree, data[:, i], sampler.k + 1, true)
        v_ij = mapslices(normalize, data[:, i] .- data[:, idx[2:end]], dims=1)
        n_i = sum(v_ij, dims=2)
        θ_ij = sum(v_ij .* n_i, dims=1)
        l_i = 1 / sampler.k * sum(θ_ij .>= 0)
        if l_i >= 1 - sampler.eps
            sample_mask[i] = true
        end
    end
    return sample_mask
end

ext: jl


text:
module FinancePlus

export simple_interest,compound_interest,future_value,present_value,present_annuity,present_annuity_due,present_value_interest_factor,future_value_interest_factor
# Write your package code here.
include("TimeValue.jl")

si = simple_interest(100.00 , 0.1 ,2)
println(si)

ci = compound_interest(100.00 , 0.1 ,2)
println(ci)

ann1 = present_annuity(50000.00 , 0.06 ,25)
println(ann1)

ann2 = present_annuity_due(50000.00 , 0.06 ,25)
println(ann2)

PVIF= present_value_interest_factor(0.05 ,5)
println(PVIF)

FVIF = future_value_interest_factor(0.05 ,5)
println(FVIF)
end


text:
include("regression.jl")
include("classification.jl")
include("multiclass.jl")
include("event.jl")

function pdf(l::Likelihood{T},y::Real,f::Real) where T
    @error "pdf not implemented"
end

function logpdf(l::Likelihood{T},y::Real,f::Real) where T
    log(pdf(l,y,f))
end

isaugmented(::Likelihood{T}) where T = false


text:
# https://adventofcode.com/2020/day/20

function parse_input(input)
    tiles = Dict()
    for chunk in split(input, "\n\n", keepempty=false)
        rows = split(chunk, '\n', keepempty=false)
        id = parse(Int, match(r"Tile (\d+):", popfirst!(rows)).captures[1])
        tile = zeros(Bool, (length(rows), length(rows[2])))
        for (i, row) in enumerate(rows)
            tile[i, :] .= collect(row) .== '#'
        end
        tiles[id] = tile
    end
    return tiles
end

"""
Sides:
```
     1
   +---+
 4 |   | 2
   +---+
     3
```
"""
function get_side(tile, side)
    if side == 1
        return tile[1, :][:]
    elseif side == 2
        return tile[:, end][:]
    elseif side == 3
        return tile[end, :][:]
    else
        return tile[:, 1][:]
    end
end

function get_sides(tiles)
    sides = Dict()
    for (id, tile) in tiles
        sides[id] = [get_side(tile, i) for i in 1:4]
    end
    return sides
end

function find_matches(tiles)
    sides = get_sides(tiles)
    matches = Dict()
    for (x_id, x_sides) in sides, (y_id, y_sides) in sides
        if x_id == y_id
            continue
        end

        for i in 1:4, j in 1:4
            rotated = all(x_sides[i] .== reverse(y_sides[j]))
            if all(x_sides[i] .== y_sides[j]) || rotated
                match = (y_id, i, j, rotated)
                if x_id in keys(matches)
                    push!(matches[x_id], match)
                else
                    matches[x_id] = [match]
                end
            end
        end
    end
    return matches
end

function part1(input)
    tiles = parse_input(input)
    matches = find_matches(tiles)
    corners = [k for (k, v) in matches if length(v) == 2]
    @assert length(corners) == 4
    return (*)(corners...)
end

example = read("data/day-20-example.txt", String)

@assert part1(example) == 20899048083289

function flip(arr::Array{T,2}; dims=1) where {T}
    n, k = size(arr)
    out = Array{T,2}(undef, n, k)
    if dims == 1
        for i in 1:n
            out[n - i + 1, :] = arr[i, :]
        end
    elseif dims == 2
        for j in 1:k
            out[:, k - j + 1, :] = arr[:, j]
        end
    else
        return flip(flip(arr, dims=1), dims=2)
    end
    return out
end

@assert flip([1 2; 3 4], dims=1) == [3 4; 1 2]
@assert flip([1 2; 3 4], dims=2) == [2 1; 4 3]
@assert flip([1 2; 3 4], dims=(1, 2)) == [4 3; 2 1]

"""
```
     1   1
   +---+---+
 4 |   |   | 2
   +---+---+
 4 |   |   | 2
   +---+---+
     3   3
```
"""
function collect_image(matches, tiles)
    n = Int(sqrt(length(matches)))
    grid = Array{Any,2}(undef, n, n)
    tile_id = nothing

    for (id, paired) in matches
        sides = sort([x[2] for x in paired])
        if length(sides) == 2 && all(sides .== [2, 3])
            grid[1, 1] = copy(tiles[id])
            tile_id = id
            break
        end
    end

    return grid
end

function brute_force_match(pattern, matches, tiles)
    for (id, _) in matches
        candidate = tiles[id]
        for transform in [
            identity,
            x -> flip(x, dims=1),
            x -> flip(x, dims=2),
            x -> flip(x, dims=(1, 2)),
            x -> collect(x'),
            x -> flip(collect(x'), dims=1),
            x -> flip(collect(x'), dims=2),
            x -> flip(collect(x'), dims=(1, 2)),
        ]
            transformed = transform(candidate)
            if all(pattern .== get_side(transformed, 4))
                return id, transformed
            end
        end
    end
end

test = read("data/day-20.txt", String)
println("Part 1: $(result1 = part1(test))")
# println("Part 2: $(result2 = part2(test))")

@assert result1 == 107399567124539
# @assert result2 ==

max_stars_repo_name: twolodzko/advent-of-code

max_issues_repo_name: twolodzko/advent-of-code

max_forks_repo_name: twolodzko/advent-of-code


text:
import JSON
import Plots

result = JSON.parsefile("build/step_response.json")

Plots.gr()

p = Plots.plot()

Plots.xlabel!(p, "t (s)")
Plots.ylabel!(p, "value")
Plots.plot!(p, [itr[1] for itr in result], [itr[2] for itr in result], lab="input")
Plots.plot!(p, [itr[1] for itr in result], [itr[3] for itr in result], lab="output")
Plots.savefig(p, "build/step_response.pdf")

max_stars_repo_name: Shimushushushu/etude-naive-controller

max_issues_repo_name: Shimushushushu/etude-naive-controller

max_forks_repo_name: Shimushushushu/etude-naive-controller


text:
function day6(N)
    nums = parse.(Int, split(readline("../data/input6.txt"), ","))
    time = zeros(Int, 9)
    @inbounds for num in nums
        time[num+1] += 1
    end

    @inbounds for _ in 1:N
        time = circshift(time, -1)
        time[7] += time[end]
    end
    return sum(time)
end

nums = parse.(Int, split(readline("../data/input6.txt"), ","))

using StaticArrays
function flash_day6(N, nums)
    M = @SMatrix [0 1 0 0 0 0 0 0 0
         0 0 1 0 0 0 0 0 0
         0 0 0 1 0 0 0 0 0
         0 0 0 0 1 0 0 0 0
         0 0 0 0 0 1 0 0 0
         0 0 0 0 0 0 1 0 0
         1 0 0 0 0 0 0 1 0
         0 0 0 0 0 0 0 0 1
         1 0 0 0 0 0 0 0 0]

    time = zeros(Int, 9)
    @inbounds for num in nums
        time[num] += 1
    end
    return sum(M^N*time)
end

function part1(nums)
    M = @SMatrix [ 252   20  210   37  120   84   45  126   11
                    56  252   20  210   37  120   84   45  126
                    210   56  252   20  210   37  120   84   45
                    165  210   56  252   20  210   37  120   84
                    121  165  210   56  252   20  210   37  120
                    330  121  165  210   56  252   20  210   37
                    57  330  121  165  210   56  252   20  210
                    210   37  120   84   45  126   11  126    9
                    20  210   37  120   84   45  126   11  126]

    time = zeros(Int, 9)
    @inbounds for num in nums
        time[num+1] += 1
    end
    return sum(M*time)
end

function part2(nums)
    M = @SMatrix [  496266131  589731885  399865906  491122368  357868865
                    655568076  496266131  589731885  399865906  491122368
                    659462321  655568076  496266131  589731885  399865906
                    697451775  659462321  655568076  496266131  589731885
                    869885915  697451775  659462321  655568076  496266131
                    757734771  869885915  697451775  659462321  655568076
                    1080854253  757734771  869885915  697451775  659462321
                    399865906  491122368  357868865  378763547  339582910
                    589731885  399865906  491122368  357868865  378763547]

    time = zeros(Int, 5)
    @inbounds for num in nums
        time[num] += 1
    end
    return sum(M*time)
end

function flash_day6_part2(nums)
    magic = (6206821033,  5617089148,  5217223242,  4726100874,  4368232009)
    return sum(magic[i] for i in nums)
end

macro generate_fast_solution(N)
    M = [0 1 0 0 0 0 0 0 0
        0 0 1 0 0 0 0 0 0
        0 0 0 1 0 0 0 0 0
        0 0 0 0 1 0 0 0 0
        0 0 0 0 0 1 0 0 0
        0 0 0 0 0 0 1 0 0
        1 0 0 0 0 0 0 1 0
        0 0 0 0 0 0 0 0 1
        1 0 0 0 0 0 0 0 0]

    magic = sum(M ^ N; dims=1)
    magic = (magic[2], magic[3], magic[4], magic[5], magic[6])
    name = Symbol("day6_", N, "_days")
    code = :($name(nums) = (magic = $magic; sum(magic[i] for i in nums)))
    return esc(code)
end

max_stars_repo_name: lucaferranti/AoC

max_issues_repo_name: lucaferranti/AoC

max_forks_repo_name: lucaferranti/AoC


text:
module HyperParametersModule
import Base: convert, eltype, promote_type, show, string, ==, *, /, +, -, ^, isless;
using GradDescent
export
    Bound,
        LeftBound,
        RightBound,
        NullBound,

    Interval,
    interval,
    checkbounds,

    HyperParameter,
    HyperParameters,
    update!,
    getvalue,
    setvalue!,
    checkvalue,
    gettheta,
    settheta!,
    checktheta,
    upperboundtheta,
    lowerboundtheta,
    setfixed!,
    setfree!,
    setparamoptimizer!

include("bound.jl")
include("interval.jl")
include("hyperparameter.jl")

end # End HyperParameter


text:
"""
$(SIGNATURES)

Takes a html string that may contain inline katex blocks `\\(...\\)` or display katex blocks
`\\[ ... \\]` and use node and katex to pre-render them to HTML.
"""
function js_prerender_katex(hs::String)::String
    # look for \(, \) and \[, \] (we know they're paired because generated from markdown parsing)
    matches = collect(eachmatch(r"\\(\(|\[|\]|\))", hs))
    isempty(matches) && return hs

    # buffer to write the JS script
    jsbuffer = IOBuffer()
    write(jsbuffer, """
            const katex = require("$(escape_string(joinpath(PATHS[:libs], "katex", "katex.min.js")))");
            """)
    # string to separate the output of the different blocks
    splitter = "_>fdsplit<_"

    # go over each match and add the content to the jsbuffer
    for i ∈ 1:2:length(matches)-1
        # tokens are paired, no nesting
        mo, mc = matches[i:i+1]
        # check if it's a display style
        display = (mo.match == "\\[")
        # this is the content without the \( \) or \[ \]
        ms = subs(hs, nextind(hs, mo.offset + 1), prevind(hs, mc.offset))
        # add to content of jsbuffer
        write(jsbuffer, """
            var html = katex.renderToString("$(escape_string(ms))", {displayMode: $display})
            console.log(html)
            """)
        # in between every block, write $splitter so that output can be split easily
        i == length(matches)-1 || write(jsbuffer, """\nconsole.log("$splitter")\n""")
    end
    return js2html(hs, jsbuffer, matches, splitter)
end

"""
$(SIGNATURES)

Takes a html string that may contain `<pre><code ... </code></pre>` blocks and use node and
highlight.js to pre-render them to HTML.
"""
function js_prerender_highlight(hs::String)::String
    # look for "
<code" and "
" these will have been automatically generated
    # and therefore the regex can be fairly strict with spaces etc
    matches = collect(eachmatch(r"
|<\/code><\/pre>", hs))
    isempty(matches) && return hs

    # buffer to write the JS script
    jsbuffer = IOBuffer()
    write(jsbuffer, """const hljs = require('$(HIGHLIGHTJS[])');""")

    # string to separate the output of the different blocks
    splitter = "_>fdsplit<_"

    # go over each match and add the content to the jsbuffer
    for i ∈ 1:2:length(matches)-1
        # tokens are paired, no nesting
        co, cc = matches[i:i+1]
        # core code
        cs = subs(hs, nextind(hs, matchrange(co).stop), prevind(hs, matchrange(cc).start))
        # cs = escape_string(cs)
        # lang
        lang = co.captures[2]
        # un-escape code string
        cs = html_unescape(cs) |> escape_string
        # add to content of jsbuffer
        write(jsbuffer, """console.log("
" + hljs.highlight("$cs", {'language': "$lang"}).value + "
");""")
        # in between every block, write $splitter so that output can be split easily
        i == length(matches)-1 || write(jsbuffer, """console.log('$splitter');""")
    end
    return js2html(hs, jsbuffer, matches, splitter)
end

"""
$(SIGNATURES)

Convenience function to run the content of `jsbuffer` with node, and lace the results with `hs`.
"""
function js2html(hs::String, jsbuffer::IOBuffer, matches::Vector{RegexMatch},
                 splitter::String)::String
    # run it redirecting the output to a buffer
    outbuffer = IOBuffer()
    run(pipeline(`$NODE -e "$(String(take!(jsbuffer)))"`, stdout=outbuffer))

    # read the buffer and split it using $splitter
    out = String(take!(outbuffer))
    parts = split(out, splitter)

    # lace everything back together
    htmls = IOBuffer()
    head, c = 1, 1
    for i ∈ 1:2:length(matches)-1
        mo, mc = matches[i:i+1]
        write(htmls, subs(hs, head, prevind(hs, mo.offset)))
        pp = strip(parts[c])
        if startswith(pp, "
"shell>")
            pp = replace(pp, r"(\(.*?\)) pkg>"=>s"\1 pkg>")
        end
        write(htmls, pp)
        head = mc.offset + length(mc.match)
        c += 1
    end
    # add the rest of the document beyond the last mathblock
    head < lastindex(hs) && write(htmls, subs(hs, head, lastindex(hs)))

    return String(take!(htmls))
end


text:
include("utils.jl")

@testset "Tangent bundle" begin
    M = Sphere(2)

    @testset "FVector" begin
        tvs = ([1.0, 0.0, 0.0], [0.0, 1.0, 0.0])
        fv_tvs = map(v -> FVector(TangentSpace, v), tvs)
        fv1 = fv_tvs[1]
        tv1s = similar(fv_tvs[1])
        @test isa(tv1s, FVector)
        @test tv1s.type == TangentSpace
        @test size(tv1s.data) == size(tvs[1])
        @test eltype(tv1s) == eltype(tvs[1])
        @test isa(fv1 + fv1, FVector)
        @test (fv1 + fv1).type == TangentSpace
        @test isa(fv1 - fv1, FVector)
        @test (fv1 - fv1).type == TangentSpace
        @test isa(-fv1, FVector)
        @test (-fv1).type == TangentSpace
        @test isa(2*fv1, FVector)
        @test (2*fv1).type == TangentSpace
    end

    types = [Vector{Float64},
             SizedVector{3, Float64},
             MVector{3, Float64},
             Vector{Float32},
             SizedVector{3, Float32},
             MVector{3, Float32},
             Vector{Double64},
             MVector{3, Double64},
             SizedVector{3, Double64}]
    for T in types
        x = convert(T, [1.0, 0.0, 0.0])
        TB = TangentBundle(M)
        @test base_manifold(TB) == M
        @testset "Type $T" begin
            pts_tb = [ProductRepr(convert(T, [1.0, 0.0, 0.0]), convert(T, [0.0, -1.0, -1.0])),
                      ProductRepr(convert(T, [0.0, 1.0, 0.0]), convert(T, [2.0, 0.0, 1.0])),
                      ProductRepr(convert(T, [1.0, 0.0, 0.0]), convert(T, [0.0, 2.0, -1.0]))]
            @inferred ProductRepr(convert(T, [1.0, 0.0, 0.0]), convert(T, [0.0, -1.0, -1.0]))
            for pt ∈ pts_tb
                @test bundle_projection(TB, pt) ≈ pt.parts[1]
            end
            test_manifold(TB,
                          pts_tb,
                          test_reverse_diff = isa(T, Vector),
                          test_tangent_vector_broadcasting = false,
                          test_project_tangent = true)
        end
    end

    @test TangentBundle{Sphere{2}} == VectorBundle{Manifolds.TangentSpaceType, Sphere{2}}
    @test CotangentBundle{Sphere{2}} == VectorBundle{Manifolds.CotangentSpaceType, Sphere{2}}

    @test base_manifold(TangentBundle(M)) == M
    @testset "spaces at point" begin
        x = [1.0, 0.0, 0.0]
        t_x = TangentSpaceAtPoint(M, x)
        ct_x = CotangentSpaceAtPoint(M, x)
        @test base_manifold(t_x) == M
        @test base_manifold(ct_x) == M
        @test t_x.fiber.M == M
        @test ct_x.fiber.M == M
        @test t_x.fiber.VS == TangentSpace
        @test ct_x.fiber.VS == CotangentSpace
        @test t_x.x == x
        @test ct_x.x == x
    end

    @testset "tensor product" begin
        TT = Manifolds.TensorProductType(TangentSpace, TangentSpace)
        @test vector_space_dimension(VectorBundleFibers(TT, Sphere(2))) == 4
        @test vector_space_dimension(VectorBundleFibers(TT, Sphere(3))) == 9
    end

end


text:
using HypothesisTests, Test
using HypothesisTests: default_tail

@testset "Jarque-Bera" begin
sim_data_h0 = 1.0 .+ [
    0.2972879845354616, 0.3823959677906078, -0.5976344767282311, -0.01044524463737564,
    -0.839026854388764, 0.31111133849833383, 2.2950878238373105, -2.2670863488005306,
    0.5299655761667461, 0.43142152642291204, 0.5837082875687786, 0.9632716050381906,
    0.45879095505371686, -0.5223367574215084, 0.40839583832475224, -0.050451229933665284,
    -0.6936536438038856, -1.7738332534080556, 0.12076596493743928, -0.7576330377205277,
    -1.7297561813906863, 0.7959486220046159, 0.6700619812560624, 0.5508523775007609,
    -0.06337459242956062, 1.3369427822509223, -0.07314863333773934, -0.7454643166407556,
    -1.2200551285346526, -0.05317733752769253, -0.1651363578028337, -2.115370248427311,
    -0.0667679865065892, 1.2224618484680618, 0.5676954597384435, 0.7634989283233873,
    0.37859887985335006, -0.645597056636597, -0.6646462443910314, -1.8030263242108941,
    -1.2890476038497443, -0.3352700795181175, 0.07046760136115796, 0.34179376737982303,
    1.7351656660543822, 1.299915964310192, 0.20636408140261675, -1.0088599555376034,
    -0.8500556703153979, 1.1294094358843267, -1.348533456491827, 1.2402909004290597,
    -0.055647093207324264, 0.7593902139454598, -0.030266886646762546, -0.36108266775609366,
    -2.0073032025325626, -0.5688591273466379, -1.1496275243969976, 1.880622174631185,
    -1.4706245413873802, -0.5352109701657319, -0.9635439598255011, -1.3851080676759138,
    0.13430023756127, -0.6161166660088192, -1.719989356319742, 0.32076927998741056,
    -1.4473728978815112, -0.4922711880695664, -0.4130873839831828, 1.001928512455953,
    1.174681325975111, -1.326970345062047, -1.593214868044067, -0.43870997575879667,
    1.225797497997867, 1.4890902780696147, 2.1594331871008534, -0.9094790180375028,
    1.1478203006125018, -0.5050721965223731, -0.2670670352181297, 1.499082393930729,
    0.7973041601753904, -0.17164044330709693, -0.46907960993009634, 0.21762402120860017,
    0.35914640760234967, 0.3200247157080577, 0.2592163280556861, 0.459695985041628,
    0.20998636004214272, 1.5231753911772645, 1.2597946036121963, -0.16528310111426064,
    0.15619449488763584, -0.8143152943766987, -1.7098682060045913, 1.4999378761140325,
    -0.1289156385300165, -0.4501839715598432
]

t = JarqueBeraTest(sim_data_h0)

@test t.n == 102
@test t.JB ≈ 1.0204197333021678
@test t.skew ≈ -0.020527653857777352
@test t.kurt ≈ 2.5117242352057993
@test pvalue(t) ≈ 0.6003695680393418
@test default_tail(t) == :right
show(IOBuffer(), t)

sim_data_h1 = [
    0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1,
    0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0,
    1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1,
    1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1
]

t = JarqueBeraTest(sim_data_h1)

@test t.n == 102
@test t.JB ≈ 17.00081983024691
@test t.skew ≈ 0.11785113019775637
@test t.kurt ≈ 1.0138888888888888
@test pvalue(t) ≈ 0.00020338498134114293
@test default_tail(t) == :right
show(IOBuffer(), t)
end


text:
using Statistics

"""
    mse(y, t; back = false)
Return the mean square error:
    mean((y .- t) .^ 2)
# Example
```
julia> loss = mse
mse (generic function with 1 method)

julia> y, t = [1, 3, 2, 5], [1, 2, 3, 4]
([1, 3, 2, 5], [1, 2, 3, 4])

julia> loss(y, t)
0.75
```
"""
function mse(y, t; back = false)
    if back == true
        return 2(y .- t)
    end
    return mean((y .- t) .^ 2)
end


text:
#   Copyright 2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
The TestTools package provides support for testing and code quality CLIs.
"""
module TestTools

# --- Package Metadata

using TOML: TOML
const VERSION = TOML.parsefile(joinpath(pkgdir(@__MODULE__), "Project.toml"))["version"]

# --- Submodules

# jlcodestyle
include("jlcodestyle/jlcodestyle.jl")

# jlcoverage
include("jlcoverage/jlcoverage.jl")

# jltest
include("jltest/jltest.jl")

# --- Package management functions

include("pkg.jl")

end  # End of TestTools module


text:
#=
control of Lotka-Volterra population dynamics
adapted from
"Moment-sum-of-squares hierarchies for set approximation and optimal control"
by Milan Korda (2016)

TODO
- add options to use standard PSD cone formulation vs WSOS cone formulation
- add random data generation
=#

import SpecialFunctions: gamma
import DynamicPolynomials
const DP = DynamicPolynomials
import SemialgebraicSets
const SAS = SemialgebraicSets
import SumOfSquares
import PolyJuMP

struct LotkaVolterraJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    deg::Int # polynomial degrees
end

function build(inst::LotkaVolterraJuMP{T}) where {T <: Float64}
    # parameters
    n = 4 # number of species
    m = 2 * n # number of control inputs (u)
    Q = 0.475
    q = 0.525
    l_x = 1.0 # cost of rho
    l_u = [-1.0, 0.5, 0.6, 0.8, 1.1, 2.0, 4.0, 6.0] # cost of sigmas
    r = [1.0, 0.6, 0.4, 0.2] # growth rate of species

    DP.@polyvar x_h[1:n]
    x_mon = DP.monomials(x_h, 0:inst.deg)
    x_o = x_h * Q .+ q
    A = [1.0 0.3 0.4 0.2; -0.2 1.0 0.4 -0.1;
        -0.1 -0.2 1.0 0.3; -0.1 -0.2 -0.3 1.0]
    M = (sum(abs, l_u) + sum(l_u)) / 2.0 + l_x
    M *= 0.01 # upper bound on the total cost
    f = r .* x_o .* (1.0 .- A * x_o)
    f_u = hcat(Matrix(-1.0I, n, n), Matrix(1.0I, n, n))
    brho = 1.0
    u_bar = 1.0 # upper bound on u
    X = SAS.@set x_h' * x_h <= 1.0 # non-extinction domain
    delta_X = SAS.@set x_h' * x_h == 1.0

    function integrate_ball_monomial(mon, n)
        as = DP.exponents(mon)
        @assert length(as) == n
        if any(isodd, as)
            return 0.0
        else
            bs = (as .+ 1) ./ 2.0
            return 2.0 * prod(gamma.(bs)) / (gamma(sum(bs)) * (sum(as) + n))
        end
    end

    integrate_ball(p, n) = sum(DP.coefficient(t) * integrate_ball_monomial(t, n)
        for t in DP.terms(p))

    model = SumOfSquares.SOSModel()
    JuMP.@variable(model, rho, PolyJuMP.Poly(x_mon))
    JuMP.@variable(model, rho_T, PolyJuMP.Poly(x_mon))
    JuMP.@variable(model, sigma[1:m], PolyJuMP.Poly(x_mon))

    JuMP.@objective(model, Min, integrate_ball(l_x * rho, n) +
        sum(integrate_ball(sigma[i] * l_u[i], n) for i in 1:m) +
        M * integrate_ball(rho_T, n))

    JuMP.@constraint(model, rho <= 0, domain = delta_X)
    JuMP.@constraint(model, rho_T + brho * rho +
        sum(DP.differentiate(rho * f[i], x_h[i]) / Q for i in 1:n) +
        sum(sum(DP.differentiate(sigma[j] * f_u[i, j], x_h[i]) / Q
            for i in 1:n) for j in 1:m) >= 1, domain = X)
    JuMP.@constraint(model, [i in 1:m], u_bar * rho >= sigma[i], domain = X)
    JuMP.@constraint(model, rho_T >= 0, domain = X)
    JuMP.@constraint(model, [i in 1:m], sigma[i] >= 0, domain = X)

    return model
end


text:
module ControlBenchmarks

################################################################################
# Main exports from this package
#
# The individual benchmarks are exported at the top of their respective files.
################################################################################
export controlbenchmark,
       parameterbounds,
       listbenchmarks,
       getbenchmark

################################################################################
# Packages we need
################################################################################
using ControlSystems
using ControlSystems: AbstractSystem, nstates, ninputs
using LazySets
using LinearAlgebra
using Printf
using PrettyTables
using Polyhedra
using SparseArrays

# Bounds for a variable/variables
"""
    Bound = LazySets.Interval

Alias to allow easy use of a `LazySets.Interval` to represent the upper and lower bound for a single
variable.
"""
const Bound = LazySets.Interval

"""
    Bounds = Vector{Bound}

Alias to allows easy use of a collection of single-variable upper and lower bounds.
"""
const Bounds = Vector{Bound}

"""
    ControlBenchmarkProblem

The abstract type that all benchmark problems must derive from.
"""
abstract type ControlBenchmarkProblem end

include( "internals.jl" )
include( "interface.jl" )

################################################################################
# The benchmark systems
################################################################################

# The IFAC benchmark systems
include( "benchmarks/ifac/prob9001-binaryDistilationColumn.jl" )
include( "benchmarks/ifac/prob9002-drumBoiler.jl" )
include( "benchmarks/ifac/prob9006-airplaneFlutter.jl" )
include( "benchmarks/ifac/prob9007-hydraulicPositioningSystem.jl" )

# Mass-Spring systems
include( "benchmarks/massSpringSystems/linearMassSpringDamper.jl" )

# Miscellaneous benchmark systems
include( "benchmarks/ballOnPlate.jl" )
include( "benchmarks/jonesMorari.jl" )

end # module


text:
function compare_to(obj::Enum, arg0::Enum)
    return jcall(obj, "compareTo", jint, (Enum,), arg0)
end

function describe_constable(obj::Enum)
    return jcall(obj, "describeConstable", Optional, ())
end

function equals(obj::Enum, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_declaring_class(obj::Enum)
    return jcall(obj, "getDeclaringClass", Class, ())
end

function get_localized_string(obj::LocalizedOptimFormats, arg0::Locale)
    return jcall(obj, "getLocalizedString", JString, (Locale,), arg0)
end

function get_source_string(obj::LocalizedOptimFormats)
    return jcall(obj, "getSourceString", JString, ())
end

function hash_code(obj::Enum)
    return jcall(obj, "hashCode", jint, ())
end

function name(obj::Enum)
    return jcall(obj, "name", JString, ())
end

function ordinal(obj::Enum)
    return jcall(obj, "ordinal", jint, ())
end

function to_string(obj::Enum)
    return jcall(obj, "toString", JString, ())
end

function value_of(::Type{Enum}, arg0::Class, arg1::JString)
    return jcall(Enum, "valueOf", Enum, (Class, JString), arg0, arg1)
end

function value_of(::Type{LocalizedOptimFormats}, arg0::JString)
    return jcall(LocalizedOptimFormats, "valueOf", LocalizedOptimFormats, (JString,), arg0)
end

function values(::Type{LocalizedOptimFormats})
    return jcall(LocalizedOptimFormats, "values", Vector{LocalizedOptimFormats}, ())
end


text:
import 'lib_src/print.jl'

def main {
    print_strln('Hello, world!')
}

max_stars_repo_name: jeffreymanzione/jeff-vm-lang

max_issues_repo_name: jeffreymanzione/jeff-vm-lang

max_forks_repo_name: jeffreymanzione/jeff-vm-lang


text:
@testset "occupancy" begin

let
    md = CuModuleFile(joinpath(@__DIR__, "ptx/dummy.ptx"))
    dummy = CuFunction(md, "dummy")

    active_blocks(dummy, 1)
    active_blocks(dummy, 1; shmem=64)

    occupancy(dummy, 1)
    occupancy(dummy, 1; shmem=64)

    launch_configuration(dummy)
    launch_configuration(dummy; shmem=64)
    launch_configuration(dummy; shmem=64, max_threads=64)

    cb_calls = 0
    launch_configuration(dummy; shmem=threads->(cb_calls += 1; 0))
    @test cb_calls > 0
end

end


text:
include("anderson.jl")
include("problem_types.jl")


text:
"""
    shift_constant(set::MOI.AbstractScalarSet,
                   offset)

Returns a new scalar set `new_set` such that `func`-in-`set` is equivalent to
`func + offset`-in-`new_set`.

## Examples

The call `shift_constant(MOI.Interval(-2, 3), 1)` is equal to
`MOI.Interval(-1, 4)`.
"""
function shift_constant(set::Union{MOI.LessThan{T},
                                   MOI.GreaterThan{T},
                                   MOI.EqualTo{T}},
                        offset::T) where T
    return typeof(set)(MOI.constant(set) + offset)
end
function shift_constant(set::MOI.Interval, offset)
    return MOI.Interval(set.lower + offset, set.upper + offset)
end

const ScalarLinearSet{T} = Union{MOI.EqualTo{T}, MOI.LessThan{T}, MOI.GreaterThan{T}}
const VectorLinearSet = Union{MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives}

vector_set_type(::Type{<:MOI.EqualTo}) = MOI.Zeros
vector_set_type(::Type{<:MOI.LessThan}) = MOI.Nonpositives
vector_set_type(::Type{<:MOI.GreaterThan}) = MOI.Nonnegatives

scalar_set_type(::Type{<:MOI.Zeros}, T::Type) = MOI.EqualTo{T}
scalar_set_type(::Type{<:MOI.Nonpositives}, T::Type) = MOI.LessThan{T}
scalar_set_type(::Type{<:MOI.Nonnegatives}, T::Type) = MOI.GreaterThan{T}

"""
    is_diagonal_vectorized_index(index::Base.Integer)

Return whether `index` is the index of a diagonal element in a
[`MOI.AbstractSymmetricMatrixSetTriangle`](@ref) set.
"""
function is_diagonal_vectorized_index(index::Base.Integer)
    perfect_square = 1 + 8index
    return isqrt(perfect_square)^2 == perfect_square
end

"""
    side_dimension_for_vectorized_dimension(n::Integer)

Return the dimension `d` such that
`MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(d))` is `n`.
"""
function side_dimension_for_vectorized_dimension(n::Base.Integer)
    return div(isqrt(1 + 8n), 2)
end


text:
#
# [MLJ.jl]: https://github.com/alan-turing-institute/MLJ.jl
# [RDatasets.jl]: https://github.com/JuliaStats/RDatasets.jl
# [DecisionTree.jl]: https://github.com/bensadeghi/DecisionTree.jl
#
# ## Preliminary steps
#
# ### Data
#
# As in "[choosing a model](choosing-a-model.html)", let's load the Iris dataset and unpack it:

using MLJ, Statistics
MLJ.color_off() #src
X, y = @load_iris;

# let's also load the `DecisionTreeClassifier`:

@load DecisionTreeClassifier

# ### MLJ Machine
#
# In MLJ, remember that a *model* is an object that only serves as a container for the hyperparameters of the model.
# A *machine* is an object wrapping both a model and data and can contain information on the *trained* model; it does *not* fit the model by itself.
# However, it does check that the model is compatible with the scientific type of the data and will warn you otherwise.

tree_model = DecisionTreeClassifier()
tree = machine(tree_model, X, y)

# A machine is used both for supervised and unsupervised model.
# In this tutorial we give an example for the supervised model first and then go on with the unsupervised case.
#
# ## Training and testing a supervised model
#
# Now that you've declared the model you'd like to consider and the data, we are left with the standard training and testing step for a supervised learning algorithm.
#
# ### Splitting the data
#
# To split the data into a *training* and *testing* set, you can use the function `partition` to obtain indices for data points that should be considered either as training or testing data:

train, test = partition(eachindex(y), 0.7, shuffle=true)
test[1:3]

# ### Fitting and testing the machine
#
# To fit the machine, you can use the function `fit!` specifying the rows to be used for the training:

fit!(tree, rows=train)

# Note that this **modifies** the machine which now contains the trained parameters of the decision tree.
# You can inspect the result of the fitting with the `fitresult` field of the model object:

tree.fitresult

# This `fitresult` will vary from model to model though classifiers will usually give out a tuple with the first element corresponding to the fitting and the second one keeping track of how classes are named (so that predictions can be appropriately named).
#
# You can now use the machine to make predictions with the `predict` function specifying rows to be used for the prediction:

ŷ = predict(tree, rows=test)
@show ŷ[1]

# Note that the output is *probabilistic*, effectively a vector with a score for each class.
# You could get the mode by using the `mode` function on `ŷ` or using `predict_mode`:

ȳ = predict_mode(tree, rows=test)
@show ȳ[1]
@show mode(ŷ[1])

# To measure the discrepancy between `ŷ` and `y` you could use the mean cross entropy:

mce = cross_entropy(ŷ, y[test]) |> mean
round(mce, digits=4)

# ## Unsupervised models
#
# Unsupervised models define a `transform` method,
# and may optionally implement an `inverse_transform` method.
# As in the supervised case, we use a machine to wrap the unsupervised model and the data:

v = [1, 2, 3, 4]
stand_model = UnivariateStandardizer()
stand = machine(stand_model, v)

# We can then fit the machine and use it to apply the corresponding *data transformation*:

fit!(stand)
w = transform(stand, v)
@show round.(w, digits=2)
@show mean(w)
@show std(w)

# In this case, the model also has an inverse transform:

vv = inverse_transform(stand, w)
sum(abs.(vv .- v))

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
struct Figure10
    lowcutoff::Vector{Float64}
    highcutoff::Float64
    cutdates::Vector{Float64}
    power::Matrix{Float64}
    Ω::Vector{Float64}
    specdates::Vector{Float64}
    function Figure10()
        @load "dat/anonomous_data.jld"
        @load "dat/wavefits.jld"
        wavetime_days = wavetime ./ 24 ./ 3600
        spectrogram, Ω, specdates = PaperAnalysis.compute_spectrogram(waveseries,wavetime_days)
        new([f.ini.cutoff for f in fits], 3.8, wavetime_days[fitstart:fitend], spectrogram, Ω, specdates)
    end
end

@recipe function f(f::Figure10)
    xlims := (f.specdates[1],f.specdates[end])
    ylims := (0,1.28π)
    xticks := (1:5,map(i->"day $i",1:5))
    @series begin
        seriestype := :heatmap
        f.specdates, f.Ω, 10 .* log10.(f.power)
    end
    @series begin
        linestyle := :dot
        seriescolor := :black
        label := false
        linewidth := 1.3
        f.cutdates, _-> f.highcutoff
    end
    @series begin
        linestyle := :dot
        seriescolor := :black
        label := false
        linewidth := 1.3
        f.cutdates, f.lowcutoff
    end
    @series begin
        linestyle := :dash
        seriescolor := :black
        label := false
        linewidth := 1.3
        seriestype := :vline
        f.cutdates[[1,end]]
    end
end


text:
using MIMEBundles
using Test
using JSON

@testset "MIMEBundle for String" begin
    bundle = MIMEBundle("foo")
    mimes = keys(bundle.data)
    @test length(mimes) == 1
    @test first(mimes) == MIME("text/plain")
    @test bundle.data[MIME("text/plain")] == "\"foo\""
end


text:
module UnionArraysTests

include("utils.jl")

function include_tests(m = @__MODULE__, dir = @__DIR__)
    for file in readdir(dir)
        if match(r"^test_.*\.jl$", file) !== nothing
            Base.include(m, joinpath(dir, file))
        end
    end
end

include_tests()

end  # module UnionArraysTests


text:
@testset "bearings.jl" begin
    p1 = rad2deg(Point(0.5, 1.0))
    p2 = rad2deg(Point(1.0, 0.5))
    p3 = rad2deg(Point(0.6, 1.0))  # above p1

    rs = RouteSection(p1,p2)
    @test bearing(p1, p2) ≈ rad2deg(5.81412807071) atol = 0.001
    @test bearing(rs) ≈ rad2deg(5.81412807071) atol = 0.001
    @test bearing(p3, p1) ≈ 180.0 atol = 0.001

    @test final_bearing(p1, p2) ≈ rad2deg(5.45864813531) atol = 0.001
    @test final_bearing(rs) ≈ rad2deg(5.45864813531) atol = 0.001
end


text:
#!/usr/bin/env julia

#
# Task worker
# Connects PULL socket to tcp://localhost:5557
# Collects workloads from ventilator via that socket
# Connects PUSH socket to tcp://localhost:5558
# Sends results to sink via that socket
#

using ZMQ

context = Context()

# Socket to receive messages on
receiver = Socket(context, PULL)
ZMQ.connect(receiver, "tcp://localhost:5557")

# Socket to send messages to
sender = Socket(context, PUSH)
ZMQ.connect(sender, "tcp://localhost:5558")

# Process tasks forever
while true
    s = bytestring(ZMQ.recv(receiver))

    # Simple progress indicator for the viewer
    write(STDOUT, ".")
    flush(STDOUT)

    # Do the work
    sleep(int(s)*0.001)

    # Send results to sink
    ZMQ.send(sender, b"")
end

max_stars_repo_name: JailbreakFox/LightWeightRepository

max_issues_repo_name: JailbreakFox/LightWeightRepository

max_forks_repo_name: JailbreakFox/LightWeightRepository


text:
## Parsing ##

function parse_meta(stream::IO)
    eventnum = Parsers.parse(Int64, stream)
    hitcount = Parsers.parse(Int64, stream)
    primarycount = Parsers.parse(Int64, stream)
    skip(stream, 1)
    return eventnum, hitcount, primarycount
end

function parse_hit(stream::IO)
    h = Hit(
        Parsers.parse(Float64, stream),
        Parsers.parse(Float64, stream),
        Parsers.parse(Float64, stream),
        Parsers.parse(Float64, stream),
        Parsers.parse(Float64, stream),
        Parsers.parse(Int64, stream),
        Parsers.parse(Int64, stream),
        Parsers.parse(Int64, stream)
    )
    skip(stream, 9)
    return h
end

## Reading events ##

struct RootHitIter <: AbstractHitIter
    hitcount::Integer
    stream::IO
end

Base.length(itr::RootHitIter) = itr.hitcount
Base.eltype(::Type{RootHitIter}) = Hit

function Base.iterate(itr::RootHitIter, i = 0)
    i == itr.hitcount && return nothing
    return parse_hit(itr.stream), i + 1
end

struct RootHitReader
    stream::IO
    ownstream::Bool
    lock::ReentrantLock
end

function RootHitReader(stream::IOStream, ownstream::Bool)
    RootHitReader(IOBuffer(mmap(stream)), ownstream, stream.lock)
end

function RootHitReader(stream::IO, ownstream::Bool)
    RootHitReader(stream, ownstream, ReentrantLock())
end

Base.IteratorSize(::Type{RootHitReader}) = Base.SizeUnknown()

Base.IteratorEltype(::Type{RootHitReader}) = Base.HasEltype()
Base.eltype(::Type{RootHitReader}) = Event{RootHitIter}

function Base.iterate(reader::RootHitReader, state = nothing)
    eof(reader) && return (close(reader); nothing)
    return read(reader, Event{Vector{Hit}}), nothing
end

function Base.read(reader::RootHitReader, T=Event{Vector{Hit}})
    lock(reader.lock) do
        enum, hitcnt, primcnt = parse_meta(reader.stream)
        return convert(T, Event(enum, hitcnt, primcnt, RootHitIter(hitcnt, reader.stream)))
    end
end

function Base.eof(reader::RootHitReader)
    lock(reader.lock) do
        eof(reader.stream)
    end
end

function Base.close(reader::RootHitReader)
    if reader.ownstream
        lock(reader.lock) do
            close(reader.stream)
        end
    end
end

loadstreaming(stream::IO) = RootHitReader(stream, false)
loadstreaming(f::Function, stream::IO) = f(loadstreaming(stream))

loadstreaming(path::AbstractString) = RootHitReader(open(path, lock=false), true)
function loadstreaming(f::Function, path::AbstractString)
    reader = loadstreaming(path)
    try
        return f(reader)
    finally
        close(reader)
    end
end

is_root_hit_file(path::AbstractString) = occursin(r"root.hits$", path)

function load(stream::IO; T=Event{Vector{Hit}})
    loadstreaming(stream) do stream
        return collect(T, stream)
    end
end

function load(path::AbstractString; T=Event{Vector{Hit}})
    is_root_hit_file(path) || error("cannot read events from $f")
    open(path, lock=false) do f
        return load(f; T=T)
    end
end


text:
include("utils.jl")
include("commands.jl")
include("inspect-vectors.jl")
using JLD2

function create_prototypes(trainfile, paramsfile)
    D = loadjson(trainfile; textkey="text", labelkey="klass")
    ncenters_per_class = 30
    kind=:fft
    kncconfig = KncProtoConfig(
        ncenters=-ncenters_per_class,
        maxiters=10,
        initial_clusters=kind,
        kernel=DirectKernel(CosineDistance()),
        centerselection=TextCentroidSelection()
    )
    best = load(paramsfile, "population") |> first |> first
    textconfig = best.textconfig # copy(best.textconfig, nlist=[3])
    config = MicroTC_Config(textconfig, best.textmodel, kncconfig)
    model = MicroTC(config, D.corpus, D.labels; tok=Tokenizer(config.textconfig))
    jldsave("H.model"; model)
    model
    # savemodel("H.model", model)
end

#model = load("H.model", "model")
#model = create_prototypes("pan21/train-es.json_train.json", "_pan21-es-microtc-64+acc+bsize=4/params")

max_stars_repo_name: sadit/text-classification-scripts

max_issues_repo_name: sadit/text-classification-scripts

max_forks_repo_name: sadit/text-classification-scripts


text:
function ziperate(v::I) where I
    maxlen = maximum(map(length,v))
    return zip(map(v) do x
               exterate(x,maxlen)
               end...)
end

max_stars_repo_name: sebastianpech/SpecialIterators

max_issues_repo_name: sebastianpech/SpecialIterators

max_forks_repo_name: sebastianpech/SpecialIterators


text:
# TODO: extend type restriction to `IndexT <: Union{<: Index, <: IndexVal}`
struct IndexSet{N, IndexT <: Index, DataT <: Tuple}
  data::DataT

  IndexSet{N,
           IndexT,
           DataT}(data) where {N,
                               IndexT,
                               DataT <: NTuple{N, IndexT}} =
    new{N, IndexT, DataT}(data)

  """
      IndexSet{Any}()

  Create a special "empty" IndexSet with data `Tuple{}` and `Any` number of indices.

  This is used as the IndexSet of an `emptyITensor()`, an ITensor with `NDTensors.Empty` storage and any number of indices.
  """
  IndexSet{Any}() = new{Any, Union{}, Tuple{}}()
end

#
# Order value type
#

@eval struct Order{N}
  (OrderT::Type{ <: Order})() = $(Expr(:new, :OrderT))
end

@doc """
    Order{N}

A value type representing the order of an ITensor.
""" Order

"""
    Order(N) = Order{N}()

Create an instance of the value type Order representing
the order of an ITensor.
"""
Order(N) = Order{N}()

"""
    IndexSet(::Function, ::Order{N})

Construct an IndexSet of length N from a function that accepts
an integer between 1:N and returns an Index.

# Examples
```julia
IndexSet(n -> Index(1, "i\$n"), Order(4))
```
"""
IndexSet(f::Function, N::Int) =
  IndexSet(ntuple(f, N))

IndexSet(f::Function, ::Order{N}) where {N} =
  IndexSet(ntuple(f, Val(N)))

# Definition to help with generic code
const Indices{N, IndexT} = Union{IndexSet{N,
                                          IndexT,
                                          NTuple{N, IndexT}},
                                 NTuple{N, IndexT}}

function (IndexSet{N, IndexT}(inds)::IndexSet{N, IndexT, NTuple{N, IndexT}}) where {N, IndexT}
  data = NTuple{N, IndexT}(inds)
  return IndexSet{N, IndexT, NTuple{N, IndexT}}(data)
end

"""
    IndexSet{N, IndexT}(inds)
    IndexSet{N, IndexT}(inds::Index...)

Construct an IndexSet of order N and element type IndexT
from a collection of indices (any collection that is convertable to a Tuple).
"""
function IndexSet{N, IndexT}(inds::Index...) where {N, IndexT}
  return IndexSet{N, IndexT}(inds)
end

"""
    IndexSet{N}(inds)
    IndexSet{N}(inds::Index...)

Construct an IndexSet of order `N` from a collection of indices
(any collection that is convertable to a Tuple).
"""
IndexSet{N}(inds) where {N} =
  IndexSet{N, eltype(inds)}(inds)

IndexSet{N}(inds::NTuple{N, IndexT}) where {N, IndexT} =
  IndexSet{N, IndexT, NTuple{N, IndexT}}(inds)

IndexSet{N}(inds::Index...) where {N} = IndexSet{N}(inds)

IndexSet{N}(is::IndexSet{N}) where {N} = is

IndexSet{0}(::Tuple{}) = IndexSet{0, Union{}, Tuple{}}(())

IndexSet{0}() = IndexSet{0}(())

"""
    IndexSet(inds)
    IndexSet(inds::Index...)

Construct an IndexSet from a collection of indices
(any collection that is convertable to a Tuple).
"""
IndexSet(inds) = IndexSet{length(inds)}(inds)

IndexSet(inds::NTuple{N, IndexT}) where {N, IndexT} = IndexSet{N, IndexT, NTuple{N, IndexT}}(inds)

IndexSet(::Tuple{}) = IndexSet{0, Union{}, Tuple{}}(())

IndexSet() = IndexSet(())

IndexSet(inds::Index...) = IndexSet(inds)

IndexSet(is::IndexSet) = is

"""
    convert(::Type{IndexSet}, t)

Convert the collection to an IndexSet,
as long as it can be converted to a Tuple.
"""
Base.convert(::Type{IndexSet}, t) = IndexSet(t)

Base.convert(::Type{IndexSet}, is::IndexSet) = is

Base.convert(::Type{IndexSet{N}}, t) where {N} = IndexSet{N}(t)

Base.convert(::Type{IndexSet{N}}, is::IndexSet{N}) where {N} = is

Base.convert(::Type{IndexSet{N,IndexT}}, t) where {N,IndexT} = IndexSet{N,IndexT}(t)

Base.convert(::Type{IndexSet{N,IndexT}}, is::IndexSet{N,IndexT}) where {N,IndexT} = is

Base.Tuple(is::IndexSet) = Tuple(data(is))

"""
    IndexSet(inds::Vector{<:Index})

Convert a Vector of indices to an IndexSet.

Warning: this is not type stable, since a Vector
is dynamically sized and an IndexSet is statically sized.
Consider using the constructor `IndexSet{N}(inds::Vector)`.
"""
IndexSet(inds::Vector{IndexT}) where {IndexT} = IndexSet(tuple(inds...))

"""
    IndexSet{N}(inds::Vector{<:Index})

Convert a Vector of indices to an IndexSet of size N.

Type stable conversion of a Vector of indices to an IndexSet
(in contrast to `IndexSet(::Vector{<:Index})`).
"""
IndexSet{N}(inds::Vector{IndexT}) where {N, IndexT} =
  IndexSet{N, IndexT, NTuple{N, IndexT}}(tuple(inds...))

"""
    not(inds::Union{IndexSet, Tuple{Vararg{<:Index}}})
    not(inds::Index...)
    !(inds::Union{IndexSet, Tuple{Vararg{<:Index}}})
    !(inds::Index...)

Represents the set of indices not in the specified
IndexSet, for use in pattern matching (i.e. when
searching for an index or priming/tagging specified
indices).
"""
not(is::IndexSet) = Not(is)
not(inds::Index...) = not(IndexSet(inds...))
not(inds::Tuple{Vararg{<:Index}}) = not(IndexSet(inds))
Base.:!(is::IndexSet) = not(is)
Base.:!(inds::Index...) = not(inds...)
Base.:!(inds::Tuple{Vararg{<:Index}}) = not(inds)

"""
    NDTensors.data(is::IndexSet)

Return the raw storage data for the indices.
Currently the storage is a Tuple.

This is mostly for internal usage, please
contact us if there is functionality you want
availabe for IndexSet.
"""
data(is::IndexSet) = is.data

# This is used in type promotion in the Tensor contraction code
Base.promote_rule(::Type{<:IndexSet},
                  ::Type{Val{N}}) where {N} = IndexSet{N}

NDTensors.ValLength(::Type{<:IndexSet{N}}) where {N} = Val{N}

NDTensors.ValLength(::IndexSet{N}) where {N} = Val(N)

Order(::IndexSet{N}) where {N} = Order(N)

Order(::Type{<:IndexSet{N}}) where {N} = Order(N)

# Convert to an Index if there is only one
# TODO: also define the function `only`
function Index(is::IndexSet)
  length(is) != 1 && error("Number of Index in IndexSet ≠ 1")
  return is[1]
end

"""
    getindex(is::IndexSet, n::Int)

Get the Index of the IndexSet in the dimension n.
"""
Base.getindex(is::IndexSet, n) = getindex(data(is), n)

"""
    getindex(is::IndexSet, v::AbstractVector)

Get the indices of the IndexSet in the dimensions
specified in the collection v, returning an IndexSet.

Warning: this function is not type stable.
"""
Base.getindex(is::IndexSet,
              v::AbstractVector) = IndexSet(getindex(data(is), v))

"""
    setindex(is::IndexSet, i::Index, n::Int)

Set the Index of the IndexSet in the dimension n.

This function is mostly for internal use, if you want to
replace the indices of an IndexSet, use the `replaceind`,
`replaceind!`, `replaceinds`, and `replaceinds!` functions,
which check for compatibility of the indices and ensure
the proper Arrow directions.
"""
function Base.setindex(is::IndexSet,
                       i::Index,
                       n::Int)
  return IndexSet(setindex(data(is), i, n))
end

"""
    length(is::IndexSet)

The number of indices in the IndexSet.
"""
Base.length(::IndexSet{N}) where {N} = N

"""
length(::Type{<:IndexSet})

The number of indices in the IndexSet type.
"""
Base.length(::Type{<:IndexSet{N}}) where {N} = N

"""
    size(is::IndexSet)

The size of the IndexSet, the same as `(length(is),)`.

Mostly for internal use for compatability with Base methods,
like for broadcasting.
"""
Base.size(is::IndexSet) = size(data(is))

"""
    axes(is::IndexSet)

The axes of the IndexSet, the same as `(Base.OneTo(length(is)),)`.

Mostly for internal use for compatability with Base methods,
like for broadcasting.
"""
Base.axes(is::IndexSet) = axes(data(is))

NDTensors.dims(is::IndexSet{N}) where {N} = dims(Tuple(is))

NDTensors.dims(is::NTuple{N,<:Index}) where {N} = ntuple(i->dim(is[i]),Val(N))

"""
    dim(is::IndexSet)

Get the product of the dimensions of the indices
of the IndexSet (the total dimension of the space).
"""
NDTensors.dim(is::IndexSet) = dim(Tuple(is))

NDTensors.dim(is::IndexSet{0}) = 1

NDTensors.dim(is::Tuple{Vararg{<:Index}}) = prod(dims(is))

"""
    dim(is::IndexSet, n::Int)

Get the dimension of the Index n of the IndexSet.
"""
NDTensors.dim(is::IndexSet, pos::Int) = dim(is[pos])

"""
    NDTensors.strides(is::IndexSet)

Get the strides of the n-dimensional tensor if it had
this IndexSet.
"""
NDTensors.strides(is::IndexSet) = Base.size_to_strides(1, dims(is)...)

"""
    NDTensors.stride(is::IndexSet. i::Int)

Get the stride of the n-dimensional tensor if it had
this IndexSet in the dimension `i`.
"""
NDTensors.stride(is::IndexSet, k::Int) = NDTensors.strides(is)[k]

"""
    dag(is::IndexSet)

Return a new IndexSet with the indices daggered (flip
all of the arrow directions).
"""
dag(is::IndexSet) = map(i -> dag(i), is)

"""
    iterate(is::IndexSet[, state])

Iterate over the indices of an IndexSet.

# Example
```jldoctest
julia> using ITensors;

julia> i = Index(2);

julia> is = IndexSet(i,i');

julia> for j in is
         println(plev(j))
       end
0
1
```
"""
Base.iterate(is::IndexSet, state) = iterate(data(is), state)

Base.iterate(is::IndexSet) = iterate(data(is))

# To allow for the syntax is[end]
Base.firstindex(::IndexSet) = 1

# To allow for the syntax is[begin]
Base.lastindex(is::IndexSet) = length(is)

"""
    eltype(::IndexSet)

Get the element type of the IndexSet.
"""
Base.eltype(is::IndexSet{<: Any, IndexT}) where {IndexT} = IndexT

Base.eltype(::Type{<: IndexSet{<: Any,
                               IndexT}}) where {IndexT} = IndexT

# Needed for findfirst (I think)
Base.keys(is::IndexSet{N}) where {N} = 1:N

# This is to help with some generic programming in the Tensor
# code (it helps to construct an IndexSet(::NTuple{N,Index}) where the
# only known thing for dispatch is a concrete type such
# as IndexSet{4})
NDTensors.similar_type(::Type{<:IndexSet},
                       ::Val{N}) where {N} = IndexSet{N}

NDTensors.similar_type(::Type{<:IndexSet},
                       ::Type{Val{N}}) where {N} = IndexSet{N}

"""
    sim(is::IndexSet)

Make a new IndexSet with similar indices.

You can also use the broadcast version `sim.(is)`.
"""
sim(is::IndexSet) = map(i -> sim(i), is)

# For generic code in NDTensors
NDTensors.sim(is::IndexSet) = sim(is)

import .NDTensors: mindim
"""
    mindim(is::IndexSet)

Get the minimum dimension of the indices in the index set.

Returns 1 if the IndexSet is empty.
"""
function mindim(is::IndexSet)
  length(is) == 0 && (return 1)
  md = dim(is[1])
  for n in 2:length(is)
    md = min(md, dim(is[n]))
  end
  return md
end

"""
    maxdim(is::IndexSet)

Get the maximum dimension of the indices in the index set.

Returns 1 if the IndexSet is empty.
"""
function maxdim(is::IndexSet)
  length(is) == 0 && (return 1)
  md = dim(is[1])
  for n ∈ 2:length(is)
    md = max(md,dim(is[n]))
  end
  return md
end

#
# Set operations
#

"""
    ==(is1::IndexSet, is2::IndexSet)

IndexSet equality (order dependent). For order
independent equality use `issetequal` or
`hassameinds`.
"""
function Base.:(==)(A::IndexSet, B::IndexSet)
  length(A) ≠ length(B) && return false
  for (a,b) in zip(A,B)
    a ≠ b && return false
  end
  return true
end

"""
    ITensors.fmatch(pattern) -> ::Function

fmatch is an internal function that
creates a function that accepts an Index.
The function returns true if the Index matches
the provided pattern, and false otherwise.

For example:
```
i = Index(2, "s")
fmatch("s")(i) == true
```
"""
fmatch(is::IndexSet) = in(is)
fmatch(is::Tuple{Vararg{<:Index}}) = fmatch(IndexSet(is))
fmatch(is::Index...) = fmatch(IndexSet(is...))

fmatch(pl::Int) = hasplev(pl)

fmatch(tags::TagSet) = hastags(tags)
fmatch(tags::AbstractString) = fmatch(TagSet(tags))

fmatch(id::IDType) = hasid(id)

fmatch(n::Not) = !fmatch(parent(n))

# Function that always returns true
ftrue(::Any) = true

fmatch(::Nothing) = ftrue

"""
    ITensors.fmatch(; inds = nothing,
                      tags = nothing,
                      plev = nothing,
                      id = nothing) -> Function

An internal function that returns a function
that accepts an Index that checks if the
Index matches the provided conditions.
"""
function fmatch(; inds = nothing,
                  tags = nothing,
                  plev = nothing,
                  id = nothing)
  return i -> fmatch(inds)(i) &&
              fmatch(plev)(i) &&
              fmatch(id)(i) &&
              fmatch(tags)(i)
end

"""
    indmatch

Checks if the Index matches the provided conditions.
"""
indmatch(i::Index; kwargs...) = fmatch(; kwargs...)(i)

# Set functions
#
# setdiff
# intersect
# symdiff
# union
# filter
#

function Base.setdiff!(f::Function,
                       r,
                       A::IndexSet,
                       Bs::IndexSet...)

  N = length(r)
  j = 1
  for a in A
    if f(a) && all(B -> a ∉ B, Bs)
      j > N && error("Too many intersects found")
      r[j] = a
      j += 1
    end
  end
  j ≤ N && error("Too few intersects found")
  return r
end

function Base.setdiff(f::Function,
                      A::IndexSet,
                      Bs::IndexSet...)
  R = eltype(A)[]
  for a ∈ A
    f(a) && all(B -> a ∉ B, Bs) && push!(R, a)
  end
  return R
end

#Base.setdiff(f::Function, is1::IndexSet, iss::IndexSet...) =
#  setdiff(Order(count(i -> f(i) && all(is -> i ∉ is, iss), is1)), f, is1, iss...)

"""
    setdiff(A::IndexSet, Bs::IndexSet...)

Output the IndexSet with Indices in `A` but not in
the IndexSets `Bs`.
"""
Base.setdiff(A::IndexSet, Bs::IndexSet...; kwargs...) =
  setdiff(fmatch(; kwargs...), A, Bs...)

Base.setdiff(O::Order, A::IndexSet, Bs::IndexSet...;
             kwargs...) =
  setdiff(fmatch(; kwargs...), O, A, Bs...)

"""
  setdiff(f::Function, ::Order{N}, A::IndexSet, B::IndexSet...)

Output the IndexSet in the set difference of `A` and `B`,
optionally filtering by the function `f`.

Specify the number of output indices with `Order(N)`.
"""
function Base.setdiff(f::Function,
                      ::Order{N},
                      A::IndexSet{<:Any, IndexT},
                      B::IndexSet...) where {N, IndexT}
  r = mutable_storage(Order{N}, IndexT)
  setdiff!(f, r, A, B...)
  return IndexSet{N, IndexT, NTuple{N, IndexT}}(Tuple(r))
end

function firstsetdiff(f::Function,
                      A::IndexSet,
                      Bs::IndexSet...)
  for a in A
    f(a) && all(B -> a ∉ B, Bs) && return a
  end
  return nothing
end

"""
    firstsetdiff(A::IndexSet, Bs::IndexSet...)

Output the first Index in `A` that is not in the IndexSets `Bs`.
Otherwise, return a default constructed Index.
"""
firstsetdiff(A::IndexSet,
             Bs::IndexSet...;
             kwargs...) = firstsetdiff(fmatch(; kwargs...), A, Bs...)

function Base.intersect(f::Function, A::IndexSet, B::IndexSet)
  R = eltype(A)[]
  for a in A
    f(a) && a ∈ B && push!(R,a)
  end
  return R
end

#Base.intersect(f::Function, is1::IndexSet, is2::IndexSet) =
#  intersect(f, Order(count(i -> f(i) && i ∈ is2, is1)), is1, is2)

mutable_storage(::Type{Order{N}},
                ::Type{IndexT}) where {N, IndexT <: Index} =
  MVector{N, IndexT}(undef)

function Base.intersect(f::Function,
                        ::Order{N},
                        A::IndexSet{<:Any, IndexT},
                        B::IndexSet) where {N, IndexT}
  R = mutable_storage(Order{N}, IndexT)
  intersect!(f, R, A, B)
  return IndexSet{N, IndexT, NTuple{N, IndexT}}(Tuple(R))
end

Base.intersect(O::Order,
               A::IndexSet,
               B::IndexSet;
               kwargs...) where {N} =
  intersect(fmatch(; kwargs...), O, A, B)

function Base.intersect!(f::Function,
                         R::AbstractVector,
                         A::IndexSet,
                         B::IndexSet)
  N = length(R)
  j = 1
  for a in A
    if f(a) && a ∈ B
      j > N && error("Too many intersects found")
      R[j] = a
      j += 1
    end
  end
  j ≤ N && error("Too few intersects found")
  return R
end

"""
    intersect(A::IndexSet, B::IndexSet; kwargs...)

    intersect(f::Function, A::IndexSet, B::IndexSet)

    intersect(:Order{N}, A::IndexSet, B::IndexSet; kwargs...)

    intersect(f::Function, ::Order{N}, A::IndexSet, B::IndexSet)

Output the IndexSet in the intersection of `A` and `B`,
optionally filtering with keyword arguments `tags`, `plev`, etc.
or by a function `f(::Index) -> Bool`.

Specify the output number of indices `N` with `Order(N)`.
"""
Base.intersect(A::IndexSet, B::IndexSet; kwargs...) =
  intersect(fmatch(; kwargs...), A, B)

function firstintersect(f::Function, A::IndexSet, B::IndexSet)
  for a in A
    f(a) && a ∈ B && return a
  end
  return nothing
end

"""
    firstintersect(A::IndexSet, B::IndexSet; kwargs...)

    firstintersect(f::Function, A::IndexSet, B::IndexSet)

Output the first Index common to `A` and `B`, optionally
filtering by tags, prime level, etc. or by a function
`f`.

If no common Index is found, return `nothing`.
"""
firstintersect(A::IndexSet, B::IndexSet; kwargs...) =
  firstintersect(fmatch(; kwargs...), A, B)

"""
    filter(f::Function, inds::IndexSet)

    filter(f::Function, ::Order{N}, inds::IndexSet)

Filter the IndexSet by the given function (output a new
IndexSet with indices `i` for which `f(i)` returns true).

Note that this function is not type stable, since the number
of output indices is not known at compile time.

To make it type stable, specify the desired order by
passing an instance of the type `Order`.
"""
Base.filter(f::Function,
            is::IndexSet) =
  IndexSet(filter(f, Tuple(is)))

Base.filter(is::IndexSet, args...; kwargs...) =
  filter(fmatch(args...; kwargs...), is)

# To fix ambiguity error with Base function
Base.filter(is::IndexSet, tags::String; kwargs...) =
  filter(fmatch(tags; kwargs...),is)

function Base.filter!(f::Function,
                      r,
                      is::IndexSet{<:Any, IndexT}) where {IndexT}
  N = length(r)
  j = 1
  for i in is
    if f(i)
      j > N && error("Too many intersects found")
      r[j] = i
      j += 1
    end
  end
  j ≤ N && error("Too few intersects found")
  return r
end

function Base.filter(f::Function,
                     O::Order{N},
                     is::IndexSet{<:Any, IndexT}) where {N, IndexT}
  #t = filter(f, Tuple(is))
  #return IndexSet{N, IndexT, NTuple{N, IndexT}}(t)
  r = mutable_storage(Order{N}, IndexT)
  filter!(f, r, is)
  return IndexSet{N, IndexT, NTuple{N, IndexT}}(Tuple(r))
end

Base.filter(O::Order, is::IndexSet, args...; kwargs...) =
  filter(ITensors.fmatch(args...; kwargs...), O, is)

"""
    getfirst(is::IndexSet)

Return the first Index in the IndexSet. If the IndexSet
is empty, return `nothing`.
"""
function getfirst(is::IndexSet)
  length(is) == 0 && return nothing
  return first(is)
end

"""
    getfirst(f::Function, is::IndexSet)

Get the first Index matching the pattern function,
return `nothing` if not found.
"""
function getfirst(f::Function, is::IndexSet)
  for i in is
    f(i) && return i
  end
  return nothing
end

getfirst(is::IndexSet,
         args...; kwargs...) = getfirst(fmatch(args...;
                                               kwargs...),is)

Base.findall(is::IndexSet,
             args...; kwargs...) = findall(fmatch(args...;
                                                  kwargs...), is)

"""
    indexin(ais::IndexSet, bis::IndexSet)

For collections of Indices, returns the first location in
`bis` for each value in `ais`, as a Tuple.
"""
function Base.indexin(ais::IndexSet,
                      bis::IndexSet)
  return ntuple(i -> findfirst(bis, ais[i]), Val(length(ais)))
end

Base.findfirst(is::IndexSet,
               args...; kwargs...) = findfirst(fmatch(args...;
                                                      kwargs...), is)

"""
    map(f, is::IndexSet)

Apply the function to the elements of the IndexSet,
returning a new IndexSet.
"""
Base.map(f::Function, is::IndexSet) = IndexSet(map(f, data(is)))

#
# Tagging functions
#

function prime(f::Function,
               is::IndexSet,
               args...)
  return map(i -> f(i) ? prime(i,args...) : i, is)
end

"""
    prime(A::IndexSet, plinc, ...)

Increase the prime level of the indices by the specified amount.
Filter which indices are primed using keyword arguments
tags, plev and id.
"""
prime(is::IndexSet,
      plinc::Integer,
      args...; kwargs...) = prime(fmatch(args...; kwargs...),
                                  is, plinc)

prime(f::Function,
      is::IndexSet) = prime(f, is, 1)

prime(is::IndexSet,
      args...;
      kwargs...) = prime(is, 1, args...; kwargs...)

"""
    adjoint(is::IndexSet)

For is' notation.
"""
Base.adjoint(is::IndexSet) = prime(is)

function setprime(f::Function,
                  is::IndexSet,
                  args...)
  return map(i -> f(i) ? setprime(i, args...) : i, is)
end

setprime(is::IndexSet,
         plev::Integer,
         args...; kwargs...) = setprime(fmatch(args...; kwargs...),
                                        is, plev)

noprime(f::Function,
        is::IndexSet,
        args...) = setprime(is, 0, args...; kwargs...)

noprime(is::IndexSet,
        args...;
        kwargs...) = setprime(is, 0, args...; kwargs...)

function _swapprime(f::Function,
                    i::Index,
                    pl1::Int,
                    pl2::Int)
  if f(i)
    if hasplev(i, pl1)
      return setprime(i, pl2)
    elseif hasplev(i, pl2)
      return setprime(i, pl1)
    end
    return i
  end
  return i
end

function swapprime(f::Function,
                   is::IndexSet,
                   pl1::Int,
                   pl2::Int)
  return map(i -> _swapprime(f, i, pl1, pl2), is)
end

swapprime(is::IndexSet,
          pl1::Int,
          pl2::Int,
          args...; kwargs...) = swapprime(fmatch(args...; kwargs...),
                                          is, pl1, pl2)

mapprime(f::Function, is::IndexSet,
         pl1::Int, pl2::Int) = mapprime(f, is, pl1 => pl2)

mapprime(is::IndexSet, pl1::Int, pl2::Int,
         args...; kwargs...) = mapprime(fmatch(args...; kwargs...),
                                        is, pl1 => pl2)

function _mapprime(i::Index,
                   rep_pls::Pair{Int, Int}...)
  for (pl1, pl2) in rep_pls
    hasplev(i, pl1) && return setprime(i, pl2)
  end
  return i
end

function mapprime(f::Function,
                  is::IndexSet,
                  rep_pls::Pair{Int, Int}...)
  return map(i -> f(i) ? _mapprime(i, rep_pls...) : i, is)
end

mapprime(is::IndexSet,
         rep_pls::Pair{Int, Int}...;
         kwargs...) = mapprime(fmatch(; kwargs...),
                               is, rep_pls...)

function addtags(f::Function,
                 is::IndexSet,
                 args...)
  return map(i -> f(i) ? addtags(i, args...) : i, is)
end

addtags(is::IndexSet,
        tags,
        args...; kwargs...) = addtags(fmatch(args...; kwargs...),
                                      is, tags)

function settags(f::Function,
                 is::IndexSet,
                 args...)
  return map(i -> f(i) ? settags(i, args...) : i, is)
end

settags(is::IndexSet,
        tags,
        args...;
        kwargs...) = settags(fmatch(args...; kwargs...),
                             is, tags)

function removetags(f::Function,
                    is::IndexSet,
                    args...)
  return map(i -> f(i) ? removetags(i, args...) : i, is)
end

removetags(is::IndexSet,
           tags,
           args...;
           kwargs...) = removetags(fmatch(args...; kwargs...),
                                   is, tags)

function _replacetags(i::Index,
                      rep_ts::Pair...)
  for (tags1, tags2) in rep_ts
    hastags(i, tags1) && return replacetags(i, tags1, tags2)
  end
  return i
end

# Version taking a list of Pairs
replacetags(f::Function,
            is::IndexSet,
            rep_ts::Pair...) =
  map(i -> f(i) ? _replacetags(i, rep_ts...) : i, is)

replacetags(is::IndexSet,
            rep_ts::Pair...;
            kwargs...) = replacetags(fmatch(; kwargs...),
                                     is, rep_ts...)

# Version taking two input TagSets/Strings
replacetags(f::Function,
            is::IndexSet,
            tags1,
            tags2) = replacetags(f, is, tags1 => tags2)

replacetags(is::IndexSet,
            tags1,
            tags2,
            args...;
            kwargs...) = replacetags(fmatch(args...; kwargs...),
                                     is, tags1 => tags2)

function _swaptags(f::Function,
                   i::Index,
                   tags1,
                   tags2)
  if f(i)
    if hastags(i, tags1)
      return replacetags(i, tags1, tags2)
    elseif hasteags(i, tags2)
      return replacetags(i, tags2, tags1)
    end
    return i
  end
  return i
end

function swaptags(f::Function,
                  is::IndexSet,
                  tags1,
                  tags2)
  return map(i -> _swaptags(f, i, tags1, tags2), is)
end

swaptags(is::IndexSet,
         tags1,
         tags2,
         args...;
         kwargs...) = swaptags(fmatch(args...; kwargs...),
                               is, tags1, tags2)

replaceinds(is::IndexSet, rep_inds::Pair{ <: Index, <: Index}...) =
  replaceinds(is, zip(rep_inds...)...)

function replaceinds(is::IndexSet, inds1, inds2)
  is1 = IndexSet(inds1)
  poss = indexin(is1, is)
  for (j, pos) in enumerate(poss)
    isnothing(pos) && continue
    i1 = is[pos]
    i2 = inds2[j]
    space(i1) != space(i2) && error("Indices must have the same spaces to be replaced")
    i2 = setdir(i2, dir(i1))
    is = setindex(is, i2, pos)
  end
  return is
end

replaceind(is::IndexSet, i1::Index, i2::Index) = replaceinds(is, (i1,), (i2,))

replaceind(is::IndexSet, rep_i::Pair{ <: Index, <: Index}) = replaceinds(is, rep_i)

function swapinds(is::IndexSet, inds1, inds2)
  return replaceinds(is, (inds1..., inds2...), (inds2..., inds1...))
end

swapind(is::IndexSet, i1::Index, i2::Index) = swapinds(is, (i1,), (i2,))

removeqns(is::IndexSet) = is

function permute(is1::IndexSet{N}, is2::IndexSet{N}) where {N}
  perm = NDTensors.getperm(is1, is2)
  return NDTensors.permute(is1, perm)
end

#
# Helper functions for contracting ITensors
#

function compute_contraction_labels(Ais::IndexSet{NA},
                                    Bis::IndexSet{NB}) where {NA,NB}
  Alabels = MVector{NA,Int}(ntuple(_->0,Val(NA)))
  Blabels = MVector{NB,Int}(ntuple(_->0,Val(NB)))

  ncont = 0
  for i = 1:NA, j = 1:NB
    if Ais[i] == Bis[j]
      Alabels[i] = Blabels[j] = -(1+ncont)
      ncont += 1
    end
  end

  u = ncont
  for i = 1:NA
    if(Alabels[i]==0) Alabels[i] = (u+=1) end
  end
  for j = 1:NB
    if(Blabels[j]==0) Blabels[j] = (u+=1) end
  end

  return (Tuple(Alabels),Tuple(Blabels))
end

function compute_contraction_labels(Cis::IndexSet{NC},
                                    Ais::IndexSet{NA},
                                    Bis::IndexSet{NB}) where {NC,NA,NB}
  Alabels,Blabels = compute_contraction_labels(Ais, Bis)
  Clabels = MVector{NC,Int}(ntuple(_->0,Val(NC)))
  for i = 1:NC
    locA = findfirst(==(Cis[i]), Ais)
    if !isnothing(locA)
      Clabels[i] = Alabels[locA]
    else
      locB = findfirst(==(Cis[i]), Bis)
      Clabels[i] = Blabels[locB]
    end
  end
  return (Tuple(Clabels),Alabels,Blabels)
end

#
# TupleTools
#

"""
    pop(is::IndexSet)

Return a new IndexSet with the last Index removed.
"""
pop(is::IndexSet) = IndexSet(NDTensors.pop(Tuple(is)))

# Overload the unexported NDTensors version
NDTensors.pop(is::IndexSet) = pop(is)

"""
    popfirst(is::IndexSet)

Return a new IndexSet with the first Index removed.
"""
popfirst(is::IndexSet) = IndexSet(NDTensors.popfirst(Tuple(is)))

# Overload the unexported NDTensors version
NDTensors.popfirst(is::IndexSet) = popfirst(is)

"""
    push(is::IndexSet, i::Index)

Make a new IndexSet with the Index i inserted
at the end.
"""
push(is::IndexSet,
     i::Index) = IndexSet(NDTensors.push(data(is), i))

# Overload the unexported NDTensors version
NDTensors.push(is::IndexSet,
             i::Index) = push(is, i)

push(is::IndexSet{0},
     i::Index) = IndexSet(i)

# Overload the unexported NDTensors version
NDTensors.push(is::IndexSet{0},
             i::Index) = push(is, i)

"""
    pushfirst(is::IndexSet, i::Index)

Make a new IndexSet with the Index i inserted
at the beginning.
"""
pushfirst(is::IndexSet,
          i::Index) = IndexSet(NDTensors.pushfirst(data(is), i))

# Overload the unexported NDTensors version
NDTensors.pushfirst(is::IndexSet,
                  i::Index) = pushfirst(is, i)

pushfirst(is::IndexSet{0},
          i::Index) = IndexSet(i)

# Overload the unexported NDTensors version
NDTensors.pushfirst(is::IndexSet{0},
                  i::Index) = pushfirst(is, i)

"""
    instertat(is1::IndexSet, is2, pos::Int)

Remove the index at pos and insert the indices
is2 starting at that position.
"""
function insertat(is1::IndexSet,
                  is2,
                  pos::Int)
  return IndexSet(NDTensors.insertat(Tuple(is1),
                                   Tuple(IndexSet(is2)),
                                   pos))
end

# Overload the unexported NDTensors version
NDTensors.insertat(is1::IndexSet,
                   is2,
                   pos::Int) = insertat(is1, is2, pos)

"""
    instertafter(is1::IndexSet, is2, pos)

Insert the indices is2 after position pos.
"""
function insertafter(is::IndexSet, I...)
  return IndexSet(NDTensors.insertafter(Tuple(is), I...))
end

# Overload the unexported NDTensors version
NDTensors.insertafter(is::IndexSet,
                    I...) = insertafter(is, I...)

function deleteat(is::IndexSet, I...)
  return IndexSet(NDTensors.deleteat(Tuple(is),I...))
end

# Overload the unexported NDTensors version
NDTensors.deleteat(is::IndexSet,
                   I...) = deleteat(is, I...)

function getindices(is::IndexSet, I...)
  return IndexSet(NDTensors.getindices(Tuple(is), I...))
end

NDTensors.getindices(is::IndexSet,
                     I...) = getindices(is, I...)

#
# QN functions
#

"""
    setdirs(is::IndexSet, dirs::Arrow...)

Return a new IndexSet with indices `setdir(is[i], dirs[i])`.
"""
function setdirs(is::IndexSet{N}, dirs) where {N}
  return IndexSet(ntuple(i -> setdir(is[i], dirs[i]), Val(N)))
end

"""
    dir(is::IndexSet, i::Index)

Return the direction of the Index `i` in the IndexSet `is`.
"""
function dir(is1::IndexSet, i::Index)
  return dir(getfirst(is1, i))
end

"""
    dirs(is::IndexSet, inds)

Return a tuple of the directions of the indices `inds` in
the IndexSet `is`, in the order they are found in `inds`.
"""
function dirs(is1::IndexSet, inds)
  return ntuple(i -> dir(is1, inds[i]), Val(length(inds)))
end

"""
    dirs(is::IndexSet)

Return a tuple of the directions of the indices `is`.
"""
function dirs(is::IndexSet{N}) where {N}
  return ntuple(i -> dir(is[i]), Val(N))
end

hasqns(is::IndexSet) = any(hasqns,is)

"""
    nblocks(::IndexSet, i::Int)

The number of blocks in the specified dimension.
"""
function NDTensors.nblocks(inds::IndexSet, i::Int)
  return nblocks(Tuple(inds),i)
end

function NDTensors.nblocks(inds::IndexSet, is)
  return nblocks(Tuple(inds),is)
end

"""
    nblocks(::IndexSet)

A tuple of the number of blocks in each
dimension.
"""
function NDTensors.nblocks(inds::IndexSet{N}) where {N}
  return ntuple(i->nblocks(inds,i),Val(N))
end

function NDTensors.nblocks(inds::NTuple{N,<:Index}) where {N}
  return nblocks(IndexSet(inds))
end

ndiagblocks(inds) = minimum(nblocks(inds))

# TODO: generic to IndexSet and BlockDims
function eachblock(inds::IndexSet)
  return CartesianIndices(nblocks(inds))
end

# TODO: turn this into an iterator instead
# of returning a Vector
function eachdiagblock(inds::IndexSet{N}) where {N}
  return [ntuple(_->i,Val(N)) for i in 1:ndiagblocks(inds)]
end

"""
    flux(inds::IndexSet, block::Tuple{Vararg{Int}})

Get the flux of the specified block, for example:
```
i = Index(QN(0)=>2, QN(1)=>2)
is = IndexSet(i, dag(i'))
flux(is, (1,1)) == QN(0)
flux(is, (2,1)) == QN(1)
flux(is, (1,2)) == QN(-1)
flux(is, (2,2)) == QN(0)
```
"""
function flux(inds::IndexSet, block)
  qntot = QN()
  for n in 1:length(inds)
    ind = inds[n]
    qntot += qn(ind, block[n])
  end
  return qntot
end

"""
    flux(inds::IndexSet, I::Int...)

Get the flux of the block that the specified
index falls in.
```
i = Index(QN(0)=>2, QN(1)=>2)
is = IndexSet(i, dag(i'))
flux(is, 3, 1) == QN(1)
flux(is, 1, 2) == QN(0)
```
"""
flux(inds::IndexSet,
     vals::Int...) = flux(inds, block(inds, vals...))

"""
    ITensors.block(inds::IndexSet, I::Int...)

Get the block that the specified index falls in.

This is mostly an internal function, and the interface
is subject to change.
```
i = Index(QN(0)=>2, QN(1)=>2)
is = IndexSet(i, dag(i'))
ITensors.block(is, 3, 1) == (2,1)
ITensors.block(is, 1, 2) == (1,1)
```
"""
block(inds::IndexSet,
      vals::Int...) = blockindex(inds, vals...)[2]

function Base.show(io::IO, is::IndexSet)
  print(io,"IndexSet{$(length(is))} ")
  for i in is
    print(io, i)
    print(io, " ")
  end
end

#
# Read and write
#

function readcpp(io::IO,::Type{<:IndexSet};kwargs...)
  format = get(kwargs,:format,"v3")
  is = IndexSet()
  if format=="v3"
    size = read(io,Int)
    function readind(io,n)
      i = readcpp(io,Index;kwargs...)
      stride = read(io,UInt64)
      return i
    end
    is = IndexSet(ntuple(n->readind(io,n),size))
  else
    throw(ArgumentError("read IndexSet: format=$format not supported"))
  end
  return is
end

function HDF5.write(parent::Union{HDF5File,HDF5Group},
                    name::AbstractString,
                    is::IndexSet)
  g = g_create(parent,name)
  attrs(g)["type"] = "IndexSet"
  attrs(g)["version"] = 1
  N = length(is)
  write(g,"length",N)
  for n=1:N
    write(g,"index_$n",is[n])
  end
end

function HDF5.read(parent::Union{HDF5File,HDF5Group},
                   name::AbstractString,
                   ::Type{<:IndexSet})
  g = g_open(parent,name)
  if read(attrs(g)["type"]) != "IndexSet"
    error("HDF5 group or file does not contain IndexSet data")
  end
  N = read(g,"length")
  it = ntuple(n->read(g,"index_$n",Index),N)
  return IndexSet(it)
end

import .NDTensors.store
@deprecate store(is::IndexSet) data(is)