text:
using RedPitayaDAQServer
using PyPlot

# obtain the URL of the RedPitaya
include("config.jl")

rp = RedPitayaCluster([URLs[1]])

dec = 64
modulus = 4800
base_frequency = 125000000
periods_per_step = 5
samples_per_period = div(modulus, dec)
periods_per_frame = 50 # about 0.5 s frame length
frame_period = dec*samples_per_period*periods_per_frame / base_frequency
slow_dac_periods_per_frame = div(50, periods_per_step)

stopADC(rp)
masterTrigger(rp, false)
clearSequence(rp)

decimation(rp, dec)
samplesPerPeriod(rp, samples_per_period)
periodsPerFrame(rp, periods_per_frame)
passPDMToFastDAC(master(rp), true)

modeDAC(rp, "STANDARD")
frequencyDAC(rp,1,1, base_frequency / modulus)

freq = frequencyDAC(rp,1,1)
println(" frequency = $(freq)")
signalTypeDAC(rp, 1 , "SINE")
amplitudeDAC(rp, 1, 1, 0.1)
phaseDAC(rp, 1, 1, 0.0 ) # Phase has to be given in between 0 and 1

ramWriterMode(rp, "TRIGGERED")
triggerMode(rp, "EXTERNAL")

# Sequence
slowDACStepsPerFrame(rp, slow_dac_periods_per_frame)
numSlowDACChan(master(rp), 1)
sleep(0.5)

fig = figure(1)
clf()

masterTrigger(rp, false)

# Lookup Sequence
clearSequence(rp)
amplitudeDAC(rp, 1, 1, 0.1)
lut = collect(range(0,0.7,length=slow_dac_periods_per_frame))
seq = ArbitrarySequence(lut, nothing, slow_dac_periods_per_frame, 2, 0, 0)
appendSequence(master(rp), seq)
success = prepareSequence(master(rp))
startADC(rp)
masterTrigger(rp, true)
sleep(0.1)

uCurrentFrame = readFrames(rp, 0, 3)
subplot(2,2, 1)
plot(vec(uCurrentFrame[:, 1, :, :]))
title("Arbitrary")
stopADC(rp)
masterTrigger(rp, false)

sleep(0.1)
# Constant Sequence
clearSequence(rp)
amplitudeDAC(rp, 1, 1, 0.1) # Amplitude is set to zero after a sequence
lut = [0.2]
seq = ConstantSequence(lut, nothing, slow_dac_periods_per_frame, 2, 0, 0)
appendSequence(master(rp), seq)
success = prepareSequence(master(rp))
startADC(rp)
masterTrigger(rp, true)
sleep(0.1)

uCurrentFrame = readFrames(rp, 0, 3)
subplot(2,2, 2)
plot(vec(uCurrentFrame[:, 1, :, :]))
title("Constant")
stopADC(rp)
masterTrigger(rp, false)

sleep(0.1)
# Pause Sequence
clearSequence(rp)
amplitudeDAC(rp, 1, 1, 0.1)
seq = PauseSequence(nothing, slow_dac_periods_per_frame, 2)
appendSequence(master(rp), seq)
success = prepareSequence(master(rp))
startADC(rp)
masterTrigger(rp, true)
sleep(0.1)

uCurrentFrame = readFrames(rp, 0, 3)
subplot(2,2, 3)
plot(vec(uCurrentFrame[:, 1, :, :]))
title("Pause")
stopADC(rp)
masterTrigger(rp, false)

sleep(0.1)
# Range Sequence
clearSequence(rp)
amplitudeDAC(rp, 1, 1, 0.1)
lut = [-1.0, 2.0/slow_dac_periods_per_frame]
seq = RangeSequence(lut, nothing, slow_dac_periods_per_frame, 2, 0, 0)
appendSequence(master(rp), seq)
success = prepareSequence(master(rp))
startADC(rp)
masterTrigger(rp, true)
sleep(0.1)

uCurrentFrame = readFrames(rp, 0, 3)
subplot(2,2, 4)
plot(vec(uCurrentFrame[:, 1, :, :]))
title("Range")
stopADC(rp)
masterTrigger(rp, false)
clearSequence(rp)

fig

max_stars_repo_name: tknopp/RedPitayaDAQServer

max_issues_repo_name: tknopp/RedPitayaDAQServer

max_forks_repo_name: tknopp/RedPitayaDAQServer


text:
# Tools to analyze parts of tensor expressions and extract information

# analyze tensor expressions and extract information
function decomposetensor(ex::Expr)
    istensor(ex) || throw(ArgumentError("not a valid tensor: $ex"))
    if ex.head == :ref || ex.head == :typed_hcat
        if length(ex.args) == 1
            return ex.args[1], Any[], Any[]
        elseif isa(ex.args[2], Expr) && ex.args[2].head == :parameters
            return ex.args[1], ex.args[3:end], ex.args[2].args
        else
            return ex.args[1], ex.args[2:end], Any[]
        end
    else #if ex.head == :typed_vcat
        if isa(ex.args[2], Expr) && (ex.args[2].head == :row || ex.args[2].head == :tuple)
            leftind = ex.args[2].args
        else
            leftind = ex.args[2:2]
        end
        if isa(ex.args[3], Expr) && (ex.args[3].head == :row || ex.args[3].head == :tuple)
            rightind = ex.args[3].args
        else
            rightind = ex.args[3:3]
        end
        return ex.args[1], leftind, rightind
    end
end

gettensorobject(ex) = decomposetensor(ex)[1]
getleftindices(ex) = decomposetensor(ex)[2]
getrightindices(ex) = decomposetensor(ex)[3]

function decomposegeneraltensor(ex)
    if istensor(ex)
        object, leftind, rightind = decomposetensor(ex)
        return (object, leftind, rightind, true, false)
    elseif isa(ex, Expr) && ex.head == :call && ex.args[1] == :+ && length(ex.args) == 2 # unary plus: pass on
        return decomposegeneraltensor(ex.args[2])
    elseif isa(ex, Expr) && ex.head == :call && ex.args[1] == :- && length(ex.args) == 2 # unary minus: flip scalar factor
        (object, leftind, rightind, α, conj) = decomposegeneraltensor(ex.args[2])
        return (object, leftind, rightind, Expr(:call, :-, α), conj)
    elseif isa(ex, Expr) && ex.head == :call && ex.args[1] == :conj && length(ex.args) == 2 # conjugation: flip conjugation flag and conjugate scalar factor
        (object, leftind, rightind, α, conj) = decomposegeneraltensor(ex.args[2])
        return (object, leftind, rightind, Expr(:call, :conj, α), !conj)
    elseif ex.head == :call && ex.args[1] == :* && length(ex.args) == 3 # scalar multiplication: multiply scalar factors
        if isscalarexpr(ex.args[2]) && isgeneraltensor(ex.args[3])
            (object, leftind, rightind, α, conj) = decomposegeneraltensor(ex.args[3])
            return (object, leftind, rightind, Expr(:call, :*, ex.args[2], α), conj)
        elseif isscalarexpr(ex.args[3]) && isgeneraltensor(ex.args[2])
            (object, leftind, rightind, α, conj) = decomposegeneraltensor(ex.args[2])
            return (object, leftind, rightind, Expr(:call, :*, α, ex.args[3]), conj)
        end
    elseif ex.head == :call && ex.args[1] == :/ && length(ex.args) == 3 # scalar multiplication: muliply scalar factors
        if isscalarexpr(ex.args[3]) && isgeneraltensor(ex.args[2])
            (object, leftind, rightind, α, conj) = decomposegeneraltensor(ex.args[2])
            return (object, leftind, rightind, Expr(:call, :/, α, ex.args[3]), conj)
        end
    elseif ex.head == :call && ex.args[1] == :\ && length(ex.args) == 3 # scalar multiplication: muliply scalar factors
        if isscalarexpr(ex.args[2]) && isgeneraltensor(ex.args[3])
            (object, leftind, rightind, α, conj) = decomposegeneraltensor(ex.args[3])
            return (object, leftind, rightind, Expr(:call, :\, ex.args[2], α), conj)
        end
    end
    throw(ArgumentError("not a valid generalized tensor expression $ex"))
end

# get left and right hand side from an expression or definition
function getlhs(ex::Expr)
    if ex.head in (:(=), :(+=), :(-=), :(:=), :(≔))
        return ex.args[1]
    else
        throw(ArgumentError("invalid assignment or definition $ex"))
    end
end
function getrhs(ex::Expr)
    if ex.head in (:(=), :(+=), :(-=), :(:=), :(≔))
        return ex.args[2]
    else
        throw(ArgumentError("invalid assignment or definition $ex"))
    end
end

# get all the tensors and tensor objects in a tensor expression (not a definition or assignment)
function gettensors(ex)
    if istensor(ex)
        Any[ex]
    elseif istensorexpr(ex)
        list = Any[]
        for e in ex.args
            append!(list, gettensors(e))
        end
        return list
    else
        return Any[]
    end
end

gettensorobjects(ex) = gettensorobject.(gettensors(ex))

# get all the existing tensor objects which are inputs (i.e. appear in the rhs of assignments and definitions)
function getinputtensorobjects(ex)
    list = Any[]
    if istensorexpr(ex)
        append!(list, gettensorobjects(ex))
    elseif isdefinition(ex)
        append!(list, gettensorobjects(getrhs(ex)))
    elseif isassignment(ex)
        if ex.head == :(+=) || ex.head == :(-=)
            lhs = getlhs(ex)
            if istensor(lhs)
                push!(list, gettensorobject(lhs))
            end
        end
        append!(list, gettensorobjects(getrhs(ex)))
    elseif isa(ex, Expr) && ex.head == :block
        for i = 1:length(ex.args)
            list2 = getinputtensorobjects(ex.args[i])
            for j = 1:i-1
                # if objects have previously been defined or assigned to, they are not inputs
                list2 = setdiff(list2, getnewtensorobjects(ex.args[j]))
                list2 = setdiff(list2, getoutputtensorobjects(ex.args[j]))
            end
            append!(list, list2)
        end
    elseif isa(ex, Expr) && ex.head in (:for, :while)
        append!(list, getinputtensorobjects(ex.args[2]))
    elseif isa(ex, Expr) && ex.head == :call && ex.args[1] == :scalar
        append!(list, gettensorobjects(ex.args[2]))
    end
    return unique!(list)
end
# get all the existing tensor objects which are outputs (i.e. appear in the lhs of assignments)
function getoutputtensorobjects(ex)
    list = Any[]
    if isassignment(ex)
        lhs = getlhs(ex)
        if istensor(lhs)
            push!(list, gettensorobject(lhs))
        end
    elseif isa(ex, Expr) && ex.head == :block
        for i = 1:length(ex.args)
            list2 = getoutputtensorobjects(ex.args[i])
            for j = 1:i-1
                # if objects have previously been defined, they are not existing outputs
                list2 = setdiff(list2, getnewtensorobjects(ex.args[j]))
            end
            append!(list, list2)
        end
    elseif isa(ex, Expr) && ex.head in (:for, :while)
        append!(list, getoutputtensorobjects(ex.args[2]))
    end
    return unique!(list)
end
# get all the existing tensor objects which are newly created (i.e. appear in the lhs of definition)
function getnewtensorobjects(ex)
    list = Any[]
    if isdefinition(ex)
        lhs = getlhs(ex)
        if istensor(lhs)
            push!(list, gettensorobject(lhs))
        end
    elseif isa(ex, Expr) && ex.head == :block
        for e in ex.args
            append!(list, getnewtensorobjects(e))
        end
    elseif isa(ex, Expr) && ex.head in (:for, :while)
        append!(list, getnewtensorobjects(ex.args[2]))
    end
    return unique!(list)
end

# for any tensor expression, get the list of uncontracted indices that would remain after evaluating that expression
function getindices(ex::Expr)
    if istensor(ex)
        _,leftind,rightind = decomposetensor(ex)
        return unique2(vcat(leftind, rightind))
    elseif ex.head == :call && (ex.args[1] == :+ || ex.args[1] == :-)
        return getindices(ex.args[2]) # getindices on any of the args[2:end] should yield the same result
    elseif ex.head == :call && ex.args[1] == :*
        indices = getindices(ex.args[2])
        for k = 3:length(ex.args)
            append!(indices, getindices(ex.args[k]))
        end
        return unique2(indices)
    elseif ex.head == :call && ex.args[1] == :/
        return getindices(ex.args[2])
    elseif ex.head == :call && ex.args[1] == :\
        return getindices(ex.args[3])
    elseif ex.head == :call && length(ex.args) == 2
        return getindices(ex.args[2])
    else
        return Any[]
    end
end
getindices(ex) = Any[]

# get all unique indices appearing in a tensor expression
function getallindices(ex::Expr)
    if istensor(ex)
        _,leftind,rightind = decomposetensor(ex)
        return unique!(vcat(leftind, rightind))
    elseif isassignment(ex) || isdefinition(ex)
        return getallindices(getrhs(ex))
    else
        return unique!(mapreduce(getallindices, vcat, ex.args))
    # elseif istensorexpr(ex)
    #     return unique!(mapreduce(getallindices, vcat, ex.args))
    # else
    #     return Any[]
    end
end
getallindices(ex) = Any[]

# # auxiliary routine
# function unique2(itr)
#     out = collect(itr)
#     i = 1
#     while i < length(out)
#         inext = findnext(isequal(out[i]), out, i+1)
#         if inext == nothing
#             i += 1
#             continue
#         end
#         while inext !== nothing
#             deleteat!(out, inext)
#             inext = findnext(isequal(out[i]), out, i+1)
#         end
#         deleteat!(out, i)
#     end
#     out
# end

ext: jl


text:
export LSRK54CarpenterKennedy, LSRK144NiegemannDiehlBusch, LSRKEulerMethod

