text:
struct FRStateDependentPolicy{S, F<:AbstractFloat} <: AbstractPolicy
    states::Array{S, 1}
    weights::Dict{S, Weights{F, F, Array{F, 1}}}
end

function Base.get(π::FRStateDependentPolicy, state_t::Array{Float64, 1}, action_t, state_tp1, action_tp1, preds_tp1)
    prob = -1.0
    prj = Int64.(floor.(state_t) .+ 1)
    if prj ∈ π.states
        prob = π.weights[prj][action_t]
    else
        prob = 0.25
    end
    return prob
end

function StatsBase.sample(rng::Random.AbstractRNG, π::FRStateDependentPolicy, state::Array{Float64, 1})
    act = -1
    prj = Int64.(floor.(state) .+ 1)
    if prj ∈ π.states
        act = StatsBase.sample(rng, 1:4, π.weights[prj])
    else
        act = rand(rng, 1:4)
    end
    return act
end

function FourCornersVariant(high_prob=0.45, low_prob=0.05)
    map = Dict(
        [3,2]=>Weights([high_prob, low_prob, low_prob, high_prob]),
        [3,9]=>Weights([high_prob, high_prob, low_prob, low_prob]),
        [10,2]=>Weights([low_prob, low_prob, high_prob, high_prob]),
        [10,9]=>Weights([low_prob, high_prob, high_prob, low_prob])
    )
    states = collect(keys(map))
    return FRStateDependentPolicy(states, map)
end

function find_same(states)
    s=0
    for p in collect(Iterators.product(states,states))
        s += if (p[1] == p[2]) 1 else 0 end
    end
    return (s - length(states))/2
end

# function gen_states(seed, num)
#     mt = Random.MersenneTwister(seed)
#     states = []
#     while length(states) < num
#         new_state = rand(mt, 1:11, 2)
#         if !(new_state in states)
#             append!(states, [new_state])
#         end
#     end
#     return states
# end

function gen_states(rng::Random.AbstractRNG, num)
    states = Array{Int64, 1}[]
    while length(states) < num
        new_state = rand(rng, 1:11, 2)
        if !(new_state in states)
            append!(states, [new_state])
        end
    end
    return states
end

gen_states(seed::Int64, num) = gen_states(Random.MersenneTwister(seed), num_states)

function RandomStateVariant(action, prob, num_states, seed::Int64)
    states = gen_states(Random.MersenneTwister(seed), num_states)
    w = Weights([act == action ? prob : (1.0-prob)/3
                 for act in 1:4])
    map = Dict([s=>w for s in states])
    return FRStateDependentPolicy(states, map)
end

function RandomStateVariant(action, prob, num_states, rng::Random.AbstractRNG)
    states = gen_states(rng, num_states)
    w = Weights([act == action ? prob : (1.0-prob)/3
                 for act in 1:4])
    map = Dict([s=>w for s in states])
    return FRStateDependentPolicy(states, map)
end

RandomStateVariant() =
    RandomStateVariant(FourRoomsParams.DOWN, 0.05, 25, 1143139448)

function random_weight_vector(rng::Random.AbstractRNG, num_actions)
    favored_action = rand(rng, 1:num_actions)
    prob = 0.5*(rand(rng) + 1.0)
    [act == favored_action ? prob : (1.0-prob)/(num_actions-1) for act in 1:num_actions]
end

function RandomStateWeightVariant(num_states::Int64, seed::Int64)
    rng = Random.MersenneTwister(seed)
    states = gen_states(rng, num_states)
    weights = Dict([state=>Weights(random_weight_vector(rng, 4)) for state in states])
    FRStateDependentPolicy(states, weights)
end

RandomStateWeightVariant() =
    RandomStateWeightVariant(25, 1143139448)


text:
__precompile__()
module SeisRemoveEQ

include("utils.jl")
include("map_removeEQ.jl")

using .Utils
using .Map_removeEQ

using Dates, Printf, JLD2, FileIO, Distributed

export seisremoveEQ

function seisremoveEQ(InputDict::Dict)

	Utils.initlogo()

	finame = InputDict["finame"]
	fodir  = InputDict["fodir"]
	foname = InputDict["foname"]

	mkpath(fodir)
	fopath=joinpath(fodir, foname*".jld2")
	InputDict["fopath"] = fopath

	#save data to fopath file
	t = jldopen(finame)
	jldopen(fopath, "w") do file
		file["info/DLtimestamplist"] = t["info/DLtimestamplist"];
		file["info/stationlist"]     = t["info/stationlist"];
		file["info/starttime"]       = t["info/starttime"];
		file["info/endtime"]         = t["info/endtime"];
		file["info/DL_time_unit"]    = t["info/DL_time_unit"];
	end

	DLtimestamplist = t["info/DLtimestamplist"];
	stationlist     = t["info/stationlist"];
	NumofTimestamp  = length(DLtimestamplist)
	JLD2.close(t)

	#print parameters
	printparams(InputDict)

	InputDict["DLtimestamplist"] = DLtimestamplist
	InputDict["stationlist"] = stationlist
	InputDict["NumofTimestamp"] = NumofTimestamp

	#evaluate memory use
	printstyled("---Start Test Process---\n"; color=:cyan, bold=true)

	max_num_of_processes_per_parallelcycle = get_memoryuse(InputDict)

	if max_num_of_processes_per_parallelcycle < 1
		error("Memory allocation is not enought (currently $MAX_MEM_PER_CPU [GB]). Please inclease MAX_MEM_PER_CPU or decrease number of stations")
	end

	printstyled("---Start removing EQ---\n"; color=:cyan, bold=true)

	# parallelize process by time stamp

	file = jldopen(fopath, "r+")

	E = []
	bt_getkurtosis = 0.0
	bt_removeeq = 0.0

	if max_num_of_processes_per_parallelcycle >= NumofTimestamp
		# one process cycle cover all time stamp
		EE = pmap(x -> map_removeEQ(x, InputDict), 1:NumofTimestamp)

		for i = 1:size(EE)[1]
			push!(E, EE[i][1])
			bt_getkurtosis += EE[i][2]
		    bt_removeeq += EE[i][3]
		end

		# save data to jld2
		for ii = 1:size(E)[1] #loop at each starttime
			for jj = 1:size(E[1])[1] #loop at each station id

			   requeststr =E[ii][jj].id
			   varname = joinpath(DLtimestamplist[ii], requeststr)
			   #save_SeisData2JLD2(fopath, varname, S[ii][jj])
			   file[varname] = E[ii][jj]
			end
		end

	else
		# processes by each max_num_of_processes_per_parallelcycle
		pitr = 1

		while pitr <=  NumofTimestamp

			E = []

		    startid1 = pitr
		    startid2 = pitr + max_num_of_processes_per_parallelcycle - 1

		    if startid1 == NumofTimestamp
		        #use one
		        EE = pmap(x -> map_removeEQ(x, InputDict), startid1:startid1)

		    elseif startid2 <= NumofTimestamp
		        #use all processors
		        EE = pmap(x -> map_removeEQ(x, InputDict), startid1:startid2)

		    else
		        #use part of processors
		        startid2 = startid1 + mod(NumofTimestamp, max_num_of_processes_per_parallelcycle) - 1
		        println(startid2)
		        EE = pmap(x -> map_removeEQ(x, InputDict), startid1:startid2)
		    end

			for i = 1:size(EE)[1]
				push!(E, EE[i][1])
				bt_getkurtosis += EE[i][2]
			    bt_removeeq += EE[i][3]
			end

		    # save data to jld2
		    for ii = 1:size(E)[1] #loop at each starttime
		        for jj = 1:size(E[1])[1] #loop at each station id

		            requeststr =E[ii][jj].id
		            varname = joinpath(DLtimestamplist[startid1+ii-1], requeststr)
		            #save_SeisData2JLD2(fopath, varname, S[ii][jj])
		            file[varname] = E[ii][jj]

		        end
		    end

		    pitr += max_num_of_processes_per_parallelcycle

		end
	end

	JLD2.close(file)

	printstyled("---Summary---\n"; color=:cyan, bold=true)
	println("time to get kurtosis =$(bt_getkurtosis)[s]")
	println("time to remove EQ    =$(bt_removeeq)[s]")
	print("\nprocess done at:")
	println(now())
	str = "EQ is successfully removed from raw data:\nOutput = $fopath\n"
	printstyled(str; color=:green, bold=true)

end

end # module

ext: jl


text:
using FiniteDifferences, ChainRulesCore

# Make FiniteDifferences work with some of the types in this package. Shame this isn't
# automated...

import FiniteDifferences: to_vec

# Make FiniteDifferences work for any user-defined type that has fields. This assumes that
# the type passed in admits a constructor that just takes all of the fields and produces a
# new object. As such, this won't work in the general case, but should work for all of the
# stuff in this package.
function to_vec(x::T) where {T}
    isempty(fieldnames(T)) && throw(error("Expected some fields. None found."))
    vecs_and_backs = map(name->to_vec(getfield(x, name)), fieldnames(T))
    vecs, backs = first.(vecs_and_backs), last.(vecs_and_backs)
    x_vec, back = to_vec(vecs)
    return x_vec, function(x′_vec)
        vecs′ = back(x′_vec)
        return T(map((back, vec)->back(vec), backs, vecs′)...)
    end
end

# My version of isapprox
function fd_isapprox(x_ad::Nothing, x_fd, rtol, atol)
    return fd_isapprox(x_fd, zero(x_fd), rtol, atol)
end
function fd_isapprox(x_ad::AbstractArray, x_fd::AbstractArray, rtol, atol)
    return all(fd_isapprox.(x_ad, x_fd, rtol, atol))
end
function fd_isapprox(x_ad::Real, x_fd::Real, rtol, atol)
    return isapprox(x_ad, x_fd; rtol=rtol, atol=atol)
end
function fd_isapprox(x_ad::NamedTuple, x_fd, rtol, atol)
    f = (x_ad, x_fd)->fd_isapprox(x_ad, x_fd, rtol, atol)
    return all([f(getfield(x_ad, key), getfield(x_fd, key)) for key in keys(x_ad)])
end
function fd_isapprox(x_ad::Tuple, x_fd::Tuple, rtol, atol)
    return all(map((x, x′)->fd_isapprox(x, x′, rtol, atol), x_ad, x_fd))
end
function fd_isapprox(x_ad::Dict, x_fd::Dict, rtol, atol)
    return all([fd_isapprox(get(()->nothing, x_ad, key), x_fd[key], rtol, atol) for
        key in keys(x_fd)])
end

function adjoint_test(
    f, ȳ, x...;
    rtol=1e-9,
    atol=1e-9,
    fdm=FiniteDifferences.Central(5, 1),
    print_results=false,
    test=true,
)
    # Compute forwards-pass and j′vp.
    y, back = Zygote.pullback(f, x...)
    adj_ad = back(ȳ)
    adj_fd = j′vp(fdm, f, ȳ, x...)

    # Check that forwards-pass agrees with plain forwards-pass.
    test && @test y ≈ f(x...)

    # Check that ad and fd adjoints (approximately) agree.
    print_results && print_adjoints(adj_ad, adj_fd, rtol, atol)
    test && @test fd_isapprox(adj_ad, adj_fd, rtol, atol)

    return adj_ad, adj_fd
end

function print_adjoints(adjoint_ad, adjoint_fd, rtol, atol)
    @show typeof(adjoint_ad), typeof(adjoint_fd)
    adjoint_ad, adjoint_fd = to_vec(adjoint_ad)[1], to_vec(adjoint_fd)[1]
    println("atol is $atol, rtol is $rtol")
    println("ad, fd, abs, rel")
    abs_err = abs.(adjoint_ad .- adjoint_fd)
    rel_err = abs_err ./ adjoint_ad
    display([adjoint_ad adjoint_fd abs_err rel_err])
    println()
end

function ensure_not_running_on_functor(f, name)
    # if x itself is a Type, then it is a constructor, thus not a functor.
    # This also catchs UnionAll constructors which have a `:var` and `:body` fields
    f isa Type && return

    if fieldcount(typeof(f)) > 0
        throw(ArgumentError(
            "$name cannot be used on closures/functors (such as $f)"
        ))
    end
end

"""
    frule_test(f, (x, ẋ)...; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), kwargs...)
# Arguments
- `f`: Function for which the `frule` should be tested.
- `x`: input at which to evaluate `f` (should generally be set to an arbitary point in the domain).
- `ẋ`: differential w.r.t. `x` (should generally be set randomly).
All keyword arguments except for `fdm` are passed to `isapprox`.
"""
function frule_test(f, (x, ẋ); rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), kwargs...)
    return frule_test(f, ((x, ẋ),); rtol=rtol, atol=atol, fdm=fdm, kwargs...)
