text:
using Cliffords, LsqFit

function clean(m::Void)
    return ""
end

function clean(m::SubString)
    return m
end

"""
Parses a GST experiment sequence. The first argument is a string describing the
    experiment sequence, while the second argument is a dictionary mapping gate
    labels to integers corresponding to indices in a pulse library. The output is
    an array of integers corresponding to indices in a pulse library.

    Experiment sequences have the following format. Gates are denoted by string
    matching the regular expression "G[^G]+", e.g.,

       "G1", "Ga", or "G\$"

    A gate sequence matches a sequence of gates as defined above,  e.g.,

       "G1GbG#GaG%Gg"

    An experiment sequence matches a sequence of gates, potentially followed by
    a sequence power, i.e., a gate sequence in parenthesis followed by a carret
    and a non-negative integer (indicating repetition), followed
    by at most another gate sequence, e.g.,

       "(GaG1G\$)^4"

       "G1Gb(GcG]\$)^4G5G10"

       "G1Gb(GcG\$)^4"

       "(GcG\$)^4G5G10"

       "G1GbG5G10"

"""
function parse_gst(str, dict)
    function parse_count(str)
        count = r"\^([0-9]*).*"
        if length(str)>1 && str[1] == '^'
            parse(Int,match(count,str).captures[1])
        else
            1
        end
    end

    if contains("{}",str)
        return [dict["{}"]]
    else
        # split on parenthesis
        parts = split(str,['(',')'])
        # check if there is repetition (in the form of "^N")
        counts = map(parse_count,parts)
        counts = vcat(counts[2:end],[counts[1]])
        # extract subsequence strings
        clean_parts = [ match(r"(\^[0-9]*|)((G[^G]+)*)",part).captures[2] for part in parts ]
        # map subsequences to arrays of gate labels
        clean_arr = map(s->split(s,'G'),clean_parts)
        # discard empty sequences
        nontrivial_arr = Vector{Int}[]
        nontrivial_counts = Int[]
        for (index,item) in enumerate(clean_arr)
            if length(item) > 1
                push!(nontrivial_arr, map(s->dict[s],item[2:end]))
                push!(nontrivial_counts, counts[index])
            end
        end
        return vcat([repmat(nontrivial_arr[i], nontrivial_counts[i]) for i in 1:length(nontrivial_arr)]...)
    end
end

begin
    global multiply

    table = Dict{Tuple{Int,Int},Int}()

    function multiply(c1,c2)
        c1_ = round(Int,c1==0 ? 1 : c1)
        c2_ = round(Int,c1==0 ? 1 : c2)
        if !((c1_,c2_) in keys(table))
            table[(c1_,c2_)] = Cliffords.localcliffordindex(localclifford(c1_)*localclifford(c2_))
        end
        return table[(c1_,c2_)]
    end
end

"""
FitResult type that stores the output of a fit.
"""
struct FitResult
    fit_params::Dict{String,Float64}
    sq_error::Float64
    Nσ::Float64
    errors::Dict{String,Float64}
    fit_curve::Function
    model_str::String
end

"""
`generic_fit(xpts, ypts, model, initial_guess, fit_params, model_string; yvars=[])`
Helper function to return a fit to a model.
# Arguments
* xpts: Independent fit variable.
* ypts: Dependent fit variable.
* model: Fit function(x, p) of data x and parameter vector p.
* initial_guess: Initial guess of fit parameters.
* fit_params::function(Array{Float64} => Dict{String, Float64}): Result dict of fit parameters.
* model_string::String: String indicating fit function.
* yvars: Variance of y points. Defaults to empty vector.
# Returns
* A FitResult object for the specified fit.
"""
function generic_fit(xpts, ypts, model, initial_guess, fit_params, model_string::String; yvars=[])
    @assert length(xpts) == length(ypts) "X and Y data length must match."
    #try
        if isempty(yvars)
          result = curve_fit(model, xpts, ypts, initial_guess)
          errors = estimate_errors(result)
          sq_error = sum(((model(xpts, result.param) - ypts) .^ 2))
        else
          @assert length(ypts) == length(yvars) "Y data and Y variance lengths must match."
          result = curve_fit(model, xpts, ypts, 1 ./ sqrt(yvars), initial_guess)
          errors = estimate_errors(result)
          sq_error = sum(((model(xpts, result.param) - ypts) .^ 2) ./ yvars)
        end
    #catch
    #    result =
    # Compute badness of fit:
    # Under the null hypothesis (that the model is valid and that the observations
    # do indeed have Gaussian statistics), the mean squared error is χ² distributed
    # with `dof` degrees of freedom. We can quantify badness-of-fit in terms of how
    # far the observed MSE is from the expected value, in units of σ = 2dof (the expected
    # standard deviation for the χ² distribution)
    dof = length(xpts)-length(initial_guess)
    Nσ = sq_error/sqrt(2*dof) - dof/sqrt(2*dof)
    return FitResult( fit_params(result.param),
              sq_error,
              Nσ,
              fit_params(errors),
              xpts->model(xpts,result.param),
              model_string)