"""
    LowStorageRungeKutta2N <: DistributedODEAlgorithm

A class of low-storage Runge-Kutta algorithms, which use only one additional
copy of the state vector ``u`` (often referred to as ``2N`` schemes).

The available implementations are:
 - [`LSRKEulerMethod`](@ref)
 - [`LSRK54CarpenterKennedy`](@ref)
 - [`LSRK144NiegemannDiehlBusch`](@ref)
"""
abstract type LowStorageRungeKutta2N <: DistributedODEAlgorithm end

"""
    LowStorageRungeKutta2NTableau

Storage for the tableau of a [`LowStorageRungeKutta2N`](@ref) algorithm.
"""
struct LowStorageRungeKutta2NTableau{Nstages, RT}
    "low storage RK coefficient vector A (rhs scaling)"
    A::NTuple{Nstages, RT}
    "low storage RK coefficient vector B (rhs add in scaling)"
    B::NTuple{Nstages, RT}
    "low storage RK coefficient vector C (time scaling)"
    C::NTuple{Nstages, RT}
end

struct LowStorageRungeKutta2NIncCache{Nstages, RT, A}
    tableau::LowStorageRungeKutta2NTableau{Nstages, RT}
    du::A
end

function cache(prob::DiffEqBase.ODEProblem, alg::LowStorageRungeKutta2N; kwargs...)
    # @assert prob.problem_type isa DiffEqBase.IncrementingODEProblem ||
    #     prob.f isa DiffEqBase.IncrementingODEFunction
    du = zero(prob.u0)
    return LowStorageRungeKutta2NIncCache(tableau(alg, eltype(du)), du)
end

nstages(::LowStorageRungeKutta2NIncCache{N}) where {N} = N

function step_u!(int, cache::LowStorageRungeKutta2NIncCache)
    tab = cache.tableau
    du = cache.du

    u = int.u
    p = int.prob.p
    t = int.t
    dt = int.dt

    for stage in 1:nstages(cache)
        #  du .= f(u, p, t + tab.C[stage]*dt) .+ tab.A[stage] .* du
        stage_time = t + tab.C[stage]*dt
        int.prob.f(du, u, p, stage_time, 1, tab.A[stage])
        u .+= (dt*tab.B[stage]) .* du
    end
end

# for Multirate
function init_inner(prob, outercache::LowStorageRungeKutta2NIncCache, dt)
    OffsetODEFunction(prob.f.f1, zero(dt), one(dt), zero(dt), outercache.du)
end
function update_inner!(innerinteg, outercache::LowStorageRungeKutta2NIncCache,
        f_slow, u, p, t, dt, stage)

    f_offset = innerinteg.prob.f
    tab = outercache.tableau
    N = nstages(outercache)

    τ0 = t+tab.C[stage]*dt
    τ1 = stage == N ? t+dt : t+tab.C[stage+1]*dt
    f_offset.α = τ0
    innerinteg.t = zero(τ0)
    innerinteg.tstop = τ1-τ0

    #  du .= f(u, p, t + tab.C[stage]*dt) .+ tab.A[stage] .* du
    f_slow(f_offset.x, u, p, τ0, 1, tab.A[stage])

    C0 = tab.C[stage]
    C1 = stage == N ? one(tab.C[stage]) : tab.C[stage+1]
    f_offset.γ = tab.B[stage] / (C1-C0)
end

"""
    LSRK54CarpenterKennedy()

The 4th-order, 5-stage [`LowStorageRungeKutta2N`])(ref) scheme from Solution
3 of [CK1994](@cite).
"""
struct LSRK54CarpenterKennedy <: LowStorageRungeKutta2N end