end

function frule_test(f, xẋs::Tuple{Any, Any}...; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), kwargs...)
    ensure_not_running_on_functor(f, "frule_test")
    xs, ẋs = collect(zip(xẋs...))
    Ω, pushforward = ChainRulesCore.frule(f, xs...)
    @test fd_isapprox(f(xs...), Ω, rtol, atol)
    dΩ_ad = pushforward(NamedTuple(), ẋs...)

    # Correctness testing via finite differencing.
    dΩ_fd = jvp(fdm, xs->f(xs...), (xs, ẋs))
    @test isapprox(
        collect(dΩ_ad),  # Use collect so can use vector equality
        collect(dΩ_fd);
        rtol=rtol,
        atol=atol,
        kwargs...
    )
end


text:
###########################
### Basic Organizations ###
###########################

listMyOrgs(c::Client) = getParsedResponse(Vector{Organization},c,Requests.get,"/user/orgs")

listUserOrgs(c::Client, user::String) = getParsedResponse(Vector{Organization},c,Requests.get,"/users/$(user)/orgs")

getOrg(c::Client,org::String) = getParsedResponse(Organization,c,Requests.get,"/orgs/$(org)")

FieldTags.@tag immutable CreateOrgOption
	orgName::String => json:"username"
	fullName::Nullable{String} => json:"full_name"
	description::Nullable{String}
	website::Nullable{String}
	location::Nullable{String}
end

function CreateOrgOption(username; kwargs...)
	kwargsDict = Dict{Symbol,Any}(k=>v for (k,v) in kwargs)

	args = Any[username]

	for field in fieldnames(CreateOrgOption)[2:end]
		push!(args,get(kwargsDict,field, fieldtype(CreateOrgOption,field)()))
	end

	CreateOrgOption(args...)
end

FieldTags.@tag immutable EditOrgOption
	fullName::Nullable{String} => json:"full_name"
	description::Nullable{String}
	website::Nullable{String}
	location::Nullable{String}
end

function EditOrgOption(;kwargs...)
	kwargsDict = Dict{Symbol,Any}(k=>v for (k,v) in kwargs)

	args = Any[]

	for field in fieldnames(EditOrgOption)[1:end]
		push!(args,get(kwargsDict,field, fieldtype(EditOrgOption,field)()))
	end

	EditOrgOption(args...)
end

editOrg(c::Client,orgname::String,opt::EditOrgOption) = getParsedResponse(Organization,c,Requests.patch,"/orgs/$(orgname)"; json = marshalJSON(opt))

################################
### Organizations Membership ###
################################

## TODO: Bug Gitea Devs for member ship API examples

addOrgMembership(c::Client,orgname::String,user::String,role::String) = ( error("Gitea API broken"); getResponse(c,Requests.put,"/orgs/$(orgname)/members/$(user)"; json = Dict("role"=>role)); nothing)


text:
function StructuredQueries._collect(src::AbstractTable, q::StructuredQueries.SelectNode)
    res = default(src)
    for h in q.helpers
        apply!(res, h, src)
    end
    return res
end

function apply!(res, h::SQ.SelectHelper, src)::Void
    res_field, f, arg_fields, = SQ.parts(h)
    row_itr = eachrow(src, arg_fields...)
    inner_eltypes = map(
        eltype,
        eltypes(src, arg_fields...)
    )

    # TODO: check for non-leaf types and follow slow path
    T = Core.Inference.return_type(f, (Tuple{inner_eltypes...},))
    res_column = NullableVector{T}()

    _apply!(res_column, h, row_itr)
    res[res_field] = res_column
    return
end

function _apply!(res_column, h::SQ.SelectHelper, row_itr)::Void
    T = eltype(res_column)
    f = h.f
    for row in row_itr
        if hasnulls(row)
            push!(res_column, Nullable{T}())
        else
            v = f(map(unsafe_get, row))
            push!(res_column, v)
        end
    end
    return
end


text:
function integrate(m1, m2, p0, e0)
    """
    General purpose integrator for eqn(5.14) of Peters, 1964 [1].

    [1] https://ui.adsabs.harvard.edu/abs/1964PhRv..136.1224P/abstract

    Params:
        m1 -- The mass of the first star  (Solar Mass)
        m2 -- The mass of the second star (Solar Mass)
        p0 -- The initial period          (days)
        e0 -- The initial eccentricity

    Returns:
        tC -- The coalescence time        (seconds)
    """
    Solar_Mass   = 1.989e30    # kg
    Solar_Radius = 696340000.0 # m
               G = 6.67e-11    # m^3 kg^-1 s^-2
               c = 299792458.0 # m/s

    convert(Float64, m1)
    convert(Float64, m2)
    convert(Float64, p0)
    convert(Float64, e0)

    ##########################################
    ############ BEGIN UNIT CHECK ############
    ##########################################
    m1 = m1 * Solar_Mass   # Kilograms     ###
    m2 = m2 * Solar_Mass   # Kilograms     ###
    p0 = p0 * 24 * 60 * 60 # Seconds       ###
    e0 = e0                # Dimensionless ###
    ##########################################
    ############# END UNIT CHECK #############
    ##########################################

    # Use Kepler's Third Law to compute the semimajor axis, in meters
    a0 = (p0^2.0 * (G * (m1+m2)) / (4.0 * pi^2.0))^(1.0/3.0) # meters

    # compute the constant c0
    c0 = a0 * (1.0-e0^2.0) * e0^(-12.0/19.0) * (1.0+(121.0 / 304.0 * e0^2.0))^(-870.0/2299.0) # meters

    # Compute the constant beta
    beta = 64.0 / 5.0 * G^3.0 * m1 * m2 * (m1 + m2) / c^5.0
    # meters^4 / s

    # Circular binary coalescence time
    tC = a0^4.0 / (4.0*beta)

    if e0 != 0
        if e0 < 0.01
            # Low eccentricity - see eqn after eqn(5.14) of Peters, 1964
            tC = c0^4.0 * e0^(48.0/19.0) / (4.0*beta)
        elseif e0 > 0.99
            # High eccentricity - see eqn after eqn after eqn(5.14) of Peters, 1964.
            tC = tC * ((768.0 / 425.0) * (1.0-e0^2.0)^3.5)
        else
            # Medium eccentricity - see eqn(5.14) of Peters, 1964
            e = 0.0
            de = e0 / 10000.0
            summand = 0.0

            while e < e0
                this_integral  = de * e^(29.0/19.0)
                this_integral *= (1.0 + (121.0/304.0) * e^2.0)^(1181.0/2299.0)
                this_integral /= (1.0 - e^2.0)^(1.5)

                summand += this_integral
                e       += de

                convert(Float64, e)
            end

            tC = (12.0/19.0) * (c0^4.0 / beta) * summand
        end
    end

    return tC
end

max_stars_repo_name: Krytic/Takahe

max_issues_repo_name: Krytic/Takahe

max_forks_repo_name: Krytic/Takahe


text:
# This script is the translation of run_example.m in
# Dylan Mikesell's DynamicWarping repo (https://github.com/dylanmikesell/DynamicWarping.)
# run this script in the terminal: '> julia DTWDT_run_example.jl'
# 2019/06 Kurama Okubo

using PlotlyJS, ORCA, FileIO, JLD2, Printf, DTWDT
#for dv/v analysis
using DataFrames, GLM

fodir = "./EXAMPLE/fig/stretched"
mkpath(fodir)

## load the data file and plot
@load "./exampledata/StretchedData.jld2"; # data with sine wave shifts of one cycle

# plot shift function
maxLag = 80; # max nuber of points to search forward and backward (can be npts, just takes longer and is unrealistic)
b      = 1; # b-value to limit strain

lvec   = (-maxLag:maxLag).*dt; # lag array for plotting below
npts   = length(u0);            # number of samples
tvec   = ( 0 : npts-1 ) .* dt; # make the time axis
stTime = st;

## compute timeshift
stbarTime, stbar, dist, dtwerror = dtwdt(u0, u1, dt,
                    dtwnorm="L2",
                    maxLag=80, #number of maxLag id to search the distance
                    b=10, # b value to controll distance calculation algorithm
                    direction=1, #direction to accumulate errors (1=forward, -1=backward, 0=double to smooth)
                    )

tvec2 = tvec + stbarTime; # make the warped time axis

# plot the distance function
function plotdistfunc()
    tr = heatmap(;x=tvec, y=lvec, z=dist, colorscale="Jet",
    zauto=true, colorbar=attr(title="Dist", titleside="right",))
    layout = Layout(
        title = "Distance function",
        xaxis=attr(title="Time [s]"),
        yaxis=attr(title="τ [s]"),
        width=800,
        height=400,
    )
    plot(tr, layout)
end

p = (plotdistfunc())
display(p)

savefig(p, fodir*"/distance_array.png")

function lineplot1()
    tr1 = scatter(;x=tvec, y=stTime, mode="lines", name="Actual")
    tr2 = scatter(;x=tvec, y=stbarTime, mode="markers", name="Estimated")
    layout = Layout(
        title = "Estimated shifts",
        xaxis=attr(title="Time [s]"),
        yaxis=attr(title="τ [s]"),
        width=800,
        height=400,
    )
    plot([tr1, tr2], layout)
end

p = (lineplot1())
display(p)
savefig(p, fodir*"/inputsignal.png")

# plot input traces
function lineplot_s1()
    tr1 = scatter(;x=tvec, y=u0, mode="lines", name="Raw")
    tr2 = scatter(;x=tvec, y=u1, mode="lines", name="Shifted",
    line = attr(dash="dot"))

    layout = Layout(
        title = "Input traces for dynamic time warping",
        xaxis=attr(title="Time [s]"),
        yaxis=attr(title="Amplitude (a.u.)"),
    )
    plot([tr1, tr2], layout)
end

function lineplot_s2()
    tr1 = scatter(;x=tvec, y=u0, mode="lines", name="Raw")
    tr2 = scatter(;x=tvec2, y=u1, mode="lines", name="Warped",
    line = attr(dash="dot"))

    layout = Layout(
        title = "Output traces for dynamic time warping",
        xaxis=attr(title="Time [s]"),
        yaxis=attr(title="Amplitude (a.u.)"),
    )
    plot([tr1, tr2], layout)
end

p = ([lineplot_s1(), lineplot_s2()])
display(p)
savefig(p, fodir*"/traces.png")

# dv/v analysis
#linear regression of stbarTime (=estimated time shift)

data = DataFrame(X=tvec, Y=stbarTime)
ols = lm(@formula(Y ~ X), data)

dtt = GLM.coef(ols)[2]
dvv = dtt # assume dt/T = dv/V

linstTime = collect(tvec) .* dtt .+ GLM.coef(ols)[1]
#plot line
# plot input traces

function lineplot_ln()
    tr1 = scatter(;x=tvec, y=stTime, mode="lines", name="Actual")
    tr2 = scatter(;x=tvec, y=stbarTime, mode="markers", name="Estimated")
    tr3 = scatter(;x=tvec, y=linstTime, mode="lines", name="Linear regression",
                    line=attr(width=5.0, color="red"))
    layout = Layout(
        title = @sprintf("true dv/v: -0.01, estimated dv/v = %4.8f", dvv),
        xaxis=attr(title="Time [s]"),
        yaxis=attr(title="τ [s]"),
        width = 800,
        height = 400,
    )
    plot([tr1, tr2, tr3], layout)
end

p = (lineplot_ln())
display(p)
savefig(p, fodir*"/linregression.png")

# plot schematic of dynamic warping
yshift0 = maximum(abs.(u0))
yshift1 = maximum(abs.(u1))

u0s = u0 .+ yshift0
u1s = u1 .- yshift1

layout = Layout(width=800, height=600,
				xaxis=attr(title="Time [s]", dtick=10.0),
				yaxis=attr(title="Amplitude (a.u.)"),
                font =attr(size=18),
                title = "Schematic of dynamic time warping")
p = plot([NaN], layout)

#plot connection
linespan = 5
for i = 1:linespan:npts
    x0 = tvec[i]
    y0 = u0s[i]
    x1 = tvec[stbar[i]+i]
    y1 = u1s[stbar[i]+i]

    trace1 = scatter(;x=[x0, x1], y=[y0, y1], mode="lines",
        line_color="rgb(200, 200, 200)",
        line=attr(dash = false,),
        showlegend=false,
        )
    addtraces!(p, trace1)
end