end

"""`fit_RB(xpts, ypts, yvars=[])`
Fit to RB decay a*(1-b)^x + c
"""
function fit_RB(xpts, ypts; yvars=[])
    RB_fit_dict(p) = Dict("a" => p[1], "b" => p[2], "c" => p[3])
    model(n, p) = p[1] * (1-p[2]).^n + p[3]
    p_guess = [0.5, .01, 0.5]
    return generic_fit(xpts, ypts, model, p_guess, RB_fit_dict, "a * exp(-n/b) + c", yvars=yvars)
end

max_stars_repo_name: matthewware/RGST-analysis

max_issues_repo_name: matthewware/RGST-analysis

max_forks_repo_name: matthewware/RGST-analysis


text:
using SparseArrays

struct Config
    centering::Bool

    function Config(;centering = true)
        new(centering)
    end
end

const defaultconfig = Config()

struct Memory
    nu::Int # number of users
    ni::Int # number of items
    bu::Vector{Float64} # user biases
    bi::Vector{Float64} # item biases
    Dui::SparseMatrixCSC{Float64,Int} # user-item residuals
    Diu::SparseMatrixCSC{Float64,Int} # item-user residuals
    Suu::SparseMatrixCSC{Float64,Int} # user-user similarities
    Sii::SparseMatrixCSC{Float64,Int} # item-item similarities
end