function tableau(::LSRK54CarpenterKennedy, RT)
    RKA = (
        RT(0),
        RT(-567301805773 // 1357537059087),
        RT(-2404267990393 // 2016746695238),
        RT(-3550918686646 // 2091501179385),
        RT(-1275806237668 // 842570457699),
    )

    RKB = (
        RT(1432997174477 // 9575080441755),
        RT(5161836677717 // 13612068292357),
        RT(1720146321549 // 2090206949498),
        RT(3134564353537 // 4481467310338),
        RT(2277821191437 // 14882151754819),
    )

    RKC = (
        RT(0),
        RT(1432997174477 // 9575080441755),
        RT(2526269341429 // 6820363962896),
        RT(2006345519317 // 3224310063776),
        RT(2802321613138 // 2924317926251),
    )

    return LowStorageRungeKutta2NTableau(RKA, RKB, RKC)
end

"""
    LSRK144NiegemannDiehlBusch()

The 4th-order, 14-stage, [`LowStorageRungeKutta2N`])(ref) scheme of
[NDB2012](@cite) with optimized stability region
"""
struct LSRK144NiegemannDiehlBusch <: LowStorageRungeKutta2N end

function tableau(::LSRK144NiegemannDiehlBusch, RT)

    RKA = (
        RT(0),
        RT(-0.7188012108672410),
        RT(-0.7785331173421570),
        RT(-0.0053282796654044),
        RT(-0.8552979934029281),
        RT(-3.9564138245774565),
        RT(-1.5780575380587385),
        RT(-2.0837094552574054),
        RT(-0.7483334182761610),
        RT(-0.7032861106563359),
        RT(0.0013917096117681),
        RT(-0.0932075369637460),
        RT(-0.9514200470875948),
        RT(-7.1151571693922548),
    )

    RKB = (
        RT(0.0367762454319673),
        RT(0.3136296607553959),
        RT(0.1531848691869027),
        RT(0.0030097086818182),
        RT(0.3326293790646110),
        RT(0.2440251405350864),
        RT(0.3718879239592277),
        RT(0.6204126221582444),
        RT(0.1524043173028741),
        RT(0.0760894927419266),
        RT(0.0077604214040978),
        RT(0.0024647284755382),
        RT(0.0780348340049386),
        RT(5.5059777270269628),
    )

    RKC = (
        RT(0),
        RT(0.0367762454319673),
        RT(0.1249685262725025),
        RT(0.2446177702277698),
        RT(0.2476149531070420),
        RT(0.2969311120382472),
        RT(0.3978149645802642),
        RT(0.5270854589440328),
        RT(0.6981269994175695),
        RT(0.8190890835352128),
        RT(0.8527059887098624),
        RT(0.8604711817462826),
        RT(0.8627060376969976),
        RT(0.8734213127600976),
    )

    LowStorageRungeKutta2NTableau(RKA, RKB, RKC)
end

"""
    LSRKEulerMethod()

An implementation of explicit Euler method using [`LowStorageRungeKutta2N`](@ref) infrastructure.
This is mainly for debugging.
"""
struct LSRKEulerMethod <: LowStorageRungeKutta2N
end
tableau(::LSRKEulerMethod, RT) =
    LowStorageRungeKutta2NTableau((RT(0),), (RT(1),), (RT(0),))


text:
module ProximalAlgorithms

const RealOrComplex{R} = Union{R, Complex{R}}
const Maybe{T} = Union{T, Nothing}

include("compat.jl")

# utilities

include("utilities/conjugate.jl")
include("utilities/fbetools.jl")
include("utilities/iterationtools.jl")
include("utilities/lbfgs.jl")

# algorithms

include("algorithms/forwardbackward.jl")
include("algorithms/zerofpr.jl")
include("algorithms/panoc.jl")
include("algorithms/douglasrachford.jl")
include("algorithms/primaldual.jl")
include("algorithms/davisyin.jl")
include("algorithms/lilin.jl")

end # module


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

module TestOptimizer

using Base.Test

using MXNet
using MXNet.mx.LearningRate
using MXNet.mx.Momentum

function test_fixed_η()
  info("Optimizer::LearningRate::Fixed")
  x = LearningRate.Fixed(.42)
  @test get(x) == .42
  update!(x)
  @test get(x) == .42
end  # function test_fixed_η

function check_η_decay(x)
  info("Optimizer::LearningRate::$x")

  η = get(x)
  @test η == 1

  for i ∈ 1:5
    update!(x)
    η′ = get(x)
    @test η′ < η
    η = η′
  end
end  # function check_η_decay

test_exp_η() = LearningRate.Exp(1) |> check_η_decay

test_inv_η() = LearningRate.Inv(1) |> check_η_decay

function test_μ_null()
  info("Optimizer::Momentum::Null")
  x = Momentum.Null()
  @test iszero(get(x))
end

function test_μ_fixed()
  info("Optimizer::Momentum::Fixed")
  x = Momentum.Fixed(42)
  @test get(x) == 42
end

@testset "Optimizer Test" begin
    @testset "LearningRate Test" begin
      test_fixed_η()
      test_exp_η()
      test_inv_η()
    end

    @testset "Momentum Test" begin
      test_μ_null()
      test_μ_fixed()
    end
end

end  # module TestOptimizer

max_stars_repo_name: Vikas-kum/incubator-mxnet

max_issues_repo_name: Vikas-kum/incubator-mxnet

max_forks_repo_name: Vikas-kum/incubator-mxnet


text:
using LinearAlgebra

@testset "Zhang-Yeung inequality" begin
    n = 4
    G = polymatroidcone(n, Polyhedra.DefaultLibrary{Float64}(lp_solver))
    # This vector is not entropic but it is a polymatroid
    invalidf = invalidfentropy(12)

    @testset "Basic tests" begin
        @test matusrentropy(1, 14) in G
        @test matusrentropy(1, 23) in G
        @test matusrentropy(1, 24) in G
        @test matusrentropy(2, 3) in G
        @test matusrentropy(2, 4) in G
        @testset "Non-entropic but polymatroid vector invalidf" begin
            @test invalidf in G
        end
    end

    # Let's cut it out !

    I = set(1)
    J = set(2)
    K = set(3)
    L = set(4)
    zhangyeungineq = 3(nonnegative(n,union(I,K)) + nonnegative(n,union(I,L)) + nonnegative(n,union(K,L))) +
                    nonnegative(n,union(J,K)) + nonnegative(n,union(J,L)) -
                    (nonnegative(n,I) + 2(nonnegative(n,K) + nonnegative(n,L)) + nonnegative(n,union(I,J)) +
                    4nonnegative(n,union(I,union(K,L))) + nonnegative(n,union(J,union(K,L))))
    @testset "Alternative expression" begin
        @test zhangyeunginequality() == zhangyeungineq
    end

    @testset "The Zhang-Yeung inequality is a new inequality..." begin
        @test !(zhangyeungineq in G)
    end
    @testset "...and it sees that invalidf is not a valid entropy" begin
        @test dot(zhangyeungineq, invalidf) == -1
    end
    # Actually, invalidf is even the certificate returned by CDD
    # to show that the Zhang-Yeung inequality is new
    #certificate = redundant(zhangyeungineq, G)[2]
    #certificate *= 2 / certificate[1]
    #@test PrimalEntropy{15,Int}(Array{Int}(certificate)) == invalidf
end


text:
include("POD_experiment/FLay02H.jl")

@testset "User limit testing" begin

@testset "Knapsack 50% limit" begin
    println("==================================")
    println("KNAPSACK 50%")
    println("==================================")

    m = Model(with_optimizer(
        Juniper.Optimizer,
        DefaultTestSolver(;traverse_strategy=:DBFS,branch_strategy=:MostInfeasible,mip_gap=0.5))
    )

    v = [10,20,12,23,42]
    w = [12,45,12,22,21]
    @variable(m, x[1:5], Bin)

    @objective(m, Max, dot(v,x))

    @NLconstraint(m, sum(w[i]*x[i]^2 for i=1:5) <= 45)

    status = solve(m)
    objval = JuMP.objective_value(m)
    println("Obj: ", objval)
    best_bound_val = JuMP.objective_bound(m)
    gap_val = getobjgap(m)

    @test status == MOI.OBJECTIVE_LIMIT

    @test best_bound_val >= objval
    @test 0.1 <= gap_val <= 0.5
end

@testset "FLay02H time limit" begin
    println("==================================")
    println("FLay02H time limit")
    println("==================================")

    m = get_FLay02H()

    set_optimizer(m, with_optimizer(
        Juniper.Optimizer,
        DefaultTestSolver(
            branch_strategy=:StrongPseudoCost,
            time_limit = 5,
            incumbent_constr = true)
    ))

    status = solve(m)

    @test status == MOI.LOCALLY_SOLVED || status == MOI.TIME_LIMIT
    @test getsolvetime(m) <= 15 # it might be a bit higher than 5s
end

end


text:
secondtoms(x::Float64) = Millisecond(round(Int, 1000x))
secondtoms(x::AbstractString) = secondtoms(parse(Float64, x))

function parsetime(x)
    xs = split(x, ':')
    n = length(xs)
    if n == 1
        secondtoms(x)
    elseif n == 2
        Millisecond(Minute(xs[1])) + secondtoms(xs[2])
    else
        Millisecond(Hour(xs[1])) + Millisecond(Minute(xs[2])) + secondtoms(xs[3])
    end
end

function goodtime(x)
    try
        return parsetime(x)
    catch
        return nothing
    end
end

const times_dialog = Dict("start" => Dialog(), "stop" => Dialog())

function get_time(vf, type)
    @label beginning
    __time = requestᵐ("""When in "$(vf.name)" did this POI $(type)?""", times_dialog[type])
    _time = goodtime(__time)
    if _time ≡ nothing
        println("Malformed time. Try again…")
        @goto beginning
    end
    if _time > duration(vf)
        println("Specified time is longer than the duration of this video file (", Time(0) + duration(vf), "). Try again…")
        @goto beginning
    end
    _time
end

isgoodvideo(x, y) = x.video == y.video && x.index ≥ y.index

function getstop(vfs, _start)
    i = if length(vfs) == 1
        1
    else
        options = getfield.(vfs, :name)
        menu = RadioMenu(options)
        if length(options) > 1
            request("In which video file did this POI stop?", menu)
        else
            1
        end
    end
    @label stop_time
    _stop = get_time(vfs[i], "stop")
    if i == 1 && _stop < _start
        println("Stoping time cannot come before starting time. Try again…")
        @goto stop_time
    end
    _stop, i
end

const poi_comment_dialog = Dialog()

function findvideoindex(vf)
    for (j, v) in enumerate(videos)
        i = findfirst(isequal(vf), files(v))
        if !isnothing(i)
            return v, i, j
        end
    end
    error("couldn't find video file: $vf")
end

function newinterval(ask_stop::Bool)
    @label newintervall
    i = if length(videofile_menu.options) > 1
        requestᵐ("In which video file did this POI start?", videofile_menu)
    else
        1
    end
    start_video = videofiles[i]
    video, i, _ = findvideoindex(start_video)
    _start = get_time(start_video, "start")
    start = mapreduce(duration, +, files(video)[1:i-1], init = _start)
    stop = if ask_stop
        _stop, j = getstop(files(video)[i:end], _start)
        i += j - 1
        mapreduce(duration, +, files(video)[1:i-1], init = _stop)
    else
        nothing
    end
    comment = requestᵐ("Comments about this specific POI?", poi_comment_dialog)
    try
        Temporal(video, AbstractPeriod(start, stop), comment)
    catch ex
        @warn "something was wrong" ex
        @goto newintervall
    end
end


text:
using ValueFunctionSolver
using Base.Test

# write your own tests here
# test indexing

I1 = 2
I2 = 3
I3 = 4
X = [i1*i2*i3 for i1=1:I1, i2=1:I2, i3=1:I3]
Xlong = X[:]
for i1=1:I1, i2=1:I2, i3=1:I3
    @test X[i1, i2, i3] == Xlong[index_vcattensor([i1, i2, i3], [I1, I2, I3])]
end

xs = [linspace(0., d*.2, d+1) for d = 1:5]
ws = [ones(length(xs[d]))/length(xs[d]) for d = 1:5]
xgrid, w = tensor_grid(xs,ws)


text:
Logger.log("""
              2018-04-19:
                            config.app_is_api is deprecated and will be removed in a future update.
                            Please delete all references to it, including in config/app.jl.
                            Failing to remove all references to config.app_is_api will cause the app to crash.
                            """, :warn)


text:
module SpinModule
	export ## spin_info.jl
	SpinInfo,
	isZeroTemperature,
	isClassical,
	getModelname,
	getStates,
	isSymmetricFactorization,
	getTemperature,
	getExternalfield,
	getEnvParameters,
	setTemperature!,
	setExternalfield!,
	setEnvParameters!,
	testSpinInfo,
		## spin_trotter.jl
	TrotterInfo,
	getTrotterparameter,
	getTrotteriteration,
	getTrotterlayers,
	getBeta,
	setTrotterparameter!,
	setTrotterIteration!,
	iteration2layer,
	testTrotterInfo,
		## spin_model.jl
	SpinModel,
	buildSpinSystem,
	getMeasureOperator,
	getFactorW,
	isPottsModel,
	getHamiltonian,
	testSpinModel,
		## spin_classical.jl
	ClassicalSpinModel,
	getBeta,
		## spin_classical_ising.jl
	IsingModel,
	initialize!,
	testIsingModel,
		## spin_classical_potts.jl
	PottsModel,
	testPottsModel,
		## spin_classical_clock.jl
	ClockModel,
	testClockModel,
		## spin_quantum.jl
	QuantumSpinModel,
	checkInfo!,
		## spin_quantum_ising.jl
	QuantumIsingModel,
	getFactorWp,
	testQuantumIsingModel

	include("spin_info.jl")
	include("spin_trotter.jl")
	include("spin_model.jl")
	include("spin_classical.jl")
	include("spin_classical_ising.jl")
	include("spin_classical_potts.jl")
	include("spin_classical_clock.jl")
	include("spin_quantum.jl")
	include("spin_quantum_ising.jl")
end


text:
using DEBmicroTrait
using CSV, DataFrames, Statistics
using JLD
using HypothesisTests

########################################
# I/O
df_isolates     = CSV.read("/Users/glmarschmann/.julia/dev/DEBmicroTrait/final_manuscript/files/isolates2traits.csv", DataFrame, missingstring="N/A")
########################################

########################################
Genome_size     = convert(Array{Float64,1}, df_isolates.Genome_size)
V_cell          = DEBmicroTrait.genome_size_to_cell_volume(Genome_size)
Min_gen_time    = df_isolates.Min_gen_time
gmax            = log(2)./Min_gen_time
Gram_stain      = convert(Array{String,1}, df_isolates.gram_stain)
########################################

########################################
γ_V_0 = DEBmicroTrait.max_specific_death_rate(Min_gen_time::Vector{Float64})
########################################

########################################
dry_mass        = 0.47*DEBmicroTrait.cell_volume_to_dry_mass(V_cell, gmax, Gram_stain)
ρ_bulk          = 1.21 # g/cm^3
Bio_0           = 1e9*1e6*ρ_bulk*dry_mass./12.011
γ_V_1           = median(Bio_0)
########################################

########################################
# I/O
save("/Users/glmarschmann/.julia/dev/DEBmicroTrait/final_manuscript/files/isolates_turnover.jld", "gV0", γ_V_0, "gV1", γ_V_1)
########################################

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait


text:
using Interpolations, Test

@testset "==" begin

    # issue #333
    knots = ([1,1.2],)
    vals = [1.0, 2.0]
    scheme = Gridded(Linear())
    itp = interpolate(knots, vals , scheme)
    @test itp == itp
    @test itp == deepcopy(itp)
    @test itp != interpolate(([1,1.3],), vals, scheme)
    @test itp != interpolate(knots, [1.0, 3.0], scheme)
    @test itp != interpolate(knots, vals, Gridded(Constant()))

    scheme = BSpline(Quadratic(Reflect(OnCell())))
    vals = [1.0, 2.0]
    itp = interpolate(vals, scheme)
    @test itp == itp
    @test itp == deepcopy(itp)
    @test itp != interpolate(vals, BSpline(Linear()))
    @test itp != interpolate(randn(2), scheme)

    vals = rand(Float64, 2,3)
    itp1 = interpolate(vals, BSpline(Quadratic(Flat(OnGrid()))))
    itp2 = interpolate(Float64.(vals), BSpline(Quadratic(Flat(OnGrid()))))
    @test itp1 == itp2
end

@testset "Core" begin
    A = reshape([0], 1, 1, 1, 1, 1)
    wis = ntuple(d->Interpolations.WeightedAdjIndex(1, (1,)), ndims(A))
    @test @inferred(A[wis...]) === 0
    wis = ntuple(d->Interpolations.WeightedAdjIndex(1, (1.0,)), ndims(A))
    @test @inferred(A[wis...]) === 0.0
    wis = ntuple(d->Interpolations.WeightedArbIndex((1,), (1,)), ndims(A))
    @test @inferred(A[wis...]) === 0
    wis = ntuple(d->Interpolations.WeightedArbIndex((1,), (1.0,)), ndims(A))
    @test @inferred(A[wis...]) === 0.0
    wis = ntuple(d->Interpolations.WeightedArbIndex((1,1), (1,0)), ndims(A))
    @test @inferred(A[wis...]) === 0
    wis = ntuple(d->Interpolations.WeightedArbIndex((1,1), (1.0,0.0)), ndims(A))
    @test @inferred(A[wis...]) === 0.0

    wi = Interpolations.WeightedAdjIndex(2, (0.2, 0.8))
    @test wi*2 === Interpolations.WeightedAdjIndex(2, (0.4, 1.6))
    @test wi/2 === Interpolations.WeightedAdjIndex(2, (0.1, 0.4))
    @test Interpolations.indextuple(wi) == (2, 3)
    wi = Interpolations.WeightedArbIndex((2, -1), (0.2, 0.8))
    @test wi*2 === Interpolations.WeightedArbIndex((2, -1), (0.4, 1.6))
    @test wi/2 === Interpolations.WeightedArbIndex((2, -1), (0.1, 0.4))
    @test Interpolations.indextuple(wi) == (2, -1)
end


text:
using Probably
using Probably: apply!

# Here's a function of a single bit.
# We're interested in whether the function is constant. On a classical computer
# we have to evaluate `f` twice to see `f(true)` and `f(false)`. On a quantum
# computer can evaluate `f` just once.

# Try uncommenting different definitions of `f` to see the result change.

# Here are some variable functions:
f(x) = !x
# f(x) = x

# Here's a constant function (it ignores its input):
# f(x) = true

# Generate the reversable unitary matrix for `f`:
U = classical((x, y) -> (x, f(x) ⊻ y), 2)

# Verify that the function behaves as expected:
x = true # or false
a = QBool(x)
b = QBool()
apply!(U, a, b)
measure(b) == f(x)

# Now we apply `f` with the inputs in superposition.
a = hadamard(QBool(false)) # a = |+⟩
b = hadamard(QBool(true))  # b = |-⟩
apply!(U, a, b)
hadamard(a)

# The output will always be `false` if `f` is constant.
measure(a)


text:
"""
    get_node(node::Node, name::String)
    get_node(node::Node, id::Int)

Get a node in an mtg by name or id. If names are not unique in the MTG, the function
will return the first it finds.

# Examples

```julia
file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))),"test","files","simple_plant.mtg")
mtg = read_mtg(file)
node_6 = get_node(mtg, "node_6")
node_6_2 = get_node(mtg, 6)
```
"""
function get_node(node::Node, name::String)
    if node.name == name
        return node
    else
        if !isleaf(node)
            for chnode in ordered_children(node)
                rnode = get_node(chnode, name)
                if rnode !== nothing
                    return rnode
                end
            end
        end
    end
end

function get_node(node::Node, id::Int)
    if node.id == id
        return node
    else
        if !isleaf(node)
            for chnode in ordered_children(node)
                rnode = get_node(chnode, id)
                if rnode !== nothing
                    return rnode
                end
            end
        end
    end
end


text:
# These can be changed, to trade off better performance for space
const global maxallowedprobe = 16
const global maxprobeshift   = 6

mutable struct HashIndices{T} <: AbstractIndices{T}
    slots::Array{UInt8,1}
    inds::Array{T,1}
    ndel::Int
    count::Int
    idxfloor::Int  # an index <= the indices of all used slots
    maxprobe::Int
end

HashIndices() = HashIndices{Any}()

"""
    HashIndices{I}()

Construct an empty `HashIndices` with indices of type `I`. This container uses hashes for
fast lookup, and is insertable. (See `isinsertable`).
"""
function HashIndices{T}(; sizehint::Int = 16) where {T}
    sz = Base._tablesz(sizehint)
    HashIndices{T}(zeros(UInt8, sz), Vector{T}(undef, sz), 0, 0, 1, 0)
end

## Constructors

"""
    HashIndices(iter)
    HashIndices{I}(iter)

Construct a `HashIndices` with indices from iterable container `iter`.
"""
function HashIndices(iter)
    if Base.IteratorEltype(iter) === Base.EltypeUnknown()
        # TODO: implement automatic widening from iterators of Base.EltypeUnkown
        iter = collect(iter)
    end

    return HashIndices{eltype(iter)}(iter)
end

function HashIndices{T}(iter) where {T}
    iter_size = Base.IteratorSize(iter)
    if iter_size isa Union{Base.HasLength, Base.HasShape}
        h = HashIndices{T}(; sizehint = length(iter)*2)
    else
        h = HashIndices{T}()
    end

    for i in iter
        insert!(h, i) # should this be `set!` or `insert!`?
    end

    return h
end

function Base.copy(h::HashIndices{T}) where {T}
    return HashIndices{T}(copy(h.slots), copy(h.inds), h.ndel, h.count, h.idxfloor, h.maxprobe)
end

Base.empty(::HashIndices, ::Type{T}) where {T} = HashIndices{T}()

## Length
Base.length(h::HashIndices) = h.count

## Token interface

istokenizable(::HashIndices) = true
tokentype(::HashIndices) = Int

@propagate_inbounds isslotempty(h::HashIndices, i::Int) = h.slots[i] == 0x0
@propagate_inbounds isslotfilled(h::HashIndices, i::Int) = h.slots[i] == 0x1
@propagate_inbounds isslotdeleted(h::HashIndices, i::Int) = h.slots[i] == 0x2 # deletion marker/tombstone

istokenassigned(h::HashIndices, i::Int) = isslotfilled(h, i)

# iteratetoken

function skip_deleted(h::HashIndices, i)
    L = length(h.slots)
    @inbounds while i <= L && !isslotfilled(h, i)
        i += 1
    end
    return i
end

@propagate_inbounds function iteratetoken(h::HashIndices{T}) where {T}
    idx = skip_deleted(h, h.idxfloor)
    h.idxfloor = idx # An optimization to skip unnecessary elements when iterating multiple times

    if idx > length(h.inds)
        return nothing
    else
        return (idx, idx + 1)
    end
end

@propagate_inbounds function iteratetoken(h::HashIndices{T}, idx::Int) where {T}
    idx = skip_deleted(h, idx)

    if idx > length(h.inds)
        return nothing
    else
        return (idx, idx + 1)
    end
end

# gettoken
function hashtoken(key, sz::Int)
    # Given key what is the hash slot? sz is a power of two
    (((hash(key)%Int) & (sz-1)) + 1)::Int
end

function gettoken(h::HashIndices{T}, key::T) where {T}
    sz = length(h.inds)
    iter = 0
    maxprobe = h.maxprobe
    token = hashtoken(key, sz)
    keys = h.inds

    @inbounds while true
        if isslotempty(h, token)
            break
        end
        if !isslotdeleted(h, token) && (key === keys[token] || isequal(key, keys[token]))
            return (true, token)
        end

        token = (token & (sz-1)) + 1
        iter += 1
        iter > maxprobe && break
    end
    return (false, 0)
end

# gettokenvalue
@propagate_inbounds function gettokenvalue(h::HashIndices, token::Int)
    return h.inds[token]
end

# insertable interface
isinsertable(::HashIndices) = true

function Base.empty!(h::HashIndices{T}) where {T}
    fill!(h.slots, 0x0) # It should be OK to reduce this back to some smaller size.
    sz = length(h.slots)
    empty!(h.inds)
    resize!(h.inds, sz)
    h.ndel = 0
    h.count = 0
    h.idxfloor = 1
    return h
end

function Base.rehash!(h::HashIndices, newsz::Int = length(h.inds))
    _rehash!(h, nothing, newsz)
    return h
end

function _rehash!(h::HashIndices{T}, oldv::Union{Nothing, Vector}, newsz::Int) where {T}
    olds = h.slots
    oldk = h.inds
    sz = length(olds)
    newsz = Base._tablesz(newsz)
    h.idxfloor = 1
    if h.count == 0
        resize!(h.slots, newsz)
        fill!(h.slots, 0)
        resize!(h.inds, newsz)
        error()
        oldv === nothing || resize!(oldv, newsz)
        h.ndel = 0
        return oldv
    end

    slots = zeros(UInt8, newsz)
    keys = Vector{T}(undef, newsz)
    vals = oldv === nothing ? nothing : Vector{eltype(oldv)}(undef, newsz)
    count = 0
    maxprobe = h.maxprobe

    for i ∈ 1:sz
        @inbounds if olds[i] == 0x1
            k = oldk[i]
            v = vals === nothing ? nothing : oldv[i]
            index0 = index = hashtoken(k, newsz)
            while slots[index] != 0
                index = (index & (newsz-1)) + 1
            end
            probe = (index - index0) & (newsz-1)
            probe > maxprobe && (maxprobe = probe)
            slots[index] = 0x1
            keys[index] = k
            vals === nothing || (vals[index] = v)
            count += 1
        end
    end

    h.slots = slots
    h.inds = keys
    h.count = count
    h.ndel = 0
    h.maxprobe = maxprobe

    return vals
end

Base.sizehint!(h::HashIndices, newsz::Int) = _sizehint!(h, nothing, newsz)

function _sizehint!(h::HashIndices{T}, values::Union{Nothing, Vector}, newsz::Int) where {T}
    oldsz = length(h.slots)
    if newsz <= oldsz
        # TODO: shrink
        # be careful: rehash!() assumes everything fits. it was only designed
        # for growing.
        return hash
    end
    # grow at least 25%
    newsz = min(max(newsz, (oldsz*5)>>2),
                Base.max_values(T))
    return _rehash!(h, values, newsz)
end

function gettoken!(h::HashIndices{T}, key::T) where {T}
    (token, _) = _gettoken!(h, nothing, key) # This will make sure a slot is available at `token` (or `-token` if it is new)

    if token < 0
        @inbounds (token, _) = _insert!(h, nothing, key, -token) # This will fill the slot with `key`
        return (false, token)
    else
        return (true, token)
    end
end

# get the index where a key is stored, or -pos if not present
# and the key would be inserted at pos
# This version is for use by insert!, set! and get!
function _gettoken!(h::HashIndices{T}, values::Union{Nothing, Vector}, key::T) where {T}
    sz = length(h.inds)
    iter = 0
    maxprobe = h.maxprobe
    token = hashtoken(key, sz)
    avail = 0
    keys = h.inds

    # Search of the key is present or if there is a deleted slot `key` could fill.
    @inbounds while true
        if isslotempty(h, token)
            if avail < 0
                return (avail, values)
            end
            return (-token, values)
        end

        if isslotdeleted(h, token)
            if avail == 0
                # found an available deleted slot, but we need to keep scanning
                # in case `key` already exists in a later collided slot.
                avail = -token
            end
        elseif key === keys[token] || isequal(key, keys[token])
            return (token, values)
        end

        token = (token & (sz-1)) + 1
        iter += 1
        iter > maxprobe && break
    end

    avail < 0 && return (avail, values)

    # The key definitely isn't present, but a slot may become available if we increase
    # `maxprobe` (up to some reasonable global limits).
    maxallowed = max(maxallowedprobe, sz>>maxprobeshift)

    @inbounds while iter < maxallowed
        if !isslotfilled(h,token)
            h.maxprobe = iter
            return (-token, values)
        end
        token = (token & (sz-1)) + 1
        iter += 1
    end

    # If we get here, then all the probable slots are filled, and the only recourse is to
    # increase the size of the hash map and try again
    values = _rehash!(h, values, h.count > 64000 ? sz*2 : sz*4)
    return _gettoken!(h, values, key)
end

@propagate_inbounds function _insert!(h::HashIndices{T}, values::Union{Nothing, Vector}, key::T, token::Int) where {T}
    h.slots[token] = 0x1
    h.inds[token] = key
    h.count += 1
    if token < h.idxfloor
        h.idxfloor = token
    end

    # TODO revisit this...
    #=
    sz = length(h.inds)
    # Rehash now if necessary
    if h.ndel >= ((3*sz)>>2) || h.count*3 > sz*2
        # > 3/4 deleted or > 2/3 full
        values = _rehash!(h, values, h.count > 64000 ? h.count*2 : h.count*4)
        (_, token) = gettoken(h, key)
    end
    =#

    return (token, values)
end

function deletetoken!(h::HashIndices{T}, token::Int) where {T}
    h.slots[token] = 0x2
    isbitstype(T) || ccall(:jl_arrayunset, Cvoid, (Any, UInt), h.inds, token-1)

    h.ndel += 1
    h.count -= 1
    return h
end

# Since deleting elements doesn't mess with iteration, we can use `unsafe_filter!``
Base.filter!(pred, h::HashIndices) = Base.unsafe_filter!(pred, h)

# The default insertable indices
Base.empty(d::AbstractDictionary, ::Type{T}) where {T} = HashIndices{T}()


text:
using SimpleTraits
using Phylo.API
using AxisArrays

@traitfn function _treehistory(tree::T, node::N1) where
    {N1, NL, N, B, T <: AbstractTree{OneTree, <:Rooted, NL, N, B};
     PreferBranchObjects{T}}
    branches = B[]
    nodestoprocess = N1[node]
    nodesprocessed = N1[]
    while !isempty(nodestoprocess)
        nextnode = pop!(nodestoprocess)
        push!(nodesprocessed, nextnode)
        if hasinbound(tree, nextnode)
            push!(branches, _getinbound(tree, nextnode))
            push!(nodestoprocess, _getparent(tree, nextnode))
        end
    end
    return branches, nodesprocessed
end

@traitfn function _treehistory(tree::T, node::N1) where
    {N1, NL, N, B, T <: AbstractTree{OneTree, <:Rooted, NL, N, B};
     !PreferBranchObjects{T}}
    branches = Int[]
    nodestoprocess = N1[node]
    nodesprocessed = N1[]
    while !isempty(nodestoprocess)
        nextnode = pop!(nodestoprocess)
        push!(nodesprocessed, nextnode)
        if hasinbound(tree, nextnode)
            push!(branches, _getinbound(tree, nextnode))
            push!(nodestoprocess, _getparent(tree, nextnode))
        end
    end
    return branches, nodesprocessed
end

@traitfn function _treefuture(tree::T, node::N1) where
    {N1, NL, N, B, T <: AbstractTree{OneTree, <:Rooted, NL, N, B};
     PreferBranchObjects{T}}
    branches = B[]
    nodestoprocess = N1[node]
    nodesprocessed = N1[]
    while !isempty(nodestoprocess)
        nextnode = pop!(nodestoprocess)
        push!(nodesprocessed, nextnode)
        append!(branches, _getoutbounds(tree, nextnode))
        append!(nodestoprocess, _getchildren(tree, nextnode))
    end
    return branches, nodesprocessed
end

@traitfn function _treefuture(tree::T, node::N1) where
    {N1, NL, N, B, T <: AbstractTree{OneTree, <:Rooted, NL, N, B};
     !PreferBranchObjects{T}}
    branches = Int[]
    nodestoprocess = N1[node]
    nodesprocessed = N1[]
    while !isempty(nodestoprocess)
        nextnode = pop!(nodestoprocess)
        push!(nodesprocessed, nextnode)
        append!(branches, _getoutbounds(tree, nextnode))
        append!(nodestoprocess, _getchildren(tree, nextnode))
    end
    return branches, nodesprocessed
end

"""
    branchhistory(tree::AbstractTree, node)

Find the branch route between a node on a tree and its root
"""
function branchhistory end
@traitfn function branchhistory(tree::T, node::N) where
    {N, T <: AbstractTree{OneTree, <:Rooted}; MatchNodeType{T, N}}
    return _treehistory(tree, node)[1]
end
@traitfn function branchhistory(tree::T, node::N) where
    {N, T <: AbstractTree{OneTree, <:Rooted, N}; !MatchNodeType{T, N}}
    return _treehistory(tree, getnode(tree, node))[1]
end

"""
    nodehistory(tree::AbstractTree, node)

Find the node route between a node on a tree and its root
"""
function nodehistory end
@traitfn function nodehistory(tree::T, node::N) where
    {N, T <: AbstractTree{OneTree, <:Rooted}; MatchNodeType{T, N}}
    return _treehistory(tree, node)[2]
end
@traitfn function nodehistory(tree::T, node::N) where
    {N, T <: AbstractTree{OneTree, <:Rooted, N}; !MatchNodeType{T, N}}
    return getnodename.(tree, _treehistory(tree, getnode(tree, node))[2])
end

"""
    branchfuture(tree::AbstractTree, node)

Find the branches between a node on a tree and its leaves
"""
function branchfuture end
@traitfn function branchfuture(tree::T, node::N) where
    {N, T <: AbstractTree{OneTree, <:Rooted}; MatchNodeType{T, N}}
    return _treefuture(tree, node)[1]
end
@traitfn function branchfuture(tree::T, node::N) where
    {N, T <: AbstractTree{OneTree, <:Rooted, N}; !MatchNodeType{T, N}}
    return _treefuture(tree, getnode(tree, node))[1]
end

"""
    nodefuture(tree::AbstractTree, node)

Find the nodes between a node on a tree and its leaves
"""
function nodefuture end
@traitfn function nodefuture(tree::T, node::N) where
    {N, T <: AbstractTree{OneTree, <:Rooted}; MatchNodeType{T, N}}
    return _treefuture(tree, node)[2]
end
@traitfn function nodefuture(tree::T, node::N) where
    {N, T <: AbstractTree{OneTree, <:Rooted, N}; !MatchNodeType{T, N}}
    return getnodename.(tree, _treefuture(tree, getnode(tree, node))[2])
end

"""
    branchroute(tree::AbstractTree, node1, node2)

Find the branch route between two nodes on a tree
"""
function branchroute end
@traitfn function branchroute(tree::T, node1::N, node2::N) where
    {RT, N, T <: AbstractTree{OneTree, RT}; !MatchNodeType{T, N}}
    return branchroute(tree, getnode(tree, node1), getnode(tree, node2))
end
@traitfn function branchroute(tree::T, node1::N, node2::N) where
    {RT, N, T <: AbstractTree{OneTree, RT}; MatchNodeType{T, N}}
    branches1, nodes1 = _treehistory(tree, node1)
    branches2, nodes2 = _treehistory(tree, node2)
    nodes1[end] == nodes2[end] ||
        return error("No route between nodes")
    common = branches1 ∩ branches2
    nodes1[end] == nodes2[end] || error("No route between nodes in tree")
    return append!([b for b in branches1 if b ∉ common],
                   [b for b in reverse(branches2) if b ∉ common])
end

"""
    noderoute(tree::AbstractTree, node1, node2)

Find the node route between two nodes on a tree
"""
function noderoute end
@traitfn function noderoute(tree::T, node1::N, node2::N) where
    {RT, N, T <: AbstractTree{OneTree, RT, N}; !MatchNodeType{T, N}}
    route = noderoute(tree, getnode(tree, node1), getnode(tree, node2))
    return [getnodename(tree, n) for n in route]
end
@traitfn function noderoute(tree::T, node1::N, node2::N) where
    {RT, N, T <: AbstractTree{OneTree}; MatchNodeType{T, N}}
    branches1, nodes1 = _treehistory(tree, node1)
    branches2, nodes2 = _treehistory(tree, node2)
    nodes1[end] == nodes2[end] || error("No route between nodes in tree")
    common = nodes1[end]
    while min(length(nodes1), length(nodes2)) > 0 && nodes1[end] == nodes2[end]
        common = nodes1[end]
        pop!(nodes1)
        pop!(nodes2)
    end
    push!(nodes1, common)
    return append!(nodes1, reverse(nodes2))
end


text:
#! /usr/bin/env julia

symlink("link_me.txt", "im_a_link.txt")

println(islink("im_a_link.txt"))
println(islink("link_me.txt"))

max_stars_repo_name: sjkelly/JuliaPlayground

max_issues_repo_name: sjkelly/JuliaPlayground

max_forks_repo_name: sjkelly/JuliaPlayground


text:
module RocketErrorIfNotperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: error_if_not()" begin

    println("Testing: operator error_if_not()")

    run_proxyshowcheck("ErrorIf", error_if_not((d) -> true))

    run_testset([
        (
            source      = from(1:42) |> error_if_not((d) -> d <= 2),
            values      = @ts([ 1, 2, e ]),
            source_type = Int
        ),
        (
            source      = from(1:42) |> error_if_not((d) -> d <= 2, (d) -> "Error $d"),
            values      = @ts([ 1, 2, e("Error 3") ]),
            source_type = Int
        ),
        (
            source      = from(1:42) |> error_if_not((d) -> d <= 2) |> catch_error((err, obs) -> of(3)),
            values      = @ts([ 1, 2, 3, c ]),
            source_type = Int
        ),
        (
            source      = completed(Int) |> error_if_not((d) -> d <= 2),
            values      = @ts(c),
            source_type = Int
        ),
        (
            source      = faulted(String, "e") |> error_if_not((d) -> d <= 2),
            values      = @ts(e("e")),
            source_type = String
        ),
        (
            source = never() |> error_if_not((d) -> d <= 2),
            values = @ts()
        )
    ])

    @testset "`error_if` should unsubscribe on check fail" begin
        source = Subject(Int)

        events = []

        subscription = subscribe!(source |> error_if_not((d) -> d <= 2), lambda(
            on_next = (d) -> push!(events, d),
            on_error = (e) -> push!(events, "Error"),
            on_complete = () -> push!(events, "Completed")
        ))

        @test events == [ ]

        next!(source, 1)

        @test events == [ 1 ]

        next!(source, 2)

        @test events == [ 1, 2 ]

        next!(source, 3)

        @test events == [ 1, 2, "Error" ]

        next!(source, 4)

        @test events == [ 1, 2, "Error" ]

        unsubscribe!(subscription)

        @test events == [ 1, 2, "Error" ]
    end

    @testset "`error_if` should not check anything after unsubscription" begin
        source = Subject(Int)

        events = []

        subscription = subscribe!(source |> error_if_not((d) -> d <= 2), lambda(
            on_next = (d) -> push!(events, d),
            on_error = (e) -> push!(events, "Error"),
            on_complete = () -> push!(events, "Completed")
        ))

        @test events == [ ]

        next!(source, 1)

        @test events == [ 1 ]

        next!(source, 2)

        @test events == [ 1, 2 ]

        unsubscribe!(subscription)

        next!(source, 3)

        @test events == [ 1, 2 ]

        next!(source, 4)

        @test events == [ 1, 2 ]
    end

end

end


text:
# redeclare enum values without a prefix
#
# this is useful when enum values from an underlying C library, typically prefixed for the
# lack of namespacing in C, are to be used in Julia where we do have module namespacing.
macro enum_without_prefix(enum, prefix)
    if isa(enum, Symbol)
        mod = __module__
    elseif Meta.isexpr(enum, :(.))
        mod = getfield(__module__, enum.args[1])
        enum = enum.args[2].value
    else
        error("Do not know how to refer to $enum")
    end
    enum = getfield(mod, enum)
    prefix = String(prefix)

    ex = quote end
    for instance in instances(enum)
        name = String(Symbol(instance))
        @assert startswith(name, prefix)
        push!(ex.args, :(const $(Symbol(name[length(prefix)+1:end])) = $(mod).$(Symbol(name))))
    end

    return esc(ex)
end


text:
# Multi Target Filter tests

## Continuous filter tests

# Create filter 1
linear_sys = LinearSystem(0.5*eye(2), 0.001*eye(2))
linear_obs = LinearObserver(eye(2), 0.001*eye(2))
initial_est = UncertainContinuousState([0.0, 0.0], 0.1*eye(2))
kf1 = KalmanFilter(linear_sys, linear_obs, initial_est)

# Create filter 2
linear_sys = LinearSystem(-eye(2), 0.001*eye(2))
linear_obs = LinearObserver(eye(2), 0.001*eye(2))
initial_est = UncertainContinuousState([1.0, 2.0], 0.1*eye(2))
kf2 = KalmanFilter(linear_sys, linear_obs, initial_est)

# Create filter 3
linear_sys = LinearSystem(-eye(2), 0.001*eye(2))
linear_obs = LinearObserver(eye(2), 0.001*eye(2))
initial_est = UncertainContinuousState([10.0, 20.0], 0.1*eye(2))
kf3 = KalmanFilter(linear_sys, linear_obs, initial_est)

# Build MTF
mtf = NearestNeighborMTF([kf1])
add!(mtf, [kf2, kf3])
@test_throws ArgumentError NearestNeighborMTF([])
@test_throws ArgumentError add!(mtf, kf2)

# Simulate data and process
srand(1);
simulator = make_simulator(mtf)
for i = 1:10
    for (j, sim_data) in enumerate(simulate(simulator, i*0.1))
        process!(mtf, sim_data[2])
        @test mtf.filter_bank[j].estimate.t == i*0.1
    end
end

## Discrete filter tests

# Create filter 1
linear_sys = LinearSystem(0.5*eye(2), 0.001*eye(2))
linear_obs = LinearObserver(eye(2), 0.001*eye(2))
initial_est = UncertainDiscreteState([0.0, 0.0], 0.1*eye(2))
kf1 = KalmanFilter(linear_sys, linear_obs, initial_est)

# Create filter 2
linear_sys = LinearSystem(-eye(2), 0.001*eye(2))
linear_obs = LinearObserver(eye(2), 0.001*eye(2))
initial_est = UncertainDiscreteState([1.0, 2.0], 0.1*eye(2))
kf2 = KalmanFilter(linear_sys, linear_obs, initial_est)

# Create filter 3
linear_sys = LinearSystem(-0.9*eye(2), 0.001*eye(2))
linear_obs = LinearObserver(eye(2), 0.001*eye(2))
initial_est = UncertainDiscreteState([10.0, 20.0], 0.1*eye(2))
kf3 = KalmanFilter(linear_sys, linear_obs, initial_est)

# Build MTF
mtf = NearestNeighborMTF([kf1, kf2, kf3])
@test_throws ArgumentError NearestNeighborMTF([])
@test_throws ArgumentError add!(mtf, kf2)

# Simulate data and process
srand(1);
simulator = make_simulator(mtf)
for i = 1:10
    for (j, sim_data) in enumerate(simulate(simulator, i))
        process!(mtf, sim_data[2])
        @test mtf.filter_bank[j].estimate.t == i
    end
end


text:
using NeuralQuantum, SparseArrays
using Test
using NeuralQuantum: LdagLSparseSuperopProblem, LRhoSparseOpProblem
Nsites = 4
T = Float64

prob_types = [LdagLSparseSuperopProblem, LdagLSparseOpProblem]
@testset "LdagL problem: $prob_T" for prob_T=prob_types
    lind = quantum_ising_lind(SquareLattice([Nsites],PBC=true), g=1.6, V=2.0, γ=1.0)
    net  = cached(NDM(T, Nsites, 2, 1))

    prob = prob_T(T, lind);

    v    = state(prob, net)

    L=liouvillian(lind)
    LdagL = L.data'*L.data
    rho   = dm(net, prob, false).data
    rhov  = vec(rho)

    ic = NeuralQuantum.MCMCSREvaluationCache(net, prob); zero!(ic)

    Clocs_ex  = (LdagL*rhov)./rhov
    for i=1:spacedimension(v)
        set_index!(v, i)
        NeuralQuantum.sample_network!(ic, prob, net, v)
    end
    Clocs_net_S = ic.Evalues

    @test Clocs_ex ≈ Clocs_net_S
end

Nsites = 4
T = Float64

prob_types = [LRhoSparseSuperopProblem, LRhoSparseOpProblem]
@testset "LdagL problem: $prob_T" for prob_T=prob_types
    lind = quantum_ising_lind(SquareLattice([Nsites],PBC=true), g=1.6, V=2.0, γ=1.0)
    net  = cached(NDM(T, Nsites, 2, 1))
    net  = (NDM(T, Nsites, 2, 1))

    prob = NeuralQuantum.LdagLSparseSuperopProblem(T, lind);
    probL = prob_T(T, lind);

    v    = state(prob, net)

    L=liouvillian(lind)
    LdagL = L.data'*L.data
    rho   = dm(net, prob, false).data
    rhov  = vec(rho)

    ic = NeuralQuantum.MCMCSREvaluationCache(net, prob); zero!(ic)
    icL = NeuralQuantum.MCMCSRLEvaluationCache(net, probL); zero!(icL)

    Clocs_ex  = (LdagL*rhov)./rhov
    for i=1:spacedimension(v)
        set_index!(v, i)
        NeuralQuantum.sample_network_wholespace!(ic, prob, net, v)
        NeuralQuantum.sample_network_wholespace!(icL, probL, net, v)
    end

    SREval = EvaluatedNetwork(SR(), net)
    SREvalL = EvaluatedNetwork(SR(), net)

    evaluation_post_sampling!(SREval, ic)
    evaluation_post_sampling!(SREvalL, icL)

    # those might fail in the future... just check that everything has the same
    # name
    @test SREval.L ≈ SREvalL.L
    @test SREval.S[1] ≈ SREvalL.S[1]
    @test all([l≈r for (l,r)=zip(SREval.F, SREvalL.F)])
end

@testset "LdagL problem: LRhoSparseOpProblem" begin
    lind = quantum_ising_lind(SquareLattice([Nsites],PBC=true), g=0.0, V=0.0, γ=1.0)
    net  = cached(NDM(T, Nsites, 2, 1))

    prob  = LRhoSparseSuperopProblem(T, lind);
    probL = LRhoSparseOpProblem(T, lind);

    v    = state(prob, net)

    L=liouvillian(lind)
    LdagL = L.data'*L.data
    rho   = dm(net, prob, false).data
    rhov  = vec(rho)

    ic  = NeuralQuantum.MCMCSRLEvaluationCache(net, prob); zero!(ic)
    icL = NeuralQuantum.MCMCSRLEvaluationCache(net, probL); zero!(icL)

    Clocs_ex  = abs.((L.data*rhov)./rhov).^2
    for i=1:spacedimension(v)
        set_index!(v, i)
        NeuralQuantum.sample_network!(ic, prob, net, v)
        NeuralQuantum.sample_network!(icL, probL, net, v)
    end

    SREval = EvaluatedNetwork(SR(), net)
    SREvalL = EvaluatedNetwork(SR(), net)

    evaluation_post_sampling!(SREval, ic)
    evaluation_post_sampling!(SREvalL, icL)

    @test SREval.L ≈ SREvalL.L
    @test ic.Evalues ≈ icL.Evalues
    @test all([l≈r for (l,r)=zip(SREval.F, SREvalL.F)])
    @test all([l≈r for (l,r)=zip(SREval.S, SREvalL.S)])
end

using SparseArrays, Random

# random
prob_types = [LdagLSparseSuperopProblem, LdagLSparseOpProblem]
@testset "LdagL random problem: $prob_T" for prob_T=prob_types
    Nsites=2
    lind = quantum_ising_lind(SquareLattice([Nsites],PBC=false), g=1.6, V=2.0, γ=1.0)
    net  = cached(NDM(T, Nsites, 2, 1))

    Ham  = SparseOperator(lind.H)
    Ham.data.=sprand(ComplexF64, size(Ham.data,1),size(Ham.data,2),0.5)
    Ham.data .= Ham.data + Ham.data'
    cops = jump_operators(lind)
    c1 = first(cops)
    c1.data = sprand(ComplexF64, size(c1.data,1),size(c1.data,2),0.1)
    prob = prob_T(T, Ham, cops);

    v    = state(prob, net)

    L=liouvillian(Ham, cops)
    LdagL = L.data'*L.data
    rho   = dm(net, prob, false).data
    rhov  = vec(rho)

    ic = NeuralQuantum.MCMCSREvaluationCache(net, prob); zero!(ic)

    Clocs_ex  = (LdagL*rhov)./rhov
    for i=1:spacedimension(v)
        set_index!(v, i)
        NeuralQuantum.sample_network!(ic, prob, net, v)
    end
    Clocs_net_S = ic.Evalues

    @test Clocs_ex ≈ Clocs_net_S
end

prob_types = [LRhoSparseSuperopProblem, LRhoSparseOpProblem]
@testset "LdagL random problem: $prob_T" for prob_T=prob_types
    Nsites=2
    lind = quantum_ising_lind(SquareLattice([Nsites],PBC=false), g=1.6, V=2.0, γ=1.0)
    net  = cached(NDM(T, Nsites, 2, 1))

    Ham  = SparseOperator(lind.H)
    Ham.data.=sprand(ComplexF64, size(Ham.data,1),size(Ham.data,2),0.5)
    Ham.data .= Ham.data + Ham.data'
    cops = jump_operators(lind)
    c1 = first(cops)
    c1.data = sprand(ComplexF64, size(c1.data,1),size(c1.data,2),0.1)
    prob = prob_T(T, Ham, cops);

    v    = state(prob, net)

    L=liouvillian(Ham, cops)
    rho   = dm(net, prob, false).data
    rhov  = vec(rho)

    ic = NeuralQuantum.MCMCSRLEvaluationCache(net, prob); zero!(ic)

    Clocs_ex  = abs.((L.data*rhov)./rhov).^2
    for i=1:spacedimension(v)
        set_index!(v, i)
        NeuralQuantum.sample_network!(ic, prob, net, v)
    end
    Clocs_net_S = ic.Evalues

    @test Clocs_ex ≈ Clocs_net_S
end


text:
using TensorFunctions
using Test
using TensorOperations

@testset "TensorFunctions.jl" begin
    @test TensorFunctions.isqnode(:( :b )) == true
    @test TensorFunctions.isqnode(:( b )) == false

    @test TensorFunctions.isqnodewithdim(:( :a )) == false
    @test TensorFunctions.isqnodewithdim(:( :b|1 )) == true
    @test TensorFunctions.isqnodewithdim(:( :b|size(foo,1) )) == true

    @test TensorFunctions.isintindexproduct(:( :a*5 )) == true
    @test TensorFunctions.isintindexproduct(:( :a*foo() )) == false
    @test TensorFunctions.isintindexproduct(:( :a )) == false
    @test TensorFunctions.isintindexproduct(:( :a*5*5 )) == false

    @test TensorFunctions.isrhspairedindex(:( (:a,:b|1,:c) )) == true
    @test TensorFunctions.isrhspairedindex(:( (:a,:b|foo(),:c) )) == true
    @test TensorFunctions.isrhspairedindex(:( (:a,:c) )) == true

    @test TensorFunctions.islhspairedindex(:( (:a,:b|1,:c) )) == false
    @test TensorFunctions.islhspairedindex(:( (:a,:b|foo(),:c) )) == false
    @test TensorFunctions.islhspairedindex(:( (:a,:c) )) == true

    @test TensorFunctions.isrhsindex(:( (:a,:b)  ))  == true
    @test TensorFunctions.isrhsindex(:( (:a,:b|bar())  ))  == true
    @test TensorFunctions.isrhsindex(:( :a*3  ))  == true
    @test TensorFunctions.isrhsindex(:( 5*:a  ))  == true
    @test TensorFunctions.isrhsindex(:( :a  ))  == true

    @test TensorFunctions.islhsindex(:( (:a,:b)  ))  == true
    @test TensorFunctions.islhsindex(:( (:a,:b|bar())  ))  == false
    @test TensorFunctions.islhsindex(:( :a*3  ))  == false
    @test TensorFunctions.islhsindex(:( 5*:a  ))  == false
    @test TensorFunctions.islhsindex(:( :a  ))  == true

    @test TensorFunctions.isrhstensor(:( foo[:a,:b,:c] )) == true
    @test TensorFunctions.isrhstensor(:( foo(bar)[:a,:b,:c] )) == true
    @test TensorFunctions.isrhstensor(:( foo(bar)[(:a,:b),:c] )) == true
    @test TensorFunctions.isrhstensor(:( foo(bar)[(:a|hoge(2),:b),:c] )) == true
    @test TensorFunctions.isrhstensor(:( foo(bar|>foo)[(:a|4,:b|hoge(3)),:c*6] )) == true
    @test TensorFunctions.isrhstensor(:( [:a,:b,:c] )) == false
    @test TensorFunctions.isrhstensor(:( [:a,(:b,:c)] )) == false
    @test TensorFunctions.isrhstensor(:( [(:a|4,:b|hoge(3)),:c*6] )) == false

    @test TensorFunctions.islhstensor(:( foo[:a,:b,:c] )) == false
    @test TensorFunctions.islhstensor(:( foo(bar)[:a,:b,:c] )) == false
    @test TensorFunctions.islhstensor(:( foo(bar)[(:a,:b),:c] )) == false
    @test TensorFunctions.islhstensor(:( foo(bar)[(:a|hoge(2),:b),:c] )) == false
    @test TensorFunctions.islhstensor(:( foo(bar|>hugo)[(:a|4,:b|hoge(3)),:c*6] )) == false
    @test TensorFunctions.islhstensor(:( [:a,:b,:c] )) == true
    @test TensorFunctions.islhstensor(:( [:a,(:b,:c)] )) == true
    @test TensorFunctions.islhstensor(:( [(:a|4,:b|hoge(3)),:c*6] )) == false

    @test TensorFunctions.isrhssimpletensor(:( foo[:a,:b,:c] )) == true
    @test TensorFunctions.isrhssimpletensor(:( foo(bar)[:a,:b,:c] )) == true
    @test TensorFunctions.isrhssimpletensor(:( foo(bar)[(:a,:b),:c] )) == false
    @test TensorFunctions.isrhssimpletensor(:( foo(bar)[(:a|hoge(2),:b),:c] )) == false
    @test TensorFunctions.isrhssimpletensor(:( foo(bar|>foo)[(:a|4,:b|hoge(3)),:c*6] )) == false
    @test TensorFunctions.isrhssimpletensor(:( [:a,:b,:c] )) == false
    @test TensorFunctions.isrhssimpletensor(:( [:a,(:b,:c)] )) == false
    @test TensorFunctions.isrhssimpletensor(:( [(:a|4,:b|hoge(3)),:c*6] )) == false

    @test TensorFunctions.islhssimpletensor(:( foo[:a,:b,:c] )) == false
    @test TensorFunctions.islhssimpletensor(:( foo(bar)[:a,:b,:c] )) == false
    @test TensorFunctions.islhssimpletensor(:( foo(bar)[(:a,:b),:c] )) == false
    @test TensorFunctions.islhssimpletensor(:( foo(bar)[(:a|hoge(2),:b),:c] )) == false
    @test TensorFunctions.islhssimpletensor(:( foo(bar|>hugo)[(:a|4,:b|hoge(3)),:c*6] )) == false
    @test TensorFunctions.islhssimpletensor(:( [:a,:b,:c] )) == true
    @test TensorFunctions.islhssimpletensor(:( [:a,(:b,:c)] )) == false
    @test TensorFunctions.islhssimpletensor(:( [(:a|4,:b|hoge(3)),:c*6] )) == false

    @test TensorFunctions.istensorproduct(:( A[:a,:b] * B[:a,:c,:b] * C[:d,:r] )) == true
    @test TensorFunctions.istensorproduct(:( A[:a,:b] * B[:a,(:c,:b)] * C[:d*1000,:r] )) == true
    @test TensorFunctions.istensorproduct(:( 1 * A[:a,:b] * B[:a,(:c,:b)] * C[:d*1000,:r] )) == false
    @test TensorFunctions.istensorproduct(:( foo() * A[:a,:b] * B[:a,(:c,:b)] * C[:d*1000,:r] )) == false
    @test TensorFunctions.istensorproduct(:( (A[:a,:b] * B[:a,(:c,:b)]) * C[:d*1000,:r] )) == false
    @test TensorFunctions.istensorproduct(:( (A[:a,:b] * B[:a,(:c,:b)])[:A,:B] * C[:d*1000,:r] )) == true

    @test TensorFunctions.issimpletensorproduct(:( A[:a,:b] * B[:a,:c,:b] * C[:d,:r] )) == true
    @test TensorFunctions.issimpletensorproduct(:( A[:a,:b] * B[:a,(:c,:b)] * C[:d*1000,:r] )) == false
    @test TensorFunctions.issimpletensorproduct(:( 1 * A[:a,:b] * B[:a,(:c,:b)] * C[:d*1000,:r] )) == false
    @test TensorFunctions.issimpletensorproduct(:( foo() * A[:a,:b] * B[:a,(:c,:b)] * C[:d*1000,:r] )) == false
    @test TensorFunctions.issimpletensorproduct(:( (A[:a,:b] * B[:a,(:c,:b)]) * C[:d*1000,:r] )) == false
    @test TensorFunctions.issimpletensorproduct(:( (A[:a,:b] * B[:a,(:c,:b)])[:A,:B] * C[:d,:r] )) == true

    A = randn(3,4,5,6,7,8,3,3,2,2); B = randn(4,5,2)
    @tensor C[a,c,e,f,h] := reshape(A,3, 2,2 ,5 ,6*7, 2,4 ,3*3,2,2)[a,b,c,d,e,f,g,h,i,i] * B[g,d,b]
    C = reshape(C,:,size(C)[3:end]...)
    C2 = @tensorfunc [(:a,:c),:e,:f,:h] <= A[:a,(:b,:c),:d,:e*2,(:f,:g),:h*2,:i,:i] * B[:g,:d,:b]
    @test isapprox(C, C2)

    A = randn(100); B= randn(5,5)
    @tensor C[a,bb,c] := reshape(A,2,5,10)[a,b,c] * B[b,bb]
    C2 = @tensorfunc [:a,:bb,:c] <= A[(:a,:b,:c|10)] * B[:b,:bb]
    @test isapprox(C, C2)

end


text:
module Filters

using SpecialFunctions
using LinearAlgebra, GaussQuadrature, KernelAbstractions
using KernelAbstractions.Extras: @unroll
using StaticArrays
using ..Grids
using ..Grids: Direction, EveryDirection, HorizontalDirection, VerticalDirection

using ...MPIStateArrays
using ...VariableTemplates: @vars, varsize, Vars, varsindices

export AbstractSpectralFilter, AbstractFilter
export ExponentialFilter,
    CutoffFilter, MassPreservingCutoffFilter, TMARFilter, BoydVandevenFilter

abstract type AbstractFilter end
abstract type AbstractSpectralFilter <: AbstractFilter end

"""
    AbstractFilterTarget

An abstract type representing variables that the filter
will act on
"""
abstract type AbstractFilterTarget end

"""
    vars_state_filtered(::AbstractFilterTarget, FT)

A tuple of symbols containing variables that the filter
will act on given a float type `FT`
"""
function vars_state_filtered end

"""
    compute_filter_argument!(::AbstractFilterTarget,
                             state_filter::Vars,
                             state::Vars,
                             state_auxiliary::Vars)

Compute filter argument `state_filter` based on `state`
and `state_auxiliary`
"""
function compute_filter_argument! end
"""
    compute_filter_result!(::AbstractFilterTarget,
                           state::Vars,
                           state_filter::Vars,
                           state_auxiliary::Vars)

Compute filter result `state` based on the filtered state
`state_filter` and `state_auxiliary`
"""
function compute_filter_result! end

number_state_filtered(t::AbstractFilterTarget, FT) =
    varsize(vars_state_filtered(t, FT))

"""
    FilterIndices(I)

Filter variables based on their indices `I` where `I` can
be a range or a list of indices

## Examples
```julia
FiltersIndices(1:3)
FiltersIndices(1, 3, 5)
```
"""
struct FilterIndices{I} <: AbstractFilterTarget
    FilterIndices(I::Integer...) = new{I}()
    FilterIndices(I::AbstractRange) = new{I}()
end
vars_state_filtered(::FilterIndices{I}, FT) where {I} =
    @vars(_::SVector{length(I), FT})

function compute_filter_argument!(
    ::FilterIndices{I},
    filter_state::Vars,
    state::Vars,
    aux::Vars,
) where {I}
    @unroll for s in 1:length(I)
        @inbounds parent(filter_state)[s] = parent(state)[I[s]]
    end
end

function compute_filter_result!(
    ::FilterIndices{I},
    state::Vars,
    filter_state::Vars,
    aux::Vars,
) where {I}
    @unroll for s in 1:length(I)
        @inbounds parent(state)[I[s]] = parent(filter_state)[s]
    end
end

"""
    spectral_filter_matrix(r, Nc, σ)

Returns the filter matrix that takes function values at the interpolation
`N+1` points, `r`, converts them into Legendre polynomial basis coefficients,
multiplies
```math
σ((n-N_c)/(N-N_c))
```
against coefficients `n=Nc:N` and evaluates the resulting polynomial at the
points `r`.
"""
function spectral_filter_matrix(r, Nc, σ)
    N = length(r) - 1
    T = eltype(r)

    @assert N >= 0
    @assert 0 <= Nc <= N

    a, b = GaussQuadrature.legendre_coefs(T, N)
    V = (N == 0 ? ones(T, 1, 1) : GaussQuadrature.orthonormal_poly(r, a, b))

    Σ = ones(T, N + 1)
    Σ[(Nc:N) .+ 1] .= σ.(((Nc:N) .- Nc) ./ (N - Nc))

    V * Diagonal(Σ) / V
end

"""
    modified_filter_matrix(r, Nc, σ)

Returns the filter matrix that takes function values at the interpolation
`N+1` points, `r`, converts them into Legendre polynomial basis coefficients,
multiplies
```math
σ((n-N_c)/(N-N_c))
```
against coefficients `n=Nc:N` and evaluates the resulting polynomial at the
points `r`. Unlike spectral_filter_matrix, this allows for the identity matrix,
to be applied.
"""
function modified_filter_matrix(r, Nc, σ)
    N = length(r) - 1
    T = eltype(r)

    @assert N >= 0
    @assert 0 <= Nc

    Nc > N && return Array{T}(I, N + 1, N + 1)

    a, b = GaussQuadrature.legendre_coefs(T, N)
    V = (N == 0 ? ones(T, 1, 1) : GaussQuadrature.orthonormal_poly(r, a, b))

    Σ = ones(T, N + 1)
    Σ[(Nc:N) .+ 1] .= σ.(((Nc:N) .- Nc) ./ (N - Nc))

    V * Diagonal(Σ) / V
end

"""
    ExponentialFilter(grid, Nc=0, s=32, α=-log(eps(eltype(grid))))

Returns the spectral filter with the filter function
```math
σ(η) = \exp(-α η^s)
```
where `s` is the filter order (must be even), the filter starts with
polynomial order `Nc`, and `alpha` is a parameter controlling the smallest
value of the filter function.
"""
struct ExponentialFilter{FM} <: AbstractSpectralFilter
    "filter matrices in all directions (tuple of filter matrices)"
    filter_matrices::FM

    function ExponentialFilter(
        grid,
        Nc = 0,
        s = 32,
        α = -log(eps(eltype(grid))),
    )
        dim = dimensionality(grid)

        # Support different filtering thresholds in different
        # directions (default behavior is to apply the same threshold
        # uniformly in all directions)
        if Nc isa Integer
            Nc = ntuple(i -> Nc, dim)
        elseif Nc isa NTuple{2} && dim == 3
            Nc = (Nc[1], Nc[1], Nc[2])
        end
        @assert length(Nc) == dim

        # Tuple of polynomial degrees (N₁, N₂, N₃)
        N = polynomialorders(grid)
        # In 2D, we assume same polynomial order in the horizontal
        @assert dim == 2 || N[1] == N[2]
        @assert iseven(s)
        @assert all(0 .<= Nc .<= N)

        σ(η) = exp(-α * η^s)

        AT = arraytype(grid)
        ξ = referencepoints(grid)
        filter_matrices =
            ntuple(i -> AT(spectral_filter_matrix(ξ[i], Nc[i], σ)), dim)
        new{typeof(filter_matrices)}(filter_matrices)
    end
end

"""
    BoydVandevenFilter(grid, Nc=0, s=32)

Returns the spectral filter using the logarithmic error function of
the form:
```math
σ(η) = 1/2 erfc(2*sqrt(s)*χ(η)*(abs(η)-0.5))
```
whenever s ≤ i ≤ N, and 1 otherwise. The function `χ(η)` is defined
as
```math
χ(η) = sqrt(-log(1-4*(abs(η)-0.5)^2)/(4*(abs(η)-0.5)^2))
```
if `x != 0.5` and `1` otherwise. Here, `s` is the filter order,
the filter starts with polynomial order `Nc`, and `alpha` is a parameter
controlling the smallest value of the filter function.

### References
 - [Boyd1996](@cite)
"""
struct BoydVandevenFilter{FM} <: AbstractSpectralFilter
    "filter matrices in all directions (tuple of filter matrices)"
    filter_matrices::FM

    function BoydVandevenFilter(grid, Nc = 0, s = 32)
        dim = dimensionality(grid)

        # Support different filtering thresholds in different
        # directions (default behavior is to apply the same threshold
        # uniformly in all directions)
        if Nc isa Integer
            Nc = ntuple(i -> Nc, dim)
        elseif Nc isa NTuple{2} && dim == 3
            Nc = (Nc[1], Nc[1], Nc[2])
        end
        @assert length(Nc) == dim

        # Tuple of polynomial degrees (N₁, N₂, N₃)
        N = polynomialorders(grid)
        # In 2D, we assume same polynomial order in the horizontal
        @assert dim == 2 || N[1] == N[2]
        @assert iseven(s)
        @assert all(0 .<= Nc .<= N)

        function σ(η)
            a = 2 * abs(η) - 1
            χ = iszero(a) ? one(a) : sqrt(-log1p(-a^2) / a^2)
            return erfc(sqrt(s) * χ * a) / 2
        end

        AT = arraytype(grid)
        ξ = referencepoints(grid)
        filter_matrices =
            ntuple(i -> AT(spectral_filter_matrix(ξ[i], Nc[i], σ)), dim)
        new{typeof(filter_matrices)}(filter_matrices)
    end
end

"""
    CutoffFilter(grid, Nc=polynomialorders(grid))

Returns the spectral filter that zeros out polynomial modes greater than or
equal to `Nc`.
"""
struct CutoffFilter{FM} <: AbstractSpectralFilter
    "filter matrices in all directions (tuple of filter matrices)"
    filter_matrices::FM

    function CutoffFilter(grid, Nc = polynomialorders(grid))
        dim = dimensionality(grid)

        # Support different filtering thresholds in different
        # directions (default behavior is to apply the same threshold
        # uniformly in all directions)
        if Nc isa Integer
            Nc = ntuple(i -> Nc, dim)
        elseif Nc isa NTuple{2} && dim == 3
            Nc = (Nc[1], Nc[1], Nc[2])
        end
        @assert length(Nc) == dim

        # Tuple of polynomial degrees (N₁, N₂, N₃)
        N = polynomialorders(grid)
        # In 2D, we assume same polynomial order in the horizontal
        @assert dim == 2 || N[1] == N[2]
        @assert all(0 .<= Nc .<= N)

        σ(η) = 0

        AT = arraytype(grid)
        ξ = referencepoints(grid)
        filter_matrices =
            ntuple(i -> AT(spectral_filter_matrix(ξ[i], Nc[i], σ)), dim)
        new{typeof(filter_matrices)}(filter_matrices)
    end
end

"""
    MassPreservingCutoffFilter(grid, Nc=polynomialorders(grid))

Returns the spectral filter that zeros out polynomial modes greater than or
equal to `Nc` while preserving the cell average value. Use this filter if the
jacobian is nonconstant.
"""
struct MassPreservingCutoffFilter{FM} <: AbstractSpectralFilter
    "filter matrices in all directions (tuple of filter matrices)"
    filter_matrices::FM

    function MassPreservingCutoffFilter(grid, Nc = polynomialorders(grid))
        dim = dimensionality(grid)

        # Support different filtering thresholds in different
        # directions (default behavior is to apply the same threshold
        # uniformly in all directions)
        if Nc isa Integer
            Nc = ntuple(i -> Nc, dim)
        elseif Nc isa NTuple{2} && dim == 3
            Nc = (Nc[1], Nc[1], Nc[2])
        end
        @assert length(Nc) == dim

        # Tuple of polynomial degrees (N₁, N₂, N₃)
        N = polynomialorders(grid)
        # In 2D, we assume same polynomial order in the horizontal
        @assert dim == 2 || N[1] == N[2]
        @assert all(0 .<= Nc)

        σ(η) = 0

        AT = arraytype(grid)
        ξ = referencepoints(grid)
        filter_matrices =
            ntuple(i -> AT(modified_filter_matrix(ξ[i], Nc[i], σ)), dim)
        new{typeof(filter_matrices)}(filter_matrices)
    end
end

"""
    TMARFilter()

Returns the truncation and mass aware rescaling nonnegativity preservation
filter.  The details of this filter are described in [Light2016](@cite)

Note this needs to be used with a restrictive time step or a flux correction
to ensure that grid integral is conserved.

## Examples

This filter can be applied to the 3rd and 4th fields of an `MPIStateArray` `Q`
with the code

```julia
Filters.apply!(Q, (3, 4), grid, TMARFilter())
```

where `grid` is the associated `DiscontinuousSpectralElementGrid`.
"""
struct TMARFilter <: AbstractFilter end

"""
    Filters.apply!(Q::MPIStateArray,
        target,
        grid::DiscontinuousSpectralElementGrid,
        filter::AbstractSpectralFilter;
        kwargs...)

Applies `filter` to `Q` given a `grid` and a custom `target`.

A `target` can be any of the following:
 - a tuple or range of indices
 - a tuple of symbols or strings of variable names
 - a colon (`:`) to apply to all variables
 - a custom [`AbstractFilterTarget`]

The following keyword arguments are supported for some filters:
- `direction`: for `AbstractSpectralFilter` controls if the filter is
  applied in the horizontal and/or vertical directions. It is assumed that the
  trailing dimension on the reference element is the vertical dimension and the
  rest are horizontal.
- `state_auxiliary`: if `target` requires auxiliary state to compute its argument or results.

# Examples

Specifying the `target` via indices:
```julia
Filters.apply!(Q, :, grid, TMARFilter())
Filters.apply!(Q, (1, 3), grid, CutoffFilter(grid); direction=VerticalDirection())
```

Speciying `target` via symbols or strings:
```julia
Filters.apply!(Q, (:ρ, "energy.ρe"), grid, TMARFilter())
Filters.apply!(Q, ("moisture.ρq_tot",), grid, CutoffFilter(grid);
               direction=VerticalDirection())
```
"""
function apply!(
    Q,
    target,
    grid::DiscontinuousSpectralElementGrid,
    filter::AbstractFilter;
    kwargs...,
)
    device = typeof(Q.data) <: Array ? CPU() : CUDADevice()
    event = Event(device)
    event =
        apply_async!(Q, target, grid, filter; dependencies = event, kwargs...)
    wait(device, event)
end

"""
    Filters.apply_async!(Q, target, grid::DiscontinuousSpectralElementGrid,
        filter::AbstractFilter;
        dependencies,
        kwargs...)

An asynchronous version of [`Filters.apply!`](@ref), returning an `Event`
object. `dependencies` should be an `Event` or tuple of `Event`s which need to
finish before applying the filter.

```julia
compstream = Filters.apply_async!(Q, :, grid, CutoffFilter(grid); dependencies=compstream)
wait(compstream)
```
"""
function apply_async! end

function apply_async!(
    Q,
    target::AbstractFilterTarget,
    grid::DiscontinuousSpectralElementGrid,
    filter::AbstractSpectralFilter;
    dependencies,
    state_auxiliary = nothing,
    direction = EveryDirection(),
)
    topology = grid.topology

    # Tuple of polynomial degrees (N₁, N₂, N₃)
    N = polynomialorders(grid)
    # In 2D, we assume same polynomial order in the horizontal
    dim = dimensionality(grid)
    # Currently only support same polynomial in both horizontal directions
    @assert N[1] == N[2]

    device = typeof(Q.data) <: Array ? CPU() : CUDADevice()

    nelem = length(topology.elems)
    # Number of Gauss-Lobatto quadrature points in each direction
    Nq = N .+ 1
    Nq1 = Nq[1]
    Nq2 = Nq[2]
    Nq3 = dim == 2 ? 1 : Nq[dim]

    nrealelem = length(topology.realelems)
    event = dependencies

    if direction isa EveryDirection || direction isa HorizontalDirection
        @assert dim == 2 || Nq1 == Nq2
        filtermatrix = filter.filter_matrices[1]
        event = kernel_apply_filter!(device, (Nq1, Nq2, Nq3))(
            Val(dim),
            Val(N),
            Val(vars(Q)),
            Val(isnothing(state_auxiliary) ? nothing : vars(state_auxiliary)),
            HorizontalDirection(),
            Q.data,
            isnothing(state_auxiliary) ? nothing : state_auxiliary.data,
            target,
            filtermatrix,
            ndrange = (nrealelem * Nq1, Nq2, Nq3),
            dependencies = event,
        )
    end
    if direction isa EveryDirection || direction isa VerticalDirection
        filtermatrix = filter.filter_matrices[end]
        event = kernel_apply_filter!(device, (Nq1, Nq2, Nq3))(
            Val(dim),
            Val(N),
            Val(vars(Q)),
            Val(isnothing(state_auxiliary) ? nothing : vars(state_auxiliary)),
            VerticalDirection(),
            Q.data,
            isnothing(state_auxiliary) ? nothing : state_auxiliary.data,
            target,
            filtermatrix,
            ndrange = (nrealelem * Nq1, Nq2, Nq3),
            dependencies = event,
        )
    end
    return event
end

function apply_async!(
    Q,
    target::AbstractFilterTarget,
    grid::DiscontinuousSpectralElementGrid,
    ::TMARFilter;
    dependencies,
)
    topology = grid.topology

    device = typeof(Q.data) <: Array ? CPU() : CUDADevice()

    dim = dimensionality(grid)
    N = polynomialorders(grid)
    # Currently only support same polynomial in both horizontal directions
    @assert dim == 2 || N[1] == N[2]
    Nqs = N .+ 1
    Nq = Nqs[1]
    Nqj = dim == 2 ? 1 : Nqs[2]

    nrealelem = length(topology.realelems)
    nreduce = 2^ceil(Int, log2(Nq * Nqj))

    event = dependencies
    event = kernel_apply_TMAR_filter!(device, (Nq, Nqj), (nrealelem * Nq, Nqj))(
        Val(nreduce),
        Val(dim),
        Val(N),
        Q.data,
        target,
        grid.vgeo,
        dependencies = event,
    )
    return event
end

function apply_async!(
    Q,
    target::AbstractFilterTarget,
    grid::DiscontinuousSpectralElementGrid,
    filter::MassPreservingCutoffFilter;
    dependencies,
    state_auxiliary = nothing,
    direction = EveryDirection(),
)
    topology = grid.topology

    device = typeof(Q.data) <: Array ? CPU() : CUDADevice()

    dim = dimensionality(grid)
    N = polynomialorders(grid)
    # Currently only support same polynomial in both horizontal directions
    @assert dim == 2 || N[1] == N[2]
    Nq = N .+ 1
    Nq1 = Nq[1]
    Nq2 = Nq[2]
    Nq3 = dim == 2 ? 1 : Nq[dim]

    nrealelem = length(topology.realelems)
    # parallel sum info
    nreduce = 2^ceil(Int, log2(Nq1 * Nq2 * Nq3))
    event = dependencies

    if direction isa EveryDirection || direction isa HorizontalDirection
        @assert dim == 2 || Nq1 == Nq2
        filtermatrix = filter.filter_matrices[1]
        event = kernel_apply_mp_filter!(device, (Nq1, Nq2, Nq3))(
            Val(nreduce),
            Val(dim),
            Val(N),
            Val(vars(Q)),
            Val(isnothing(state_auxiliary) ? nothing : vars(state_auxiliary)),
            HorizontalDirection(),
            Q.data,
            isnothing(state_auxiliary) ? nothing : state_auxiliary.data,
            target,
            filtermatrix,
            grid.vgeo,
            ndrange = (nrealelem * Nq1, Nq2, Nq3),
            dependencies = event,
        )
    end
    if direction isa EveryDirection || direction isa VerticalDirection
        filtermatrix = filter.filter_matrices[end]
        event = kernel_apply_mp_filter!(device, (Nq1, Nq2, Nq3))(
            Val(nreduce),
            Val(dim),
            Val(N),
            Val(vars(Q)),
            Val(isnothing(state_auxiliary) ? nothing : vars(state_auxiliary)),
            VerticalDirection(),
            Q.data,
            isnothing(state_auxiliary) ? nothing : state_auxiliary.data,
            target,
            filtermatrix,
            grid.vgeo,
            ndrange = (nrealelem * Nq1, Nq2, Nq3),
            dependencies = event,
        )
    end
    return event
end

function apply_async!(
    Q,
    indices::Union{Colon, AbstractRange, Tuple{Vararg{Integer}}},
    grid::DiscontinuousSpectralElementGrid,
    filter::AbstractFilter;
    kwargs...,
)
    if indices isa Colon
        indices = 1:size(Q, 2)
    end
    apply_async!(Q, FilterIndices(indices...), grid, filter; kwargs...)
end

function apply_async!(
    Q,
    vs::Tuple,
    grid::DiscontinuousSpectralElementGrid,
    filter::AbstractFilter;
    kwargs...,
)
    apply_async!(
        Q,
        FilterIndices(varsindices(vars(Q), vs)...),
        grid,
        filter;
        kwargs...,
    )
end

const _M = Grids._M

@doc """
    kernel_apply_filter!(::Val{dim}, ::Val{N}, direction,
                         Q, state_auxiliary, target, filtermatrix
                        ) where {dim, N}

Computational kernel: Applies the `filtermatrix` to `Q` given a
custom target `target` while preserving the cell average.

The `direction` argument is used to control if the filter is applied in the
horizontal and/or vertical reference directions.
""" kernel_apply_filter!
@kernel function kernel_apply_filter!(
    ::Val{dim},
    ::Val{N},
    ::Val{vars_Q},
    ::Val{vars_state_auxiliary},
    direction,
    Q,
    state_auxiliary,
    target::AbstractFilterTarget,
    filtermatrix,
) where {dim, N, vars_Q, vars_state_auxiliary}
    @uniform begin
        FT = eltype(Q)

        Nqs = N .+ 1
        Nq1 = Nqs[1]
        Nq2 = Nqs[2]
        Nq3 = dim == 2 ? 1 : Nqs[dim]

        if direction isa EveryDirection
            filterinξ1 = filterinξ2 = true
            filterinξ3 = dim == 2 ? false : true
        elseif direction isa HorizontalDirection
            filterinξ1 = true
            filterinξ2 = dim == 2 ? false : true
            filterinξ3 = false
        elseif direction isa VerticalDirection
            filterinξ1 = false
            filterinξ2 = dim == 2 ? true : false
            filterinξ3 = dim == 2 ? false : true
        end

        nstates = varsize(vars_Q)
        nfilterstates = number_state_filtered(target, FT)
        nfilteraux =
            isnothing(state_auxiliary) ? 0 : varsize(vars_state_auxiliary)

        # ugly workaround around problems with @private
        # hopefully will be soon fixed in KA
        l_Q2 = MVector{nstates, FT}(undef)
        l_Qfiltered2 = MVector{nfilterstates, FT}(undef)
    end

    s_Q = @localmem FT (Nq1, Nq2, Nq3, nfilterstates)
    l_Q = @private FT (nstates,)
    l_Qfiltered = @private FT (nfilterstates,)
    l_aux = @private FT (nfilteraux,)

    e = @index(Group, Linear)
    i, j, k = @index(Local, NTuple)

    @inbounds begin
        ijk = i + Nq1 * ((j - 1) + Nq2 * (k - 1))

        @unroll for s in 1:nstates
            l_Q[s] = Q[ijk, s, e]
        end

        @unroll for s in 1:nfilteraux
            l_aux[s] = state_auxiliary[ijk, s, e]
        end

        fill!(l_Qfiltered2, -zero(FT))

        compute_filter_argument!(
            target,
            Vars{vars_state_filtered(target, FT)}(l_Qfiltered2),
            Vars{vars_Q}(l_Q[:]),
            Vars{vars_state_auxiliary}(l_aux[:]),
        )

        @unroll for fs in 1:nfilterstates
            l_Qfiltered[fs] = zero(FT)
        end

        @unroll for fs in 1:nfilterstates
            s_Q[i, j, k, fs] = l_Qfiltered2[fs]
        end

        if filterinξ1
            @synchronize
            @unroll for n in 1:Nq1
                @unroll for fs in 1:nfilterstates
                    l_Qfiltered[fs] += filtermatrix[i, n] * s_Q[n, j, k, fs]
                end
            end

            if filterinξ2 || filterinξ3
                @synchronize
                @unroll for fs in 1:nfilterstates
                    s_Q[i, j, k, fs] = l_Qfiltered[fs]
                    l_Qfiltered[fs] = zero(FT)
                end
            end
        end

        if filterinξ2
            @synchronize
            @unroll for n in 1:Nq2
                @unroll for fs in 1:nfilterstates
                    l_Qfiltered[fs] += filtermatrix[j, n] * s_Q[i, n, k, fs]
                end
            end

            if filterinξ3
                @synchronize
                @unroll for fs in 1:nfilterstates
                    s_Q[i, j, k, fs] = l_Qfiltered[fs]
                    l_Qfiltered[fs] = zero(FT)
                end
            end
        end

        if filterinξ3
            @synchronize
            @unroll for n in 1:Nq3
                @unroll for fs in 1:nfilterstates
                    l_Qfiltered[fs] += filtermatrix[k, n] * s_Q[i, j, n, fs]
                end
            end
        end

        @unroll for s in 1:nstates
            l_Q2[s] = l_Q[s]
        end

        compute_filter_result!(
            target,
            Vars{vars_Q}(l_Q2),
            Vars{vars_state_filtered(target, FT)}(l_Qfiltered[:]),
            Vars{vars_state_auxiliary}(l_aux[:]),
        )

        # Store result
        ijk = i + Nq1 * ((j - 1) + Nq2 * (k - 1))
        @unroll for s in 1:nstates
            Q[ijk, s, e] = l_Q2[s]
        end

        @synchronize
    end
end

@kernel function kernel_apply_TMAR_filter!(
    ::Val{nreduce},
    ::Val{dim},
    ::Val{N},
    Q,
    target::FilterIndices{I},
    vgeo,
) where {nreduce, dim, N, I}
    @uniform begin
        FT = eltype(Q)

        Nqs = N .+ 1
        Nq1 = Nqs[1]
        Nq2 = dim == 2 ? 1 : Nqs[2]
        Nq3 = Nqs[end]

        nfilterstates = number_state_filtered(target, FT)
        nelemperblock = 1
    end

    l_Q = @private FT (nfilterstates, Nq1)
    l_MJ = @private FT (Nq1,)

    s_MJQ = @localmem FT (Nq1 * Nq2, nfilterstates)
    s_MJQclipped = @localmem FT (Nq1 * Nq2, nfilterstates)

    e = @index(Group, Linear)
    i, j = @index(Local, NTuple)

    @inbounds begin
        # loop up the pencil and load Q and MJ
        @unroll for k in 1:Nq3
            ijk = i + Nq1 * ((j - 1) + Nq2 * (k - 1))

            @unroll for sf in 1:nfilterstates
                s = I[sf]
                l_Q[sf, k] = Q[ijk, s, e]
            end

            l_MJ[k] = vgeo[ijk, _M, e]
        end

        @unroll for sf in 1:nfilterstates
            MJQ, MJQclipped = zero(FT), zero(FT)

            @unroll for k in 1:Nq3
                MJ = l_MJ[k]
                Qs = l_Q[sf, k]
                Qsclipped = Qs ≥ 0 ? Qs : zero(Qs)

                MJQ += MJ * Qs
                MJQclipped += MJ * Qsclipped
            end

            ij = i + Nq1 * (j - 1)

            s_MJQ[ij, sf] = MJQ
            s_MJQclipped[ij, sf] = MJQclipped
        end
        @synchronize

        @unroll for n in 11:-1:1
            if nreduce ≥ 2^n
                ij = i + Nq1 * (j - 1)
                ijshift = ij + 2^(n - 1)
                if ij ≤ 2^(n - 1) && ijshift ≤ Nq1 * Nq2
                    @unroll for sf in 1:nfilterstates
                        s_MJQ[ij, sf] += s_MJQ[ijshift, sf]
                        s_MJQclipped[ij, sf] += s_MJQclipped[ijshift, sf]
                    end
                end
                @synchronize
            end
        end

        @unroll for sf in 1:nfilterstates
            qs_average = s_MJQ[1, sf]
            qs_clipped_average = s_MJQclipped[1, sf]

            r = qs_average > 0 ? qs_average / qs_clipped_average : zero(FT)

            s = I[sf]
            @unroll for k in 1:Nq3
                ijk = i + Nq1 * ((j - 1) + Nq2 * (k - 1))

                Qs = l_Q[sf, k]
                Q[ijk, s, e] = Qs ≥ 0 ? r * Qs : zero(Qs)
            end
        end
    end
end

"""
    kernel_apply_mp_filter!(::Val{dim}, ::Val{N}, direction,
                         Q, state_auxiliary, target, filtermatrix
                        ) where {dim, N}
Computational kernel: Applies the `filtermatrix` to `Q` given a
custom target `target`.
The `direction` argument is used to control if the filter is applied in the
"""
@kernel function kernel_apply_mp_filter!(
    ::Val{nreduce},
    ::Val{dim},
    ::Val{N},
    ::Val{vars_Q},
    ::Val{vars_state_auxiliary},
    direction,
    Q,
    state_auxiliary,
    target::AbstractFilterTarget,
    filtermatrix,
    vgeo,
) where {nreduce, dim, N, vars_Q, vars_state_auxiliary}
    @uniform begin
        FT = eltype(Q)

        Nqs = N .+ 1
        Nq1 = Nqs[1]
        Nq2 = Nqs[2]
        Nq3 = dim == 2 ? 1 : Nqs[dim]

        if direction isa EveryDirection
            filterinξ1 = filterinξ2 = true
            filterinξ3 = dim == 2 ? false : true
        elseif direction isa HorizontalDirection
            filterinξ1 = true
            filterinξ2 = dim == 2 ? false : true
            filterinξ3 = false
        elseif direction isa VerticalDirection
            filterinξ1 = false
            filterinξ2 = dim == 2 ? true : false
            filterinξ3 = dim == 2 ? false : true
        end

        nstates = varsize(vars_Q)
        nfilterstates = number_state_filtered(target, FT)
        nfilteraux =
            isnothing(state_auxiliary) ? 0 : varsize(vars_state_auxiliary)

        # ugly workaround around problems with @private
        # hopefully will be soon fixed in KA
        u_Q = MVector{nstates, FT}(undef)
        u_Qfiltered = MVector{nfilterstates, FT}(undef)
    end

    l_Q = @localmem FT (Nq1, Nq2, Nq3, nfilterstates) # element local
    l_MQᴮ = @localmem FT (Nq1 * Nq2 * Nq3, nstates) # before applying filter
    l_MQᴬ = @localmem FT (Nq1 * Nq2 * Nq3, nstates) # after applying filter
    l_M = @localmem FT (Nq1 * Nq2 * Nq3) # local mass matrix

    p_Q = @private FT (nstates,)
    p_Qfiltered = @private FT (nfilterstates,) # scratch space for storing mat mul
    p_aux = @private FT (nfilteraux,)

    e = @index(Group, Linear)
    i, j, k = @index(Local, NTuple)
    ijk = @index(Local, Linear)

    @inbounds begin

        @unroll for s in 1:nstates
            p_Q[s] = Q[ijk, s, e]
        end

        @unroll for s in 1:nfilteraux
            p_aux[s] = state_auxiliary[ijk, s, e]
        end

        # Load mass matrix and pre-filtered mass weighted quantities to shared memory
        l_M[ijk] = vgeo[ijk, _M, e]
        @unroll for s in 1:nstates
            l_MQᴮ[ijk, s] = l_M[ijk] * p_Q[s]
        end

        fill!(u_Qfiltered, -zero(FT))

        compute_filter_argument!(
            target,
            Vars{vars_state_filtered(target, FT)}(u_Qfiltered),
            Vars{vars_Q}(p_Q[:]),
            Vars{vars_state_auxiliary}(p_aux[:]),
        )

        @unroll for fs in 1:nfilterstates
            p_Qfiltered[fs] = zero(FT)
        end

        @unroll for fs in 1:nfilterstates
            l_Q[i, j, k, fs] = u_Qfiltered[fs]
        end

        if filterinξ1
            @synchronize
            @unroll for n in 1:Nq1
                @unroll for fs in 1:nfilterstates
                    p_Qfiltered[fs] += filtermatrix[i, n] * l_Q[n, j, k, fs]
                end
            end

            if filterinξ2 || filterinξ3
                @synchronize
                @unroll for fs in 1:nfilterstates
                    l_Q[i, j, k, fs] = p_Qfiltered[fs]
                    p_Qfiltered[fs] = zero(FT)
                end
            end
        end

        if filterinξ2
            @synchronize
            @unroll for n in 1:Nq2
                @unroll for fs in 1:nfilterstates
                    p_Qfiltered[fs] += filtermatrix[j, n] * l_Q[i, n, k, fs]
                end
            end

            if filterinξ3
                @synchronize
                @unroll for fs in 1:nfilterstates
                    l_Q[i, j, k, fs] = p_Qfiltered[fs]
                    p_Qfiltered[fs] = zero(FT)
                end
            end
        end

        if filterinξ3
            @synchronize
            @unroll for n in 1:Nq3
                @unroll for fs in 1:nfilterstates
                    p_Qfiltered[fs] += filtermatrix[k, n] * l_Q[i, j, n, fs]
                end
            end
        end

        # work around for not being able to `Vars` `@private` arrays
        @unroll for s in 1:nstates
            u_Q[s] = p_Q[s]
        end

        compute_filter_result!(
            target,
            Vars{vars_Q}(u_Q),
            Vars{vars_state_filtered(target, FT)}(p_Qfiltered[:]),
            Vars{vars_state_auxiliary}(p_aux[:]),
        )
        # Store result and post-filtered mass weighted quantities
        @unroll for s in 1:nstates
            p_Q[s] = u_Q[s]
            l_MQᴬ[ijk, s] = l_M[ijk] * p_Q[s]
        end

        @synchronize
        @unroll for n in 11:-1:1
            if nreduce ≥ (1 << n)
                ijkshift = ijk + (1 << (n - 1))
                if ijk ≤ (1 << (n - 1)) && ijkshift ≤ Nq1 * Nq2 * Nq3
                    l_M[ijk] += l_M[ijkshift]
                    @unroll for s in 1:nstates
                        l_MQᴮ[ijk, s] += l_MQᴮ[ijkshift, s]
                        l_MQᴬ[ijk, s] += l_MQᴬ[ijkshift, s]
                    end
                end
                @synchronize
            end
        end

        @synchronize
        M⁻¹ = 1 / l_M[1]
        # Reset the element average and store result
        @unroll for s in 1:nstates
            Q[ijk, s, e] = p_Q[s] + M⁻¹ * (l_MQᴮ[1, s] - l_MQᴬ[1, s])
        end

        @synchronize
    end

end

end # end of module