#plot two traces
tr1 = scatter(;x=tvec, y=u0s, mode="lines", name="Raw",)
tr2 = scatter(;x=tvec, y=u1s, mode="lines", name="Shifted",)
addtraces!(p, tr1)
addtraces!(p, tr2)

deletetraces!(p, 1)
savefig(p, fodir*"/schematic.png")

println("Type something to exit this example.");
readline();


text:
module Sarkses

## =============================================================================
module SkewSuffixArrays

## Implementation of skew algorithm for construction of suffix array
## based on code provided in:
## "Kärkkäinen, J., & Sanders, P. (2003, June). Simple linear work suffix
##  array construction. In International Colloquium on Automata, Languages,
##  and Programming (pp. 943-955). Springer, Berlin, Heidelberg."

function radixPass!(a::Array{Int64}, b::Array{Int64}, r::Array{Int64},
                    rShift::Int64, n::Int64, K::Int64)::Nothing
    Kp1 = K+1
    ## count occurrences
    c = zeros(Int64, Kp1)
    @inbounds for i = 1:n
        c[r[rShift + a[i] + 1] + 1] += 1
    end
    ## exclusive prefix sums
    sum = 0
    for i = 1:Kp1
        t = c[i]
        c[i] = sum
        sum += t
    end
    ## sort
    @inbounds for i = 1:n
        b[c[r[rShift + a[i] + 1] + 1] + 1] = a[i]
        c[r[rShift + a[i] + 1] + 1] += 1
    end
    return nothing
end

function suffixArray!(s::Array{Int64}, SA::Array{Int64}, n::Int64, K::Int64)::Nothing
    n0=(n+2)÷3; n1=(n+1)÷3; n2=n÷3; n02=n0+n2
    s12 = zeros(Int64, n02+3)
    SA12 = zeros(Int64, n02+3)
    s0 = Array{Int64}(undef, n0)
    SA0 = Array{Int64}(undef, n0)
    ##
    ## generate positions of mod 1 and mod 2 suffixes
    ## the "+(n0-n1)" adds a dummy mod 1 suffix if n%3 == 1
    let j = 1
        @inbounds for i = 0:(n+(n0-n1)-1)
            if i%3 != 0
                s12[j] = i
                j += 1
            end
        end
    end
    ##
    ## lsb radix sort the mod 1 and mod 2 triples
    radixPass!(s12, SA12, s, 2, n02, K)
    radixPass!(SA12, s12, s, 1, n02, K)
    radixPass!(s12, SA12, s, 0, n02, K)
    ##
    ## find lexicographic names of triples
    let name = 0, c0 = -1, c1 = -1, c2 = -1
        @inbounds for i = 1:n02
            if s[SA12[i]+1] != c0 || s[SA12[i]+2] != c1 || s[SA12[i]+3] != c2
                name += 1
                c0 = s[SA12[i]+1]
                c1 = s[SA12[i]+2]
                c2 = s[SA12[i]+3]
            end
            if (SA12[i] % 3) == 1
                ## left half
                s12[SA12[i]÷3 + 1] = name
            else
                ## right half
                s12[SA12[i]÷3 + n0+1] = name
            end
        end
        ##
        ## recurse if names are not yet unique
        if name < n02
            suffixArray!(s12, SA12, n02, name)
            ## store unique names in s12 using the suffix array
            @inbounds for i = 1:n02
                s12[SA12[i]+1] = i
            end
        else
            ## generate the suffix array of s12 directly
            @inbounds for i = 1:n02
                SA12[s12[i]] = i-1
            end
        end
    end
    ##
    ## stably sort the mod 0 suffixes from SA12 by their first character
    let j = 1
        @inbounds for i = 1:n02
            if SA12[i] < n0
                s0[j] = 3 * SA12[i]
                j += 1
            end
        end
    end
    radixPass!(s0, SA0, s, 0, n0, K)
    ##
    ## merge sorted SA0 suffixes and sorted SA12 suffixes
    let p=0; t=n0-n1; k=0
        while k < n
            ## pos of current offset 12 suffix
            i = (SA12[t+1] < n0 ? SA12[t+1] * 3 + 1 : (SA12[t+1] - n0) * 3 + 2)
            ## pos of current offset 0 suffix
            j = SA0[p+1]
            ##
            if (SA12[t+1] < n0 ?
                (s[i+1], s12[SA12[t+1]+n0+1]) < (s[j+1], s12[j÷3+1]) :
                (s[i+1], s[i+2], s12[SA12[t+1]-n0+2]) < (s[j+1], s[j+2], s12[j÷3+n0+1]))
                ## suffix from SA12 is smaller
                SA[k+1] = i
                t += 1
                if t == n02
                    ## done --- only SA0 suffixes left
                    k += 1
                    while p < n0
                        SA[k+1] = SA0[p+1]
                        p += 1
                        k += 1
                    end
                end
            else
                SA[k+1] = j
                p += 1
                if p == n0
                    ## done --- only SA12 suffixes left
                    k += 1
                    while t < n02
                        SA[k+1] = (SA12[t+1] < n0 ?
                                   SA12[t+1] * 3 + 1 :
                                   (SA12[t+1] - n0) * 3 + 2)
                        t += 1
                        k += 1
                    end
                end
            end
            k += 1
        end
    end
    return nothing
end

end # module SkewSuffixArrays

## =============================================================================
module FMIndices

## Implementation of Burrows-Wheeler transform for purpose of exact
## k-mer matching based on backward search algorithm described in:
## "Ferragina, P., & Manzini, G. (2000, November). Opportunistic
##  data structures with applications. In Proceedings 41st Annual
##  Symposium on Foundations of Computer Science (pp. 390-398). IEEE."

## Implementation of inexact k-mer matching done following procedure
## building on above reference described in:
## "Heng Li, Richard Durbin (2009, July). Fast and accurate short read
##  alignment with Burrows–Wheeler transform. Bioinformatics, Volume 25,
##  Issue 14 (pp. 1754–1760). Oxford University Press."

function fillFMIndexCO!(bwt::Array{UInt8},
                        c::Dict{UInt8,Int64},
                        o::Dict{UInt8,Array{Int64}},
                        seq::Array{UInt8},
                        sa::Array{Int64},
                        blockSize::Int64)::Nothing
    oSize = ((length(bwt)-1) ÷ blockSize) + 1
    @inbounds for i=1:length(bwt)
        if sa[i] > 1
            ichar = seq[ sa[i]-1 ]
        else
            ichar = b"!"[1]
        end
        bwt[i] = ichar
        if !(ichar in keys(c))
            c[ichar] = 0
        end
        c[ichar] += 1
        if !(ichar in keys(o))
            o[ichar] = zeros(Int64, oSize)
        end
        o[ichar][1 + ((i-1) ÷ blockSize)] += 1
    end
    cumul = 0
    cKeys = sort(collect(keys(c)))
    @inbounds for ichar in cKeys
        cumul += c[ichar]
        c[ichar] = cumul - c[ichar]
        ochar = o[ichar]
        for i=2:oSize
            ochar[i] += ochar[i-1]
        end
    end
    return nothing
end

mutable struct FMIndex
    seq::Array{UInt8}
    sa::Array{Int64}
    bwt::Array{UInt8}
    c::Dict{UInt8,Int64}
    o::Dict{UInt8,Array{Int64}}

    blockSize::Int64

    function FMIndex(seq::Array{UInt8}, sa::Array{Int64})
        this = new()
        this.seq = seq
        this.sa = sa
        this.bwt = Array{UInt8}(undef, length(this.seq))
        this.c = Dict{UInt8,Int64}()
        this.o = Dict{UInt8,Int64}()
        this.blockSize = 100
        fillFMIndexCO!(this.bwt, this.c, this.o,
                       this.seq, this.sa, this.blockSize)
        return this
    end
end

function handleBlocking(bwt::Array{UInt8},
                        c::Dict{UInt8,Int64},
                        o::Dict{UInt8,Array{Int64}},
                        i::Int64,
                        ichar::UInt8,
                        blockSize::Int64)::Int64
    iRounded = ((i-1) ÷ blockSize)
    iNew = c[ichar] + (iRounded > 0 ? o[ichar][iRounded] : 0) + 1
    for j=(1+(blockSize*iRounded)):(i-1)
        if bwt[j] == ichar
            iNew += 1
        end
    end
    return iNew
end

function findKmer(this::FMIndex, kmer::Array{UInt8})::Tuple{Int64,Int64}
    k = length(kmer)
    start = 1
    end_ = length(this.bwt) + 1
    for i=k:-1:1
        ichar = kmer[i]
        start = handleBlocking(
            this.bwt, this.c, this.o, start, ichar, this.blockSize
        )
        end_ = handleBlocking(
            this.bwt, this.c, this.o, end_, ichar, this.blockSize
        )
    end
    return (start, end_)
end

function calculateD(rev::FMIndex, kmer::Array{UInt8})::Array{Int64}
    k = length(kmer)
    d = zeros(k)
    start = 1
    end_ = length(rev.bwt) + 1
    z = 0
    for i=1:k
        ichar = kmer[i]
        start = handleBlocking(rev.bwt, rev.c, rev.o, start, ichar,
                               rev.blockSize)
        end_ = handleBlocking(rev.bwt, rev.c, rev.o, end_, ichar,
                              rev.blockSize)
        if end_ <= start
            start = 1
            end_ = length(rev.bwt) + 1
            z += 1
        end
        d[i] = z
    end
    return d
end

function findKmer(this::FMIndex, kmer::Array{UInt8},
                  nMismatch::Int64, d::Array{Int64},
                  i::Int64, start::Int64, end_::Int64;
                  indels::Bool=false)::Array{Tuple{Int64,Int64}}
    if i < 1
        if nMismatch >= 0
            return [(start, end_)]
        else
            return Array{Tuple{Int64,Int64}}(undef, 0)
        end
    end
    if nMismatch < d[i]
        return Array{Tuple{Int64,Int64}}(undef, 0)
    end
    I = Array{Tuple{Int64,Int64}}(undef, 0)
    if indels
        append!(I, findKmer(
            this, kmer, nMismatch-1, d, i-1, start, end_, indels=true
        ))
    end
    ichar = kmer[i]
    for bchar in sort(collect(keys(this.c)))
        bstart = handleBlocking(this.bwt, this.c, this.o, start, bchar,
                                this.blockSize)
        bend = handleBlocking(this.bwt, this.c, this.o, end_, bchar,
                              this.blockSize)
        if bend > bstart
            if indels
                append!(I, findKmer(
                    this, kmer, nMismatch-1, d, i, bstart, bend, indels=true
                ))
            end
            if bchar == ichar
                append!(I, findKmer(
                    this, kmer, nMismatch, d, i-1, bstart, bend, indels=indels
                ))
            else
                append!(I, findKmer(
                    this, kmer, nMismatch-1, d, i-1, bstart, bend, indels=indels
                ))
            end
        end
    end
    return I
end

function findKmer(this::FMIndex, rev::FMIndex, kmer::Array{UInt8},
                  nMismatch::Int64; indels::Bool=false)::Array{Tuple{Int64,Int64}}
    d = calculateD(rev, kmer)
    return findKmer(this, kmer, nMismatch, d,
                    length(kmer), 1, length(this.bwt)+1, indels=indels)
end

end # module FMIndices

using CodecZlib
using CSV
using DataFrames
using DataStructures
using LinearAlgebra
using Random
using Statistics
using TranscodingStreams

function multifind(needles::Array{<:AbstractString},
                   haystack::AbstractArray{<:AbstractString}
                   )::Array{Int64}
    out = -1 * ones(Int64, length(needles))
    ln = length(needles); lh = length(haystack)
    nsp = sortperm(needles); hsp = sortperm(haystack)
    i = 1
    j = 1
    while (i <= ln) && (j <= lh)
        if haystack[hsp[j]] > needles[nsp[i]]
            i += 1
        elseif haystack[hsp[j]] == needles[nsp[i]]
            out[nsp[i]] = hsp[j]
            j += 1
        else
            j += 1
        end
    end
    return out
end

function c1idx(df::DataFrame, idx::Array{<:AbstractString})::DataFrame
    ndx = multifind(idx, df[:, 1])
    return df[ndx[ndx .> 0], :]
end

function c1idx(od::OrderedDict, idx::Array)::OrderedDict
    return OrderedDict([k => od[k] for k in idx])
end