function memorize(
    Rui::SparseMatrixCSC{Float64,Int},
    config::Config = defaultconfig,
)::Memory
    """construct memory for recommendation"""
    nu, ni = size(Rui)
    Riu = sparse(Rui')
    if config.centering
        bu = biases(Riu)
        bi = biases(Rui)
        Dui = centering(Rui, bu)
        Diu = centering(Riu, bi)
    else
        bu = zeros(nu)
        bi = zeros(ni)
        Dui = Rui
        Diu = Riu
    end
    Sii = cossim(Dui)
    Suu = cossim(Diu)
    Memory(nu, ni, bu, bi, Dui, Diu, Suu, Sii)
end

function itembased_scores(
    m::Memory,
    users::Vector{Int},
    target_items::Union{Nothing,Vector{Int}} = nothing,
)::Matrix{Float64}
    """user-item scores based on item similarities"""
    items = target_items == nothing ? collect(1:m.ni) : target_items
    scores(m.bu, m.Dui, m.Sii, users, items)
end

function userbased_scores(
    m::Memory,
    users::Vector{Int},
    target_items::Union{Nothing,Vector{Int}} = nothing,
)::Matrix{Float64}
    """user-item scores based on user similarities"""
    items = target_items == nothing ? collect(1:m.ni) : target_items
    scores(m.bi, m.Diu, m.Suu, items, users)'
end

function itembased_rankings(
    m::Memory,
    k::Int,
    users::Vector{Int},
    target_items::Union{Nothing,Vector{Int}} = nothing,
)::Tuple{Matrix{Int}, Matrix{Float64}}
    """rank top-k items for users based on item similarities and return them with scores"""
    items = target_items == nothing ? collect(1:m.ni) : target_items
    @assert k <= length(items)
    scores = itembased_scores(m, users, items)
    perms = topkperm(scores, k)
    items[perms], selectbyrow(scores, perms)
end

function userbased_rankings(
    m::Memory,
    k::Int,
    users::Vector{Int},
    target_items::Union{Nothing,Vector{Int}} = nothing,
)::Tuple{Matrix{Int}, Matrix{Float64}}
    """rank top-k items for users based on user similarities and return them with scores"""
    items = target_items == nothing ? collect(1:m.ni) : target_items
    @assert k <= length(items)
    scores = userbased_scores(m, users, items)
    perms = topkperm(scores, k)
    items[perms], selectbyrow(scores, perms)
end

function scores(
    b::Vector{Float64},
    D::SparseMatrixCSC{Float64,Int},
    S::SparseMatrixCSC{Float64,Int},
    xs::Vector{Int},
    ys::Vector{Int},
)::Matrix{Float64}
    """calc scores of users x items"""
    b = b[xs]
    D = D[xs, :]
    S = S[:, ys]

    is, js, vs = findnz(D)
    I = sparse(is, js, ones(length(vs)))
    N = abs.(S)
    W = 1. ./ (I * N)
    W[isinf.(W)] .= 0.
    W = sparse(W)

    b .+ W .* (D * S)
end

function biases(R::SparseMatrixCSC{Float64,Int})::Vector{Float64}
    """biases for each columns"""
    sum(R, dims = 1)[:] ./ mapslices(nnz, R, dims = 1)[:]
end

function centering(
    R::SparseMatrixCSC{Float64,Int},
    b::Vector{Float64},
)::SparseMatrixCSC{Float64,Int}
    """centering each values"""
    xs, ys, vs = findnz(R)
    sparse(xs, ys, vs - b[xs])
end

function cossim(D::SparseMatrixCSC{Float64,Int})::SparseMatrixCSC{Float64,Int}
    """calc cosine similarity between each columns"""
    # normalizing weights for each columns
    w = 1. ./ sqrt.(sum(D .^ 2, dims = 1)[:])
    w[isinf.(w)] .= 0.
    w = sparsevec(w)

    # normalizing D
    wD = w' .* D

    # calc similarities
    wD' * wD
end

function topkperm(V::Matrix{Float64}, k::Int)::Matrix{Int}
    """select top-k indices by each row of matrix"""
    mapslices(vs -> partialsortperm(vs, 1:k, rev = true), V, dims = 2)
end

function selectbyrow(M::Matrix{Float64}, I::Matrix{Int})::Matrix{Float64}
    """select matrix contents using indices by each row"""
    vcat((M[i, :][I[i, :]]' for i in 1:size(I, 1))...)
end

ext: jl


text:
#!/usr/bin/env julia
# QT CLUSTERING IMPLEMENTATION
# AUTHOR: FREDERIK GADE

using LinearAlgebra: dot, mul!

const SAVE_CLUSTER_VALS = (2, 3, 6, 9, 11, 23, 47, 106, 235, 551, 699)
const HELP_MESSAGE = """usage: julia $PROGRAM_FILE inputfile threshold

QT (Quality Threshold) Clustering is an algorithm that groups multi-dimensional
vectors into high quality clusters. Quality is ensured by finding large cluster
whose diameter does not exceed a given user-defined diameter threshold.
This method prevents dissimilar vectors from being forced under the same
cluster and ensures that only good quality clusters will be formed.

arguments:
    inputfile\tInput file name (tab-separated list of floats w/ optional name
\t\tin first column)
    threshold\tQuality threshold of cluster diameter
"""
const NO_ARGS = "No arguments supplied. Use -h/--help for help."
const ONE_ARG = "Only one argument was supplied. Use -h/--help for help."
const TOO_MANY_ARGS = "Too many arguments supplied. Use -h/--help for help."

##############################
### INPUT/OUTPUT FUNCTIONS ###
##############################
function tab_split(ln)
    split(ln, "\t")
end

function dp(lnsplit::Array{SubString{String},1},
        offset::Int64, i::Int64)
    datapoint = [parse(Float64, x) for x in lnsplit[offset:end]]
    if offset === 1
        return string("Point", lpad(i, 4, "0")), datapoint
    else
        return string(lnsplit[1]), datapoint
    end
end

function read_points(f, num_points)
    i = 1
    firstline = tab_split(readline(f))
    if tryparse(Float64, firstline[1]) === nothing
        dim = length(firstline) - 1
        offset = 2
    else
        dim = length(firstline)
        offset = 1
    end

    datapoints = Matrix{Float64}(undef, dim, num_points)
    names = Vector{String}(undef, num_points)

    names[i], datapoints[:, i] = dp(firstline, offset, i)
    for ln in eachline(f)
        i += 1
        names[i], datapoints[:, i] = dp(tab_split(ln), offset, i)
    end

    names, datapoints
end

function validate_input(filename, _threshold)
    percentage = false
    num_points = countlines(filename)
    threshold = tryparse(Float64, _threshold)
    if threshold === nothing
        if _threshold[end] == '%'
            percentage = true
            _threshold = chop(_threshold)
            threshold = tryparse(Float64, _threshold)
            (threshold === nothing) &&
                throw(ArgumentError("Threshold has to be a number/percentage"))
            if threshold < 0
                throw(ArgumentError("Threshold has to be positive"))
            elseif threshold > 100
                @warn "Percentage given is over 100%."
            end
        end
    end
    threshold, percentage, num_points
end

function write_output(outfile::IO, i::Int64, best_cluster::BitSet, names::Array{String, 1},
    datapoints::Array{Float64,2})
    println(outfile, "-> Cluster ", i)
    for c in best_cluster
        print(outfile, names[c], "\t")
        println(outfile, join(view(datapoints, :, c), "\t"))
    end
end

###############################
### PREPROCESSING FUNCTIONS ###
###############################

# Distance matrix calculation is copied/modified from JuliaStats/Distances.jl
# https://github.com/JuliaStats/Distances.jl/blob/master/LICENSE.md
function sumsq_percol(a::Array{Float64,2})
    n = size(a, 2)
    r = Vector{Float64}(undef, n)
    @simd for j in 1:n
        aj = view(a, :, j)
        r[j] = dot(aj, aj)
    end
    return r
end

function _get_dists!(r::Array{Float64,2}, a::Array{Float64,2})
    m, n = size(a)
    mul!(r, a', a)
    sa2 = sumsq_percol(a)
    @inbounds for j = 1:n
        for i = 1:(j - 1)
            @inbounds r[i, j] = r[j, i]
        end
        r[j, j] = 0
        sa2j = sa2[j]
        @simd for i = (j + 1):n
            r[i, j] = sqrt(max(sa2[i] + sa2j - 2 * r[i, j], 0))
        end
    end
    r
end

function get_dists(datapoints::Array{Float64,2})
    n = size(datapoints, 2)
    r = Matrix{Float64}(undef, n, n)
    _get_dists!(r, datapoints)
    diameter = maximum(r)
    r, diameter
end
# End distance matrix calculation

function get_neighbours(dists::Array{Float64,2}, threshold::Float64,
        num_points::Int64)
    neighbours = [BitSet() for i = 1:num_points]
    for i = 1:num_points
        distsi = view(dists, :, i)
        @inbounds @simd for j = i+1:num_points
            @inbounds if distsi[j] <= threshold
                push!(neighbours[i], j)
                push!(neighbours[j], i)
            end
        end
    end
    neighbours
end

################################
### CORE ALGORITHM FUNCTIONS ###
################################

function get_best_cluster(unclustered, dists, neighbours, threshold,
        candidate_clusters, candidate_diameters, calculated_dict)
    best_cluster_size = zero(Int64)
    best_diameter = zero(eltype(candidate_diameters))
    best_cluster = BitSet()
    for seed in unclustered
        if !isempty(candidate_clusters[seed])
            candidate, candidate_diameter =
                candidate_clusters[seed], candidate_diameters[seed]
            candidate_size = length(candidate)
        else
            candidate, candidate_diameter =
                generate_candidate(seed, dists, neighbours, threshold,
                calculated_dict, candidate_clusters, candidate_diameters)
            candidate_clusters[seed] = candidate
            candidate_diameters[seed] = candidate_diameter
            candidate_size = length(candidate)
        end

        if candidate_size > best_cluster_size ||
                (candidate_size == best_cluster_size
                && candidate_diameter < best_diameter)
            best_cluster = candidate
            best_cluster_size = candidate_size
            best_diameter = candidate_diameter
        end
    end
    best_cluster
end

function generate_candidate(seed, dists, neighbours,
        threshold::Float64, calculated_dict, candidate_clusters, candidate_diameters)
    candidate_cluster = BitSet([seed])
    last_added = seed
    seed_neighbours = collect(copy(neighbours[seed]))
    diameter_cache = zeros(eltype(dists), length(seed_neighbours))
    candidate_diameter = zero(eltype(dists))

    while !isempty(seed_neighbours)
        if length(candidate_cluster) in SAVE_CLUSTER_VALS
            found_seed = find_cluster(candidate_cluster, calculated_dict, seed)
            if found_seed !== 0
                return candidate_clusters[found_seed], candidate_diameters[found_seed]
            end
        end

        candidate_point_index = update_diameter_cache!(diameter_cache, dists,
            seed_neighbours, last_added, threshold)
        candidate_point_value = seed_neighbours[candidate_point_index]

        if diameter_cache[candidate_point_index] < threshold
            candidate_diameter = diameter_cache[candidate_point_index]
            @inbounds push!(candidate_cluster, candidate_point_value)
            @inbounds popat!(seed_neighbours, candidate_point_index)
            @inbounds popat!(diameter_cache, candidate_point_index)
            last_added = candidate_point_value
        else
            return candidate_cluster, candidate_diameter
        end
    end

    if candidate_diameter === zero(Float64)
        return candidate_cluster, Inf
    end

    return candidate_cluster, candidate_diameter
end

function find_cluster(candidate_cluster::BitSet, calculated_dict, seed::Int64)
    calculated_index = hash(candidate_cluster)
    found_seed = get(calculated_dict, calculated_index, 0)
    if found_seed === 0
        calculated_dict[calculated_index] = seed
    end
    found_seed
end

function update_diameter_cache!(diameter_cache::Array{Float64,1},
        dists::Array{Float64,2}, seed_neighbours::Array{Int64,1},
        last_added::Int64, threshold::Float64)
    min_index = 1
    min_dist = threshold + 1
    @inbounds for (i, (cached_dist, new_dist)) in enumerate(zip(diameter_cache,
            view(dists, seed_neighbours, last_added)))
        @inbounds diameter_cache[i] = cached_dist > new_dist ? cached_dist : new_dist
        @inbounds (diameter_cache[i] < min_dist) &&
            ((min_index, min_dist) = (i, diameter_cache[i]))
    end
    min_index
end

function update_data!(unclustered::Array{Int64,1}, neighbours, candidate_clusters, best)
    setdiff!(unclustered, best)
    for unclustered_point in unclustered
        setdiff!(neighbours[unclustered_point], best)
        if !isdisjoint(candidate_clusters[unclustered_point], best)
            candidate_clusters[unclustered_point] = BitSet()
        end
    end
    nothing
end

######################
### MAIN ALGORITHM ###
######################
function QT(filename::String, _threshold::String)
    threshold, percentage, num_points = validate_input(filename, _threshold)
    names, datapoints = open(f -> read_points(f, num_points), filename)
    dists, diameter = get_dists(datapoints)

    if percentage
        threshold = threshold/100 * diameter
    elseif threshold > diameter
        @warn "Given threshold ($threshold) is greater than the diameter of the dataset ($diameter)."
    end

    neighbours = get_neighbours(dists, threshold, num_points)

    unclustered = collect(1:num_points)
    candidate_clusters = [BitSet() for i = 1:num_points]
    candidate_diameters = Vector{Float64}(undef, num_points)

    cluster_num = 1
    outfile = open("$filename.out", "w")

    while !isempty(unclustered)
        calculated_dict = Dict{UInt64,Int64}()
        best_cluster = get_best_cluster(unclustered, dists, neighbours,
            threshold, candidate_clusters, candidate_diameters, calculated_dict)

        write_output(outfile, cluster_num, best_cluster, names, datapoints)
        update_data!(unclustered, neighbours, candidate_clusters, best_cluster)

        cluster_num += 1
    end

    close(outfile)
end

#####################
### MAIN FUNCTION ###
#####################
function main()
    arg_len = length(ARGS)
    try
        if arg_len == 2
            @time QT(ARGS[1], ARGS[2])
            exit(0)
        elseif arg_len == 1
            if ARGS[1] in ("-h", "--help")
                println(HELP_MESSAGE)
                exit(0)
            else
                println(ONE_ARG)
            end
        elseif arg_len == 0
            println(NO_ARGS)
        else
            println(TOO_MANY_ARGS)
        end
        exit(1)
    catch e
        msg = sprint(showerror, e)
        println(msg)
        exit(1)
    end
end

main()


text:
using Documenter
using Autotuning

# The DOCSARGS environment variable can be used to pass additional arguments to make.jl.
# This is useful on CI, if you need to change the behavior of the build slightly but you
# can not change the .travis.yml or make.jl scripts any more (e.g. for a tag build).
if haskey(ENV, "DOCSARGS")
    for arg in split(ENV["DOCSARGS"])
        push!(ARGS, arg)
    end
end

makedocs(
    format = Documenter.HTML(prettyurls = !("local" in ARGS)),
    modules = [Autotuning],
    sitename = "Autotuning.jl",
    authors = "Pedro Bruel and contributors.",
    linkcheck = !("skiplinks" in ARGS),
    strict = !("strict=false" in ARGS),
    doctest = ("doctest=only" in ARGS) ? :only : true,
    pages = [
        "Home" => "index.md",
        "Library" => Any[
            "Public" => "lib/public.md"
        ],
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#