function padAndUnstack(stacked::DataFrame,
                       rowkey::Union{Integer,Symbol}=:id,
                       colkey::Union{Integer,Symbol}=:bin,
                       value::Union{Integer,Symbol}=:score;
                       side::Symbol=:left,
                       fill::Union{Nothing,Float64}=nothing,
                       n::Union{Nothing,Int64}=nothing
                       )::DataFrame
    if n == nothing
        n = maximum(stacked[:, colkey])
    end
    stacked = deepcopy(stacked)
    stacked = stacked[:, [rowkey, colkey, value]]
    uniqueIds = unique(stacked[:, rowkey])
    toAdds = Array{DataFrame}(undef, 0)
    idToInds = Dict{String,Array{Int64}}()
    for id in uniqueIds
        idToInds[id] = Array{Int64}(undef, 0)
    end
    for i in 1:size(stacked, 1)
        push!(idToInds[stacked[i, 1]], i)
    end
    for id in uniqueIds
        rmatch = idToInds[id]
        maxBin = maximum(stacked[rmatch, colkey])
        if maxBin < n
            if side == :left
                stacked[rmatch, colkey] .+= (n - maxBin)
            end
            toAdd = DataFrame(r=[id for i in 1:(n-maxBin)], c=1:(n-maxBin))
            if side != :left
                toAdd.c .+= maxBin
            end
            vVal = if fill == nothing
                mean(stacked[rmatch, value])
            else
                fill
            end
            toAdd[:, :v] = [vVal for i in 1:size(toAdd, 1)]
            rename!(toAdd, [names(toAdd)[i] =>
                            Symbol(names(stacked)[i]) for i in 1:3])
            push!(toAdds, toAdd)
        end
    end
    if length(toAdds) > 0
        stacked = vcat(stacked, vcat(toAdds...))
    end
    stacked = stacked[sortperm(stacked[:, colkey]), :]
    return unstack(stacked, rowkey, colkey, value)
end

function readFasta(filename::String)::OrderedDict{String,String}
    out = OrderedDict{String,String}()
    curName = ""
    sb = IOBuffer()
    fgz = (endswith(lowercase(filename), ".gz") ||
           endswith(lowercase(filename), ".gzip"))
    fh = (fgz ?
          TranscodingStream(GzipDecompressor(), open(filename, "r")) :
          open(filename, "r"))
    while !eof(fh)
        curLine = rstrip(readline(fh))
        if startswith(curLine, ">")
            if curName != ""
                out[curName] = String(take!(sb))
            end
            curName = curLine[2:length(curLine)]
        else
            write(sb, curLine)
        end
    end
    if curName != "" && !(curName in keys(out))
        out[curName] = String(take!(sb))
    end
    close(fh)
    close(sb)
    return out
end

## =============================================================================
abstract type AbSarks end

function concatenateSeqs!(this::AbSarks,
                          seqs::AbstractDict{<:AbstractString,<:AbstractString}
                          )::Nothing
    sb = IOBuffer()
    this.bounds = Array{Int64}(undef, length(this.blocks)+1)
    this.bounds[1] = 1
    let i = 2
        for block in this.blocks
            blockSeq = String(seqs[block])
            write(sb, blockSeq)
            write(sb, "\$")
            this.bounds[i] = this.bounds[i-1] + length(blockSeq) + 1
            i += 1
        end
    end
    this.catSeq = take!(sb)
    close(sb)
    return nothing
end

function calcSuffixArray!(this::AbSarks)::Nothing
    chars = convert.(Int64, [this.catSeq[i] for i in 1:length(this.catSeq)])
    append!(chars, [0, 0, 0])
    SA = zeros(Int64, length(chars))
    SkewSuffixArrays.suffixArray!(chars, SA, length(this.catSeq), maximum(chars))
    this.sa = SA[1:length(this.catSeq)] .+ 1
    this.saInv = Array{Int64}(undef, length(this.sa))
    for i = 1:length(this.sa)
        this.saInv[ this.sa[i] ] = i
    end
    return nothing
end

## -----------------------------------------------------------------------------
function reset!(this::AbSarks, halfWindow::Int64,
                spatialLength::Int64=0, smooth::Bool=true)::Nothing
    this.halfWindow = halfWindow
    this.spatialLength = spatialLength
    windowGini!(this)
    if smooth
        window!(this)
    else
        this.windowed = Array{Float32}(undef, 0)
    end
    if this.spatialLength > 1
        spatialGini!(this)
        if smooth
            spatialWindow!(this)
        else
            this.spatialWindowed = Array{Float32}(undef, 0)
        end
    else
        this.spatGini = Array{Float32}(undef, 0)
        this.spatialWindowed = Array{Float32}(undef, 0)
    end
    return nothing
end

function resetSpatial!(this::AbSarks, spatialLength::Int64,
                       smooth::Bool=true)::Nothing
    this.spatialLength = spatialLength
    if this.spatialLength > 1
        spatialGini!(this)
        if smooth
            spatialWindow!(this)
        else
            this.spatialWindowed = Array{Float32}(undef, 0)
        end
    else
        this.spatGini = Array{Float32}(undef, 0)
        this.spatialWindowed = Array{Float32}(undef, 0)
    end
    return nothing
end

## -----------------------------------------------------------------------------
function sourceBlock(this::AbSarks, s::Int64,
                     mindex::Int64=1, maxdex::Int64=-1)::Int64
    if maxdex == -1
        maxdex = length(this.bounds) + 1
    elseif maxdex == (mindex + 1)
        return mindex
    end
    middex = (mindex + maxdex) ÷ 2
    if s < this.bounds[middex]
        return sourceBlock(this, s, mindex, middex)
    else
        return sourceBlock(this, s, middex, maxdex)
    end
end

function sourceBlock(this::AbSarks, s0::Array{Int64})::Array{Int64}
    maxdex = length(this.bounds) + 1
    return [sourceBlock(this, s, 1, maxdex) for s in s0]
end

function fillSourceBlock!(b::Array{Int64},
                          blocks::Array{String},
                          bounds::Array{Int64},
                          saInv::Array{Int64})::Nothing
    @inbounds for block in 1:length(blocks)
        nextBlock = block + 1
        for s in bounds[block]:(bounds[nextBlock]-1)
            b[ saInv[s] ] = block
        end
    end
end

function sourceBlock(this::AbSarks)::Array{Int64}
    b = Array{Int64}(undef, length(this.sa))
    fillSourceBlock!(b, this.blocks, this.bounds, this.saInv)
    return b
end

## -----------------------------------------------------------------------------
function giniImpurity(counts::Array{Int64})::Float64
    total = sum(counts)
    out = sum(convert(Array{Float64}, counts) .* (total .- counts))
    return(out / (total * total))
end

function fillSpatGini!(spatGini::Array{Float32},
                       spatialLength::Int64, saLen::Int64,
                       windGini::Array{Float32}, saInv::Array{Int64},
                       runningSum::Float64, windSize::Float64)::Float64
    @inbounds for s=2:(saLen+1-spatialLength)
        runningSum -= windGini[ saInv[s-1] ]
        runningSum += windGini[ saInv[s+spatialLength-1] ]
        spatGini[ saInv[s] ] = convert(Float32, runningSum / windSize)
    end
    return runningSum
end

function spatialGini!(this::AbSarks)::Nothing
    windSize = convert(Float64, this.spatialLength)
    this.spatGini = Array{Float32}(undef, length(this.sa))
    runningSum = 0.0
    for s = 1:this.spatialLength
        runningSum += this.windGini[ this.saInv[s] ]
    end
    this.spatGini[ this.saInv[1] ] =
            convert(Float32, runningSum / windSize)
    runningSum = fillSpatGini!(this.spatGini, this.spatialLength,
                               length(this.sa), this.windGini,
                               this.saInv, runningSum, windSize)
    for s = (length(this.sa)+2-this.spatialLength):length(this.sa)
        this.spatGini[ this.saInv[s] ] =
                convert(Float32, runningSum / windSize)
    end
    return nothing
end

## -----------------------------------------------------------------------------
function fillSpatialWindowed!(spatialWindowed::Array{Float32},
                              spatialLength::Int64, saLen::Int64,
                              windowed::Array{Float32}, saInv::Array{Int64},
                              runningSum::Float64, windSize::Float64)::Float64
    @inbounds for s = 2:(saLen+1-spatialLength)
        runningSum -= windowed[ saInv[s-1] ]
        runningSum += windowed[ saInv[s+spatialLength-1] ]
        spatialWindowed[ saInv[s] ] = convert(Float32, runningSum / windSize)
    end
    return runningSum
end

function spatialWindow!(this::AbSarks)::Nothing
    windSize = convert(Float64, this.spatialLength)
    this.spatialWindowed = Array{Float32}(undef, length(this.sa))
    runningSum = 0.0
    for s = 1:this.spatialLength
        runningSum += this.windowed[ this.saInv[s] ]
    end
    this.spatialWindowed[ this.saInv[1] ] =
            convert(Float32, runningSum / windSize)
    runningSum = fillSpatialWindowed!(this.spatialWindowed, this.spatialLength,
                                      length(this.sa), this.windowed,
                                      this.saInv, runningSum, windSize)
    for s = (length(this.sa)+2-this.spatialLength):length(this.sa)
        this.spatialWindowed[ this.saInv[s] ] =
                convert(Float32, runningSum / windSize)
    end
    return nothing
end

## -----------------------------------------------------------------------------
function findKmer(this::AbSarks, kmer::Array{UInt8})::Tuple{Int64,Int64}
    return FMIndices.findKmer(this.fm, kmer)
end

function kmer(this::AbSarks, s::Int64, k::Int64, k0::Int64=0)::Array{UInt8}
    return this.catSeq[max(1, s+k0):min(length(this.catSeq), s+k-1)]
end

function kmers(this::AbSarks, s::Array{Int64}, k::Int64,
               k0::Int64=0)::Array{Array{UInt8}}
    return [kmer(this, sel, k, k0) for sel in s]
end

function prefixAgreeSum(this::AbSarks, i::Int64, kmax::Int64)::Float64
    kmer_ = kmer(this, this.sa[i], kmax)
    kmax = min(kmax, length(kmer_))
    agreeSum = 0
    wstart = i - this.halfWindow
    wend = i + this.halfWindow + 1
    let kstart=i, kend=i+1
        for k = kmax:-1:1
            kwin = Sarkses.findKmer(this, kmer_[1:k])
            kwin1 = max(wstart, kwin[1])
            kwin2 = min(wend, kwin[2])
            agreeSum += k * ((kstart - kwin1) +
                             (kwin2 - kend))
            kstart = kwin1
            kend = kwin2
        end
    end
    return convert(Float64, agreeSum) / (2.0 * this.halfWindow)
end

## -----------------------------------------------------------------------------
function filterGuts(sa::Array{Int64}, saInv::Array{Int64},
                    windowed::Array{Float32}, windGini::Array{Float32},
                    spatialWindowed::Array{Float32}, spatGini::Array{Float32},
                    theta::Float64=-Inf, minGini::Float64=-Inf,
                    spatialTheta::Float64=-Inf, minSpatialGini::Float64=-Inf;
                    peakify::Bool=true)::Array{Int64}
    if minGini >= 1.0
        minGini = 1.0 - (1.0 - median(windGini)) * minGini
    end
    if minSpatialGini >= 1.0
        minSpatialGini = 1.0 - (1.0 - median(windGini)) * minSpatialGini
    end
    pos = Array{Int64}(undef, 0)
    keep = false
    iLeft = -1
    iRight = -1
    @inbounds for i = 1:length(sa)
        if windowed[i] >= theta && windGini[i] >= minGini
            keep = true
            s = sa[i]
            if spatialTheta > -Inf
                if (spatialWindowed[i] < spatialTheta ||
                    spatGini[i] < minSpatialGini)
                    keep = false
                end
            end
            if keep && peakify
                iLeft = saInv[1]
                if s > 1
                    iLeft = saInv[s-1]
                end
                iRight = saInv[length(saInv)]
                if s < length(saInv)
                    iRight = saInv[s+1]
                end
                if ((windowed[iLeft] > windowed[i]) ||
                    (windowed[iRight] > windowed[i]))
                    keep = false
                end
            end
            if keep
                push!(pos, i)
            end
        end
    end
    return pos
end

function filter(this::AbSarks, theta::Float64=-Inf, minGini::Float64=-Inf,
                spatialTheta::Float64=-Inf, minSpatialGini::Float64=-Inf;
                peakify::Bool=true)::Array{Int64}
    return filterGuts(this.sa, this.saInv, this.windowed, this.windGini,
                      this.spatialWindowed, this.spatGini, theta, minGini,
                      spatialTheta, minSpatialGini, peakify=peakify)
end

function filter(this::AbSarks,
                filters::Array{Dict{String,Float64},1},
                thresholds::Array{Float64,2}=Array{Float64,2}(undef, 0, 0);
                peakify::Bool=true)::Array{Array{Int64,1},1}
    hw0 = this.halfWindow
    sl0 = this.spatialLength
    out = Array{Array{Int64,1},1}(undef, 0)
    for f=1:length(filters)
        filt = filters[f]
        hw = convert(Int64, filt["halfWindow"])
        ming = ("minGini" in keys(filt) ? filt["minGini"] : -Inf)
        sl = convert(Int64, filt["spatialLength"])
        minsg = ("minSpatialGini" in keys(filt) ? filt["minSpatialGini"] : -Inf)
        theta = -Inf
        spatialTheta = -Inf
        if size(thresholds)[1] > 0
            theta = thresholds[f, 1]
            spatialTheta = thresholds[f, 2]
        elseif "theta" in keys(filt)
            theta = filt["theta"]
        elseif "spatialTheta" in keys(filt)
            spatialTheta = filt["spatialTheta"]
        end
        if hw != this.halfWindow
            reset!(this, hw, sl, true)
        elseif sl != this.spatialLength
            resetSpatial!(this, sl, true)
        end
        push!(out, filter(
            this, theta, ming,
            spatialTheta, minsg, peakify=peakify
        ))
    end
    if this.halfWindow != hw0
        reset!(this, hw0, sl0, true)
    elseif this.spatialLength != sl0
        resetSpatial!(this, sl0, true)
    end
    return out
end

## -----------------------------------------------------------------------------
function kmerPeaks(this::AbSarks,
                   filters::Array{Dict{String,Float64},1},
                   thresholds::Array{Float64,2},
                   iFilt::Array{Array{Int64,1},1};
                   kmax::Int64=12)::DataFrame
    hw0 = this.halfWindow
    sl0 = this.spatialLength
    out = DataFrame(
        i=Array{Int64}(undef,0), s=Array{Int64}(undef,0), kmer=Array{String}(undef,0),
        khat=Array{Float64}(undef,0), block=Array{String}(undef,0),
        wi=Array{Int64}(undef,0), gini=Array{Float32}(undef,0),
        spatialGini=Array{Float32}(undef,0), score=Array{Float64}(undef,0),
        windowed=Array{Float32}(undef,0), spatialWindowed=Array{Float32}(undef,0),
        kmax=Array{Int64}(undef,0),
        halfWindow=Array{Int64}(undef,0), minGini=Array{Float64}(undef,0),
        theta=Array{Float64}(undef,0), spatialLength=Array{Int64}(undef,0),
        minSpatialGini=Array{Float64}(undef,0), spatialTheta=Array{Float64}(undef,0)
    )
    hasCatScores = :catScores in fieldnames(typeof(this))
    for f = 1:length(iFilt)
        filt = filters[f]
        hw = convert(Int64, filt["halfWindow"])
        ming = ("minGini" in keys(filt) ? filt["minGini"] : -Inf)
        sl = convert(Int64, filt["spatialLength"])
        minsg = ("minSpatialGini" in keys(filt) ? filt["minSpatialGini"] : -Inf)
        if hw != this.halfWindow
            reset!(this, hw, sl, true)
        elseif sl != this.spatialLength
            resetSpatial!(this, sl, true)
        end
        for i in iFilt[f]
            s = this.sa[i]
            khat = prefixAgreeSum(this, i, kmax)
            km = kmer(this, s, kmax)[1:convert(Int64, round(khat))]
            blockIndex = sourceBlock(this, s)
            block = this.blocks[blockIndex]
            wi = s - this.bounds[blockIndex] + 1
            gini = this.windGini[i]
            sg = (length(this.spatGini) > 0 ? this.spatGini[i] : -Inf32)
            score = (hasCatScores ? this.catScores[i] : this.scores[blockIndex])
            sw = (length(this.spatialWindowed) > 0 ?
                  this.spatialWindowed[i] : -Inf32)
            append!(out, Dict(
                "i"=>i, "s"=>s, "kmer"=>String(km), "khat"=>khat,
                "block"=>block, "wi"=>wi, "gini"=>gini, "spatialGini"=>sg,
                "score"=>score, "windowed"=>this.windowed[i],
                "spatialWindowed"=>sw, "kmax"=>kmax,
                "halfWindow"=>hw, "minGini"=>ming,
                "theta"=>thresholds[f, 1], "spatialLength"=>sl,
                "minSpatialGini"=>minsg, "spatialTheta"=>thresholds[f, 2]
            ))
        end
    end
    reset!(this, hw0, sl0, true)
    return out
end

function spatialSubPeaks(this::AbSarks,
                         iFilt::Array{Int64,1},
                         theta::Float64)::Array{Int64}
    if this.spatialLength <= 1
        return [this.sa[i] for i in iFilt]
    end
    subPeaks = Set{Int64}()
    for i in iFilt
        sRangeMax = min(this.sa[i] + this.spatialLength - 1, length(this.saInv))
        union!(subPeaks, collect(this.sa[i]:sRangeMax))
    end
    filteredOut = Set{Int64}()
    for s in subPeaks
        if this.windowed[this.saInv[s]] < theta
            push!(filteredOut, s)
        elseif (!((s-1) in subPeaks)) &&
               (this.windowed[this.saInv[s-1]] >= theta)
            push!(filteredOut, s)
        end
    end
    setdiff!(subPeaks, filteredOut)
    return sort(collect(subPeaks))
end

function mergeKmerIntervals(this::AbSarks,
                            subPeaks::Array{Int64};
                            kmax::Int64=12)::Array{Tuple{Int64,Int64},1}
    subPeaks = sort(subPeaks)
    left = 0; right = 0; sLast = -1
    out = Array{Tuple{Int64,Int64},1}(undef, 0)
    for s in subPeaks
        khat = convert(Int64, round(prefixAgreeSum(this, this.saInv[s], kmax)))
        if (s-1) == sLast
            right = min(max(right, s+khat-1), length(this.catSeq))
        else
            if left >= 1
                push!(out, (left, right))
            end
            left = s
            right = min(s+khat-1, length(this.catSeq))
        end
        sLast = s
    end
    if left >= 1
        push!(out, (left, right))
    end
    return out
end

function multiMergeKmerIntervals(this::AbSarks,
                                 iFilt::Array{Array{Int64,1},1},
                                 thresholds::Array{Float64,2},
                                 filters::Array{Dict{String,Float64},1};
                                 kmax::Int64)::Array{Array{Tuple{Int64,Int64},1},1}
    out = Array{Array{Tuple{Int64,Int64},1},1}(undef, 0)
    hw0 = this.halfWindow
    sl0 = this.spatialLength
    for f = 1:length(filters)
        filt = filters[f]
        hw = convert(Int64, filt["halfWindow"])
        sl = convert(Int64, filt["spatialLength"])
        if hw != this.halfWindow
            reset!(this, hw, sl, true)
        elseif sl != this.spatialLength
            resetSpatial!(this, sl, true)
        end
        subPeaks = spatialSubPeaks(this, iFilt[f], thresholds[f, 2])
        push!(out, mergeKmerIntervals(this, subPeaks, kmax=kmax))
    end
    reset!(this, hw0, sl0, true)
    return out
end

function mergedKmerSubPeaks(
        this::AbSarks,
        filters::Array{Dict{String,Float64},1},
        thresholds::Array{Float64,2};
        peakify::Bool=true,
        kmax::Int64=12
    )::DataFrame
    eyes = filter(this, filters, thresholds, peakify=peakify)
    mergedKmerIntervals = multiMergeKmerIntervals(
        this, eyes, thresholds, filters, kmax=kmax
    )
    hw0 = this.halfWindow
    sl0 = this.spatialLength
    out = DataFrame(
        i=Array{Int64}(undef,0), s=Array{Int64}(undef,0), kmer=Array{String}(undef,0),
        khat=Array{Float64}(undef,0), block=Array{String}(undef,0),
        wi=Array{Int64}(undef,0), gini=Array{Float32}(undef,0),
        spatialGini=Array{Float32}(undef,0), score=Array{Float64}(undef,0),
        windowed=Array{Float32}(undef,0), spatialWindowed=Array{Float32}(undef,0),
        kmax=Array{Int64}(undef,0),
        halfWindow=Array{Int64}(undef,0), minGini=Array{Float64}(undef,0),
        theta=Array{Float64}(undef,0), spatialLength=Array{Int64}(undef,0),
        minSpatialGini=Array{Float64}(undef,0), spatialTheta=Array{Float64}(undef,0)
    )
    for f = 1:length(filters)
        filt = filters[f]
        hw = convert(Int64, filt["halfWindow"])
        ming = ("minGini" in keys(filt) ? filt["minGini"] : -Inf)
        sl = convert(Int64, filt["spatialLength"])
        minsg = ("minSpatialGini" in keys(filt) ? filt["minSpatialGini"] : -Inf)
        theta = -Inf
        spatialTheta = -Inf
        if size(thresholds)[1] > 0
            theta = thresholds[f, 1]
            spatialTheta = thresholds[f, 2]
        elseif "theta" in keys(filt)
            theta = filt["theta"]
        elseif "spatialTheta" in keys(filt)
            spatialTheta = filt["spatialTheta"]
        end
        if hw != this.halfWindow
            reset!(this, hw, sl, true)
        elseif sl != this.spatialLength
            resetSpatial!(this, sl, true)
        end
        for sInterval in mergedKmerIntervals[f]
            i = this.saInv[sInterval[1]]
            km = this.catSeq[sInterval[1]:sInterval[2]]
            khat = prefixAgreeSum(this, i, kmax)
            blockIndex = sourceBlock(this, sInterval[1])
            block = this.blocks[blockIndex]
            wi = sInterval[1] - this.bounds[blockIndex] + 1
            gini = this.windGini[i]
            sg = (length(this.spatGini) > 0 ? this.spatGini[i] : -Inf32)
            score = this.scores[blockIndex]
            sw = (length(this.spatialWindowed) > 0 ?
                  this.spatialWindowed[i] : -Inf32)
            append!(out, Dict(
                "i"=>i, "s"=>sInterval[1], "kmer"=>String(km), "khat"=>khat,
                "block"=>block, "wi"=>wi, "gini"=>gini, "spatialGini"=>sg,
                "score"=>score, "windowed"=>this.windowed[i],
                "spatialWindowed"=>sw, "kmax"=>kmax,
                "halfWindow"=>hw, "minGini"=>ming,
                "theta"=>thresholds[f, 1], "spatialLength"=>sl,
                "minSpatialGini"=>minsg, "spatialTheta"=>thresholds[f, 2]
            ))
        end
    end
    reset!(this, hw0, sl0, true)
    return out
end

function mergedKmers(this::AbSarks,
                     mergedKmerIntervals::Array{Tuple{Int64,Int64},1}
                     )::Array{String}
    merged = Set{String}()
    for interval in mergedKmerIntervals
        push!(merged, this.catSeq[interval[1], interval[2]])
    end
    return sort(collect(merged))
end

function mergedKmers(this::AbSarks,
                     mergedKmerIntervals::Array{Array{Tuple{Int64,Int64},1}}
                     )::Array{Array{String}}
    return [mergedKmers(this, mergedKmerInts)
            for mergedKmerInts in mergedKmerIntervals]
end

## =============================================================================
mutable struct Sarks <: AbSarks
    halfWindow::Int64
    spatialLength::Int64

    scores::Array{Float64}
    blocks::Array{String}
    blockPosition::Dict{String,Int64}
    catSeq::Array{UInt8}
    bounds::Array{Int64}
    sa::Array{Int64}
    saInv::Array{Int64}
    windGini::Array{Float32}
    spatGini::Array{Float32}
    windowed::Array{Float32}
    spatialWindowed::Array{Float32}

    fm::FMIndices.FMIndex

    origScores::Array{Float64}
    charFreqs::Matrix{Float64}
    compositionPreBeta::Matrix{Float64}

    function Sarks(fasta::String, scoreFile::String,
                   halfWindow::Int64, spatialLength::Int64=0;
                   adjustForCharFreqs::Bool=false)
        this = new()
        this.halfWindow = halfWindow
        this.spatialLength = spatialLength
        seqs = readFasta(fasta)
        scoreMap = readScores(scoreFile)
        this.blocks = sort(intersect(
            collect(keys(scoreMap)),
            collect(keys(seqs))
        ))
        this.blockPosition = Dict{String,Int64}()
        for i in 1:length(this.blocks)
            this.blockPosition[this.blocks[i]] = i
        end
        this.origScores = [scoreMap[this.blocks[i]] for i in 1:length(this.blocks)]
        concatenateSeqs!(this, seqs)
        cf = characterFrequencies(this.bounds, this.catSeq)
        for i = 1:size(cf, 2)
            cf[:, i] .-= mean(cf[:, i])
        end
        if adjustForCharFreqs
            this.charFreqs = cf
            this.compositionPreBeta = inv(transpose(cf) * cf) * transpose(cf)
            deltaScores = (this.origScores .- mean(this.origScores))
            hatScores = cf * (this.compositionPreBeta * deltaScores)
            this.scores = deltaScores - hatScores
        else
            this.charFreqs = zeros(0, 0)
            this.compositionPreBeta = zeros(0, 0)
            this.scores = this.origScores
        end
        calcSuffixArray!(this)
        this.fm = FMIndices.FMIndex(this.catSeq, this.sa)
        windowGini!(this)
        window!(this)
        if this.spatialLength > 1
            spatialGini!(this)
            spatialWindow!(this)
        else
            this.spatGini = Array{Float32}(undef, 0)
            this.spatialWindowed = Array{Float32}(undef, 0)
        end
        return this
    end

    function Sarks(seqs::AbstractDict{<:AbstractString,<:AbstractString},
                   scoreMap::AbstractDict{<:AbstractString,Float64},
                   halfWindow::Int64, spatialLength::Int64=0;
                   adjustForCharFreqs::Bool=false)
        this = new()
        this.halfWindow = halfWindow
        this.spatialLength = spatialLength
        this.blocks = sort(intersect(
            collect(String.(keys(scoreMap))),
            collect(String.(keys(seqs)))
        ))
        this.blockPosition = Dict{String,Int64}()
        for i in 1:length(this.blocks)
            this.blockPosition[this.blocks[i]] = i
        end
        this.origScores = [scoreMap[this.blocks[i]] for i in 1:length(this.blocks)]
        concatenateSeqs!(this, seqs)
        cf = characterFrequencies(this.bounds, this.catSeq)
        for i = 1:size(cf, 2)
            cf[:, i] .-= mean(cf[:, i])
        end
        if adjustForCharFreqs
            this.charFreqs = cf
            this.compositionPreBeta = inv(transpose(cf) * cf) * transpose(cf)
            deltaScores = (this.origScores .- mean(this.origScores))
            hatScores = cf * (this.compositionPreBeta * deltaScores)
            this.scores = deltaScores - hatScores
        else
            this.charFreqs = zeros(0, 0)
            this.compositionPreBeta = zeros(0, 0)
            this.scores = this.origScores
        end
        calcSuffixArray!(this)
        this.fm = FMIndices.FMIndex(this.catSeq, this.sa)
        windowGini!(this)
        window!(this)
        if this.spatialLength > 1
            spatialGini!(this)
            spatialWindow!(this)
        else
            this.spatGini = Array{Float32}(undef, 0)
            this.spatialWindowed = Array{Float32}(undef, 0)
        end
        return this
    end

    function Sarks(sarks::Sarks)
        this = new()
        this.halfWindow = sarks.halfWindow
        this.spatialLength = sarks.spatialLength
        this.blocks = sarks.blocks
        this.blockPosition = sarks.blockPosition
        this.scores = sarks.scores
        this.catSeq = sarks.catSeq
        this.bounds = sarks.bounds
        this.origScores = sarks.origScores
        this.charFreqs = sarks.charFreqs
        this.compositionPreBeta = sarks.compositionPreBeta
        this.sa = sarks.sa
        this.saInv = sarks.saInv
        this.fm = sarks.fm
        this.windGini = sarks.windGini
        this.spatGini = sarks.spatGini
        this.windowed = sarks.windowed
        this.spatialWindowed = sarks.spatialWindowed
        return this
    end

    function Sarks(sarks::Sarks, permutation::Array{Int64})
        this = new()
        this.halfWindow = sarks.halfWindow
        this.spatialLength = sarks.spatialLength
        this.blocks = sarks.blocks
        this.blockPosition = sarks.blockPosition
        this.scores = [sarks.scores[ permutation[b] ]
                       for b in 1:length(sarks.scores)]
        this.catSeq = sarks.catSeq
        this.bounds = sarks.bounds
        this.origScores = [sarks.origScores[ permutation[b] ]
                           for b = 1:length(sarks.origScores)]
        this.charFreqs = sarks.charFreqs
        this.compositionPreBeta = sarks.compositionPreBeta
        if size(this.charFreqs, 1) > 0
            deltaScores = (this.origScores .- mean(this.origScores))
            hatScores = this.charFreqs * (this.compositionPreBeta * deltaScores)
            this.scores = deltaScores - hatScores
        else
            this.scores = this.origScores
        end
        this.sa = sarks.sa
        this.saInv = sarks.saInv
        this.fm = sarks.fm
        this.windGini = sarks.windGini
        this.spatGini = sarks.spatGini
        window!(this)
        if this.spatialLength > 1
            spatialWindow!(this)
        else
            this.spatialWindowed = Array{Float32}(undef, 0)
        end
        return this
    end
end

function readScores(filename::String)::DataStructures.OrderedDict{String,Float64}
    ordKeys = Array{String}(undef, 0)
    unord = Dict{String,Float64}()
    fgz = (endswith(lowercase(filename), ".gz") ||
           endswith(lowercase(filename), ".gzip"))
    fh = (fgz ?
          TranscodingStream(GzipDecompressor(), open(filename, "r")) :
          open(filename, "r"))
    while !eof(fh)
        tokens = split(rstrip(readline(fh)), "\t")
        score = tryparse(Float64, tokens[2])
        if score !== nothing
            unord[tokens[1]] = score
            push!(ordKeys, tokens[1])
        end
    end
    out = DataStructures.OrderedDict{String,Float64}()
    for key in ordKeys
        out[key] = unord[key]
    end
    return out
end

function characterFrequencies(bounds::Array{Int64},
                              catSeq::Array{UInt8})::Matrix{Float64}
    alphabet = setdiff(unique(catSeq), UInt8('$'))
    charFreqs = zeros(length(bounds)-1, length(alphabet)-1)
    @inbounds for i = 1:size(charFreqs, 1)
        @inbounds for j = 1:size(charFreqs, 2)
            charFreqs[i, j] =
                    sum(catSeq[bounds[i]:(bounds[i+1]-1)] .== alphabet[j]) /
                    (bounds[i+1] - bounds[i] - 1)
            ## the -1 above accounts for the end-of-sequence marker ($)
        end
    end
    return charFreqs
end

## -----------------------------------------------------------------------------
function sourceScore(this::Sarks, s::Int64)::Float64
    return this.scores[sourceBlock(this, s)]
end

function sourceScore(this::Sarks, s::Array{Int64})::Array{Float64}
    return [this.scores[b] for b in sourceBlock(this, s)]
end

function fillSourceScore!(ss::Array{Float64},
                          scores::Array{Float64},
                          srcBlk::Array{Int64})::Nothing
    @inbounds for i = 1:length(srcBlk)
        ss[i] = scores[ srcBlk[i] ]
    end
    return nothing
end

function sourceScore(this::Sarks)::Array{Float64}
    srcBlk = sourceBlock(this)
    out = Array{Float64}(undef, length(srcBlk))
    fillSourceScore!(out, this.scores, srcBlk)
    return out
end

## -----------------------------------------------------------------------------
function blockCounts(this::Sarks, i::Int64, block::Array{Int64})::Array{Int64}
    out = zeros(Int64, length(this.blocks))
    for j = -this.halfWindow:this.halfWindow
        out[ block[i+j] ] += 1
    end
    return out
end

function fillWindGini!(windGini::Array{Float32},
                       halfWindow::Int64, saLen::Int64,
                       block::Array{Int64}, bCounts::Array{Int64},
                       runningGini::Float64)::Nothing
    total = (2 * halfWindow) + 1
    total2 = convert(Float64, total * total)
    dGini = 0
    oldb = -1
    newb = -1
    @inbounds for i = (halfWindow+2):(saLen-halfWindow)
        oldb = block[i - 1 - halfWindow]
        newb = block[i + halfWindow]
        dGini = -bCounts[oldb] * (total - bCounts[oldb])
        if bCounts[newb] > 0
            dGini -= bCounts[newb] * (total - bCounts[newb])
        end
        bCounts[oldb] -= 1
        bCounts[newb] += 1
        if bCounts[oldb] > 0
            dGini += bCounts[oldb] * (total - bCounts[oldb])
        end
        dGini += bCounts[newb] * (total - bCounts[newb])
        runningGini += (convert(Float64, dGini) / total2)
        windGini[i] = convert(Float32, runningGini)
    end
end

function windowGini!(this::Sarks)::Nothing
    block = sourceBlock(this)
    this.windGini = Array{Float32}(undef, length(this.sa))
    bCounts = blockCounts(this, this.halfWindow+1, block)
    runningGini = giniImpurity(bCounts)
    this.windGini[this.halfWindow+1] = convert(Float32, runningGini)
    fillWindGini!(this.windGini, this.halfWindow, length(this.sa),
                  block, bCounts, runningGini)
    for i = 1:this.halfWindow
        this.windGini[i] = this.windGini[this.halfWindow+1]
        this.windGini[length(this.sa) + 1 - i] =
                this.windGini[length(this.sa) - this.halfWindow]
    end
    return nothing
end

## -----------------------------------------------------------------------------
function fillWindowed!(windowed::Array{Float32},
                       halfWindow::Int64, saLen::Int64,
                       scores::Array{Float64}, srcBlk::Array{Int64},
                       runningSum::Float64, windSize=Float64)::Float64
    @inbounds for i = (halfWindow+2):(saLen-halfWindow)
        runningSum -= scores[ srcBlk[i-1-halfWindow] ]
        runningSum += scores[ srcBlk[i+halfWindow] ]
        windowed[i] = convert(Float32, runningSum / windSize)
    end
    return runningSum
end

function window!(this::Sarks)::Nothing
    windSize = convert(Float64, (2 * this.halfWindow) + 1)
    srcBlk = sourceBlock(this)
    this.windowed = Array{Float32}(undef, length(this.sa))
    runningSum = 0.0
    for i = 1:((2*this.halfWindow)+1)
        runningSum += this.scores[ srcBlk[i] ]
    end
    this.windowed[this.halfWindow + 1] = convert(Float32, runningSum / windSize)
    runningSum = fillWindowed!(this.windowed, this.halfWindow,
                               length(this.sa), this.scores, srcBlk,
                               runningSum, windSize)
    for i = 1:this.halfWindow
        this.windowed[i] = this.windowed[this.halfWindow+1]
        this.windowed[length(this.sa) + 1 - i] =
                convert(Float32, runningSum / windSize)
    end
    return nothing
end

## -----------------------------------------------------------------------------
function giniMask(this::Sarks,
                  minGini::Float64,
                  minSpatialGini::Float64)::AbstractArray{Bool}
    if minGini >= 1
        minGini = 1 - (1-median(this.windGini))*minGini
    end
    if minSpatialGini >= 1
        minSpatialGini = 1 - (1-median(this.windGini))*minSpatialGini
    end
    mask = trues(length(this.windGini))
    if minGini > 0
        for i = 1:length(mask)
            if this.windGini[i] < minGini
                mask[i] = false
            end
        end
    end
    if minSpatialGini > 0 && this.spatialLength > 1
        for i = 1:length(mask)
            if this.spatGini[i] < minSpatialGini
                mask[i] = false
            end
        end
    end
    return mask
end

function permutationDistribution(
        this::Sarks,
        reps::Int64,
        filters::Array{Dict{String,Float64},1};
        seed::Int64=0,
        permutations::Array{Int64,2}=Array{Int64,2}(undef, 0, 0)
    )::Array{Float32,3}
    rng = (seed != 0 ? MersenneTwister(seed) : MersenneTwister())
    if size(permutations)[1] == 0
        permutations = Array{Int64,2}(undef, reps, length(this.scores))
        for r = 1:reps
            permutations[r, :] = randperm(rng, length(this.scores))
        end
    end
    hw0 = this.halfWindow
    sl0 = this.spatialLength
    permResults = Array{Float32,3}(undef, reps, length(filters), 2)
    for f = 1:length(filters)
        filt = filters[f]
        hw = convert(Int64, filt["halfWindow"])
        ming = ("minGini" in keys(filt) ? filt["minGini"] : -Inf)
        sl = convert(Int64, filt["spatialLength"])
        minsg = ("minSpatialGini" in keys(filt) ? filt["minSpatialGini"] : -Inf)
        if hw != this.halfWindow
            reset!(this, hw, sl, false)
        elseif sl != this.spatialLength
            resetSpatial!(this, sl, false)
        end
        fmask = giniMask(this, ming, minsg)
        Threads.@threads for r = 1:reps
            rsarks = Sarks(this, permutations[r, :])
            if length(rsarks.spatialWindowed) == 0
                permResults[r, f, 1] = maximum(rsarks.windowed[fmask])
                permResults[r, f, 2] = -Inf32
            else
                permResults[r, f, 1] = -Inf32
                permResults[r, f, 2] = maximum(rsarks.spatialWindowed[fmask])
            end
        end
    end
    reset!(this, hw0, sl0, true)
    return permResults
end

function thresholdsFromPermutations(
        permDists::Array{Float32,3}, nSigma::Float64)::Array{Float64,2}
    thresholds = Array{Float64,2}(undef, size(permDists)[2], 2)
    for f in 1:size(permDists)[2]
        if any(permDists[:, f, 1] .> -Inf32)
            thresholds[f, 1] = mean(permDists[:, f, 1]) +
                               nSigma * std(permDists[:, f, 1])
            thresholds[f, 2] = -Inf
        else
            thresholds[f, 1] = -Inf
            thresholds[f, 2] = mean(permDists[:, f, 2]) +
                               nSigma * std(permDists[:, f, 2])
        end
    end
    return thresholds
end

## =============================================================================
mutable struct MicroSarks <: AbSarks
    halfWindow::Int64
    spatialLength::Int64
    subblockLength::Int64

    catScores::Array{Float32}
    blocks::Array{String}
    blockPosition::Dict{String,Int64}
    catSeq::Array{UInt8}
    bounds::Array{Int64}
    sa::Array{Int64}
    saInv::Array{Int64}
    windGini::Array{Float32}
    spatGini::Array{Float32}
    windowed::Array{Float32}
    spatialWindowed::Array{Float32}

    fm::FMIndices.FMIndex

    function MicroSarks(fasta::String, scoreFile::String,
                        halfWindow::Int64, spatialLength::Int64=0;
                        subblockLength::Int64=100,
                        shift::Int64=0)
        this = new()
        this.halfWindow = halfWindow
        this.spatialLength = spatialLength
        this.subblockLength = subblockLength
        seqs = readFasta(fasta)
        scores = DataFrame()
        if endswith(lowercase(scoreFile), r".gz(ip)?")
            scores = CSV.read(TranscodingStream(GzipDecompressor(),
                                                open(scoreFile, "r")),
                              DataFrame, delim="\t", comment="#",
                              header=false, missingstring="NA", copycols=true)
        else
            scores = CSV.read(scoreFile, DataFrame, delim="\t", comment="#",
                              header=false, missingstring="NA", copycols=true)
        end
        scores[!, 2] = Int64.(scores[:, 2]) .+ shift
        scores[!, 3] = Float64.(scores[:, 3])
        this.blocks = sort(collect(intersect(Set(scores[:, 1]), keys(seqs))))
        this.blockPosition = Dict{String,Int64}()
        for i in 1:length(this.blocks)
            this.blockPosition[this.blocks[i]] = i
        end
        toRemove = setdiff(keys(seqs), this.blocks)
        for k in toRemove
            delete!(seqs, k)
        end
        if shift != 0
            scores = scores[scores[:, 2] .> 0, :]
            for k in keys(seqs)
                kst = max(1, shift+1)
                ken = length(seqs[k]) + min(0, shift)
                seqs[k] = seqs[k][kst:ken]
            end
        end
        concatenateSeqs!(this, seqs)
        calcSuffixArray!(this)
        this.fm = FMIndices.FMIndex(this.catSeq, this.sa)
        concatenateScores!(this, scores[:, 1], scores[:, 2], scores[:, 3])
        windowGini!(this)
        window!(this)
        if this.spatialLength > 1
            spatialGini!(this)
            spatialWindow!(this)
        else
            this.spatGini = Array{Float32}(undef, 0)
            this.spatialWindowed = Array{Float32}(undef, 0)
        end
        return this
    end

    function MicroSarks(seqs::AbstractDict{<:AbstractString,<:AbstractString},
                        scores::DataFrame,
                        halfWindow::Int64, spatialLength::Int64=0;
                        subblockLength::Int64=100,
                        shift::Int64=0)
        this = new()
        this.halfWindow = halfWindow
        this.spatialLength = spatialLength
        this.subblockLength = subblockLength
        this.blocks = String.(sort(collect(intersect(Set(scores[:, 1]), keys(seqs)))))
        this.blockPosition = Dict{String,Int64}()
        for i in 1:length(this.blocks)
            this.blockPosition[this.blocks[i]] = i
        end
        toRemove = setdiff(keys(seqs), this.blocks)
        for k in toRemove
            delete!(seqs, k)
        end
        if shift != 0
            scores = deepcopy(scores)
            scores[:, 2] .+= shift
            scores = scores[scores[:, 2] .> 0, :]
            for k in keys(seqs)
                kst = max(1, shift+1)
                ken = length(seqs[k]) + max(0, shift)
                seqs[k] = seqs[k][kst:ken]
            end
        end
        concatenateSeqs!(this, seqs)
        calcSuffixArray!(this)
        this.fm = FMIndices.FMIndex(this.catSeq, this.sa)
        concatenateScores!(this, scores[:, 1], scores[:, 2], scores[:, 3])
        windowGini!(this)
        window!(this)
        if this.spatialLength > 1
            spatialGini!(this)
            spatialWindow!(this)
        else
            this.spatGini = Array{Float32}(undef, 0)
            this.spatialWindowed = Array{Float32}(undef, 0)
        end
        return this
    end

    function MicroSarks(sarks::MicroSarks)
        this = new()
        this.halfWindow = sarks.halfWindow
        this.spatialLength = sarks.spatialLength
        this.subblockLength = sarks.subblockLength
        this.blocks = sarks.blocks
        this.blockPosition = sarks.blockPosition
        this.bounds = sarks.bounds
        this.catSeq = sarks.catSeq
        this.catScores = sarks.catScores
        this.sa = sarks.sa
        this.saInv = sarks.saInv
        this.fm = sarks.fm
        this.windGini = sarks.windGini
        this.spatGini = sarks.spatGini
        this.windowed = sarks.windowed
        this.spatialWindowed = sarks.spatialWindowed
        return this
    end
end

## -----------------------------------------------------------------------------
function concatenateScores!(this::MicroSarks,
                            seq::AbstractArray{<:AbstractString},
                            pos::Array{Int64},
                            score::Array{Float64}
                            )::Nothing
    this.catScores = zeros(Float32, length(this.catSeq))
    for r in 1:length(seq)
        rb = this.blockPosition[String(seq[r])]
        rs = this.bounds[rb] + pos[r] - 1
        if rs < this.bounds[rb+1]
            this.catScores[ this.saInv[rs] ] = score[r]
        end
    end
    return nothing
end

## -----------------------------------------------------------------------------
function fillWindGini!(windGini::Array{Float32},
                       halfWindow::Int64, sa::Array{Int64},
                       bounds::Array{Int64},
                       block::Array{Int64}, bCounts::Dict{Int64,Int64},
                       subblockLength::Int64, runningGini::Float64)::Nothing
    total = (2 * halfWindow) + 1
    total2 = convert(Float64, total * total)
    dGini = 0
    oldb = -1
    newb = -1
    saLen = length(sa)
    nBlocks = length(bounds) - 1
    @inbounds for i = (halfWindow+2):(saLen-halfWindow)
        oldb = block[i - 1 - halfWindow]
        oldsb = 1 + ((sa[i - 1 - halfWindow] - bounds[oldb]) ÷ subblockLength)
        oldkey = (oldsb * nBlocks) + oldb
        newb = block[i + halfWindow]
        newsb = 1 + ((sa[i + halfWindow] - bounds[newb]) ÷ subblockLength)
        newkey = (newsb * nBlocks) + newb
        dGini = -bCounts[oldkey] * (total - bCounts[oldkey])
        if !(newkey in keys(bCounts))
            bCounts[newkey] = 0
        end
        if bCounts[newkey] > 0
            dGini -= bCounts[newkey] * (total - bCounts[newkey])
        end
        bCounts[oldkey] -= 1
        bCounts[newkey] += 1
        if bCounts[oldkey] > 0
            dGini += bCounts[oldkey] * (total - bCounts[oldkey])
        end
        dGini += bCounts[newkey] * (total - bCounts[newkey])
        runningGini += (convert(Float64, dGini) / total2)
        windGini[i] = convert(Float32, runningGini)
    end
end

function blockCounts(this::MicroSarks, i::Int64, block::Array{Int64})::Dict{Int64,Int64}
    out = Dict{Int64,Int64}()
    nBlocks = length(this.bounds) - 1
    for j = -this.halfWindow:this.halfWindow
        jb = block[i+j]
        jsb = 1 + ((this.sa[i+j] - this.bounds[jb]) ÷ this.subblockLength)
        jkey = (jsb * nBlocks) + jb
        if !(jkey in keys(out))
            out[jkey] = 0
        end
        out[jkey] += 1
    end
    return out
end

function windowGini!(this::MicroSarks)::Nothing
    block = sourceBlock(this)
    this.windGini = Array{Float32}(undef, length(this.sa))
    bCounts = blockCounts(this, this.halfWindow+1, block)
    runningGini = giniImpurity(collect(values(bCounts)))
    this.windGini[this.halfWindow+1] = convert(Float32, runningGini)
    fillWindGini!(this.windGini, this.halfWindow, this.sa, this.bounds,
                  block, bCounts, this.subblockLength, runningGini)
    for i = 1:this.halfWindow
        this.windGini[i] = this.windGini[this.halfWindow+1]
        this.windGini[length(this.sa) + 1 - i] =
                this.windGini[length(this.sa) - this.halfWindow]
    end
    return nothing
end

## -----------------------------------------------------------------------------
function fillWindowed!(windowed::Array{Float32},
                       halfWindow::Int64, saLen::Int64,
                       scores::Array{Float32},
                       runningSum::Float64, windSize=Float64)::Float64
    @inbounds for i = (halfWindow+2):(saLen-halfWindow)
        runningSum -= scores[i-(halfWindow+1)]
        runningSum += scores[i+halfWindow]
        windowed[i] = convert(Float32, runningSum / windSize)
    end
    return runningSum
end

function window!(this::MicroSarks)::Nothing
    windSize = convert(Float64, (2 * this.halfWindow) + 1)
    srcBlk = sourceBlock(this)
    this.windowed = Array{Float32}(undef, length(this.sa))
    runningSum = 0.0
    for i = 1:((2*this.halfWindow)+1)
        runningSum += this.catScores[i]
    end
    this.windowed[this.halfWindow + 1] = convert(Float32, runningSum / windSize)
    runningSum = fillWindowed!(this.windowed, this.halfWindow,
                               length(this.sa), this.catScores,
                               runningSum, windSize)
    for i = 1:this.halfWindow
        this.windowed[i] = this.windowed[this.halfWindow+1]
        this.windowed[length(this.sa) + 1 - i] =
                convert(Float32, runningSum / windSize)
    end
    return nothing
end

## =============================================================================
function rescore(sarks::Sarks,
                 newscores::AbstractDict{<:AbstractString,Float64}
                 )::Sarks
    this = Sarks(sarks)
    this.scores = [newscores[b] for b in this.blocks]
    this.origScores = [newscores[b] for b in this.blocks]
    if size(this.charFreqs, 1) > 0
        deltaScores = (this.origScores .- mean(this.origScores))
        hatScores = this.charFreqs * (this.compositionPreBeta * deltaScores)
        this.scores = deltaScores - hatScores
    else
        this.scores = this.origScores
    end
    Sarkses.window!(this)
    if this.spatialLength > 1
        Sarkses.spatialWindow!(this)
    else
        this.spatialWindowed = Array{Float32}(undef, 0)
    end
    return this
end

function multisarks(this::Sarks,
                    scores::Array{OrderedDict{<:AbstractString,Float64}}
                    )::Array{Sarks}
    out = Array{Sarks}(undef, length(scores))
    Threads.@threads for i in 1:length(scores)
        out[i] = rescore(this, scores[i])
    end
    return out
end

## 210530 DW: not clear centerDicts still needed...
function centerDicts(x::Array{OrderedDict{<:AbstractString,Float64}}
                     )::Array{OrderedDict{String,Float64}}
    out = Array{OrderedDict{String,Float64}}(undef, length(x))
    for i in 1:length(x)
        outi = OrderedDict{String,Float64}()
        xi = x[i]
        ximu = mean(values(xi))
        for k in keys(xi)
            outi[String(k)] = xi[k] - ximu
        end
        out[i] = outi
    end
    return out
end

function multiWindCovar(sarkses::Array{Sarks},
                        minGini::Float64=-Inf
                        )::Matrix{Float64}
    if minGini > 0.0
        gm = giniMask(sarkses[1], minGini, -Inf)
    end
    out = zeros(length(sarkses), length(sarkses))
    for i in 1:length(sarkses)
        for j in i:length(sarkses)
            if minGini > 0.0
                out[i, j] = sum(
                    (sarkses[i].windowed[gm] .- mean(sarkses[i].windowed[gm])) .*
                    (sarkses[j].windowed[gm] .- mean(sarkses[j].windowed[gm]))
                )
            else
                out[i, j] = sum(
                    (sarkses[i].windowed .- mean(sarkses[i].windowed)) .*
                    (sarkses[j].windowed .- mean(sarkses[j].windowed))
                )
            end
            if j > i
                out[j, i] = out[i, j]
            end
        end
    end
    return out
end

function eigenScores(sarkses::Array{Sarks},
                     minGini::Float64=-Inf;
                     retainFirst::Bool=false
                     )::Array{OrderedDict{String,Float64}}
    scovar = multiWindCovar(sarkses, minGini)
    ev = zeros(size(scovar, 1), size(scovar, 2))
    if retainFirst
        ev[1, end] = 1.0
        ev[2:end, 1:(end-1)] = eigen(scovar[2:end, 2:end]).vectors
    else
        ev = eigen(scovar).vectors
    end
    out = Array{OrderedDict{String,Float64}}(undef, length(sarkses))
    for i in 1:length(out)
        evi = ev[:, length(sarkses)-i+1]
        out[i] = OrderedDict([
            k => sum([evi[j] * sarkses[j].scores[b] for j in 1:length(evi)])
            for (b, k) in enumerate(sarkses[1].blocks)
        ]);
    end
    return out
end

## =============================================================================
function spatialMeansInner(bounds::Array{Int64},
                           saInv::Array{Int64},
                           windowed::Array{Float32},
                           filt::AbstractArray{Bool})::Array{Float64}
    blockMeans = zeros(length(bounds)-1)
    activeBlock = 1
    blockLength = 0.0
    @inbounds for s in 1:length(saInv)
        if s >= bounds[activeBlock+1]
            blockMeans[activeBlock] /= blockLength
            activeBlock += 1
            blockLength = 0.0
        end
        if filt[saInv[s]]
            blockMeans[activeBlock] += windowed[saInv[s]]
            blockLength += 1.0
        end
    end
    blockMeans[activeBlock] /= blockLength
    return blockMeans
end

function spatialMeans(this::Sarks,
                      minGini::Float64=-Inf)::OrderedDict{String,Float64}
    if minGini >= 1.0
        minGini = 1.0 - (1.0 - median(this.windGini)) * minGini
    end
    blockMeans = spatialMeansInner(this.bounds, this.saInv, this.windowed,
                                   this.windGini .>= minGini)
    out = OrderedDict{String,Float64}()
    for (b, block) in enumerate(this.blocks)
        out[block] = blockMeans[b]
    end
    return out
end

mutable struct SarksScorer
    sarks::Sarks
    suffixLen::Int64
    spatial::Bool
    function SarksScorer(sarks::Sarks,
                         suffixLen::Int64;
                         spatial::Bool=false)
        this = new()
        this.sarks = sarks
        this.suffixLen = suffixLen
        this.spatial = spatial
        return this
    end
end

function score(this::SarksScorer, seq::Array{UInt8},
               minGini::Float64=-Inf; naOmit=true)::Array{Float64}
    if minGini >= 1.0
        minGini = 1.0 - (1.0 - median(this.sarks.windGini)) * minGini
    end
    scores = zeros(min(length(seq), length(seq)-this.suffixLen+1))
    Threads.@threads for s = 1:length(scores)
        suffixEnd = min(s+this.suffixLen-1, length(seq))
        suffix = seq[s:suffixEnd]
        limits = Sarkses.findKmer(this.sarks, suffix)
        lower = min(limits[1], limits[2]-1)
        upper = min(max(limits[1], limits[2]-1), length(this.sarks.sa))
        num = 0.0
        den = 0.0
        for i = lower:upper
            if this.sarks.windGini[i] >= minGini
                if this.spatial
                    num += this.sarks.spatialWindowed[i]
                else
                    num += this.sarks.windowed[i]
                end
                den += 1.0
            end
        end
        scores[s] = num / den
    end
    if naOmit
        return scores[.!isnan.(scores)]
    else
        return scores
    end
end

function scoreWithGini(this::SarksScorer,
                       seq::Array{UInt8})::Tuple{Array{Float64},Array{Float64}}
    scores = zeros(min(length(seq), length(seq)-this.suffixLen+1))
    ginis = zeros(length(scores))
    Threads.@threads for s = 1:length(scores)
        suffixEnd = min(s+this.suffixLen-1, length(seq))
        suffix = seq[s:suffixEnd]
        limits = Sarkses.findKmer(this.sarks, suffix)
        lower = min(limits[1], limits[2]-1)
        upper = min(max(limits[1], limits[2]-1), length(this.sarks.sa))
        sScore = 0.0
        sGini = 0.0
        for i = lower:upper
            if this.spatial
                sScore += this.sarks.spatialWindowed[i]
                sGini += this.sarks.spatGini[i]
            else
                sScore += this.sarks.windowed[i]
                sGini += this.sarks.windGini[i]
            end
        end
        scores[s] = sScore / (1+upper-lower)
        ginis[s] = sGini / (1+upper-lower)
    end
    return (scores, ginis)
end

function sarksPositionalScore(this::Sarks,
                              seqs::OrderedDict{<:AbstractString,Array{UInt8,1}};
                              k::Int64=25)::DataFrame
    spat = max(this.spatialLength, 1)
    scorer = SarksScorer(this, k, spatial=(spat > 1))
    nSeqBins = OrderedDict{String,Int64}()
    for sk in keys(seqs)
        nSeqBins[String(sk)] = Int64(floor((length(seqs[sk])+1-k) / spat))
    end
    nTotalBins = sum(values(nSeqBins))
    id = Array{String}(undef, nTotalBins)
    bin = Array{Int64}(undef, nTotalBins)
    score = Array{Float64}(undef, nTotalBins)
    gini = Array{Float64}(undef, nTotalBins)
    offset = 0
    offsets = OrderedDict([sk => 0 for sk in keys(seqs)])
    for sk in keys(seqs)
        offsets[sk] = offset
        offset += nSeqBins[String(sk)]
    end
    seqKeys = collect(values(keys(seqs)))
    Threads.@threads for sk in seqKeys
        skOffset = offsets[sk]
        skScores = scoreWithGini(scorer, seqs[sk])
        for i in 1:nSeqBins[String(sk)]
            iOffset = i + skOffset
            id[iOffset] = String(sk)
            bin[iOffset] = i
            score[iOffset] = skScores[1][1 + (i-1)*spat]
            gini[iOffset] = skScores[2][1 + (i-1)*spat]
        end
    end
    out = DataFrame(id=id, bin=bin, score=score, gini=gini)
    return out
end

function sarksPositionalScore(this::Sarks,
                              seqs::OrderedDict{<:AbstractString,<:AbstractString};
                              k::Int64=25)::DataFrame
    return sarksPositionalScore(
        this,
        OrderedDict([k => Array{UInt8}(v*"\$") for (k, v) in seqs]),
        k = k
    )
end

function sarksPositionalScoreMat(this::Sarks,
                                 seqs::OrderedDict{<:AbstractString,<:AbstractString};
                                 k::Int64=25,
                                 side::Symbol=:left,
                                 fill::Symbol=:global
                                 )::Matrix
    out = sarksPositionalScore(this, seqs, k=k)
    if fill == :individual
        return Matrix(c1idx(padAndUnstack(out, side=side),
                            collect(keys(seqs)))[:, 2:end])
    else
        return Matrix(c1idx(padAndUnstack(out,
                                          side = side,
                                          fill = mean(this.scores)),
                            collect(keys(seqs)))[:, 2:end])
    end
end

function predictSarks(sarks::Sarks,
                      seqs::OrderedDict{<:AbstractString,<:AbstractString};
                      k::Int64=25,
                      side::Symbol=:left,
                      fill::Symbol=:global
                      )::Tuple{Matrix{Float32},Matrix{Float32}}
    posScores = sarksPositionalScore(
        sarks,
        OrderedDict([id => Array{UInt8}(seqs[id]*"\$") for id in keys(seqs)]),
        k = k
    )
    if fill == :individual
        posScoreMat = Matrix{Float32}(c1idx(
            padAndUnstack(posScores, side=side),
            collect(keys(seqs))
        )[:, 2:end])
        posGiniMat = Matrix{Float32}(c1idx(
            padAndUnstack(posScores, :id, :bin, :gini, side=side),
            collect(keys(seqs))
        )[:, 2:end])
        return (posScoreMat, posGiniMat)
    else
        posScoreMat = Matrix{Float32}(c1idx(
            padAndUnstack(posScores,
                          side = side,
                          fill = mean(sarks.scores)),
            collect(keys(seqs))
        )[:, 2:end])
        posGiniMat = Matrix{Float32}(c1idx(
            padAndUnstack(posScores, :id, :bin, :gini,
                          side = side,
                          fill = Float64(mean(sarks.windGini))),
            collect(keys(seqs))
        )[:, 2:end])
        return (posScoreMat, posGiniMat)
   end
end

## =============================================================================
mutable struct SarksBoost
    learningRate::Float64
    minGini::Float64
    sarkses::Array{Sarks}
    coef::Array{Float64}
    intercept::Float64
end

function simpleLinReg(x::Vector{Float64}, y::Vector{Float64})::Tuple{Float64, Float64}
    ymu = mean(y)
    dx = x .- mean(x)
    dy = y .- ymu
    slope = sum(dx .* dy) / sum(dx .* dx)
    intercept = ymu - slope * mean(x)
    return (intercept, slope)
end

function boostSarks(iterations::Int64,
                    seqs::AbstractDict{<:AbstractString,<:AbstractString},
                    scoreMap::AbstractDict{<:AbstractString,Float64},
                    halfWindow::Int64,
                    minGini::Float64=-Inf;
                    shrinkage::Float64=0.0,
                    spatialLength::Int64=0,
                    adjustForCharFreqs::Bool=false)::SarksBoost
    lr = 1.0 - shrinkage
    sarks1 = Sarks(seqs, scoreMap, halfWindow, spatialLength,
                   adjustForCharFreqs = adjustForCharFreqs)
    ordScores = [scoreMap[b] for b in sarks1.blocks]
    ordSeqs = [String(seqs[b]) for b in sarks1.blocks]
    resids0 = ordScores .- mean(ordScores)
    spatMeans1 = collect(values(spatialMeans(sarks1, minGini)))
    sarks1Mod = simpleLinReg(spatMeans1, resids0)
    out = SarksBoost(lr, minGini,
                     [sarks1], [lr*sarks1Mod[2]],
                     mean(ordScores) + lr*sarks1Mod[1])
    if iterations > 1
        pseudoResids = resids0 - lr * (sarks1Mod[1] .+ sarks1Mod[2]*spatMeans1)
        for iter = 2:iterations
            sarksIter = rescore(
                sarks1,
                OrderedDict([b => pseudoResids[i]
                             for (i, b) in enumerate(sarks1.blocks)])
            )
            spatMeansIter = collect(values(spatialMeans(sarksIter, minGini)))
            iterMod = simpleLinReg(spatMeansIter, pseudoResids)
            push!(out.sarkses, sarksIter)
            push!(out.coef, lr*iterMod[2])
            out.intercept += lr*iterMod[1]
            pseudoResids -= lr * (iterMod[1] .+ iterMod[2]*spatMeansIter)
        end
    end
    return out
end

function sarksPositionalScoreMat(this::SarksBoost,
                                 seqs::OrderedDict{<:AbstractString,<:AbstractString};
                                 k::Int64)::Matrix
    out = this.coef[1] * sarksPositionalScoreMat(this.sarkses[1], seqs, k=k);
    for i = 2:length(this.sarkses)
        out += this.coef[i] * sarksPositionalScoreMat(this.sarkses[i], seqs, k=k);
    end
    out .+= this.intercept
    return out
end

end # module Sarkses