text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

function parallelize_batch_tile_loops(statements)
    new_body = []
    for batch_loop in statements
        if isa(batch_loop, Expr) && batch_loop.head == :for
            batch_body = batch_loop.args[2]
            batch_body.args = fuse_loops(batch_body.args)
            batch_range = batch_loop.args[1].args[2]
            num_micro_batches = cld(batch_range.args[2], MICRO_BATCH_SIZE)
            batch_loopvar = batch_loop.args[1].args[1]
            statement = batch_body.args[1]
            tile_loops = []
            AstWalk(batch_loop, get_tile_loops, tile_loops)
            loopvars = [batch_loopvar]
            body = batch_body
            # batch_range = :((micro_batch - 1) * $MICRO_BATCH_SIZE + 1:
            #                 min(micro_batch * $MICRO_BATCH_SIZE, $(batch_range.args[2])))
            ranges = [batch_range]
            if length(tile_loops) > 0 && length(batch_body.args) == 1
                for loop in tile_loops
                    push!(ranges, loop.args[1].args[2])
                    push!(loopvars, loop.args[1].args[1])
                end
                body = tile_loops[end].args[2]
            end
            pre = []
            for var in loopvars
                push!(pre, :($var = _REMOVE_THIS_LINE()))
            end
            unique = gensym("loop_ranges")
            start_vars = Union{Symbol,Int}[]
            stop_vars = Union{Symbol,Int}[]
            for (i, range) in enumerate(ranges)
                # range_var = symbol(unique,:range,i)
                # push!(pre, :($range_var = $range))
                start_var = symbol(unique,:start,i)
                push!(pre, :($start_var = $(range.args[1])))
                push!(start_vars, start_var)
                stop_var = symbol(unique,:stop,i)
                push!(pre, :($stop_var = $(range.args[2])))
                push!(stop_vars, stop_var)
            end
            push!(new_body, quote# for micro_batch = 1:$num_micro_batches
                $(pre...)
                $(mk_parallel_loophead(loopvars, start_vars, stop_vars;schedule="schedule(static,1)"))
                $(body.args...)
                $(Expr(:parallel_loopend, length(loopvars)))
                # parallel_for($(ranges...)) do $(loopvars...)
                #     $body
                # end
            # end))
            end)
        else
            push!(new_body, batch_loop)
        end
    end
    new_body
end

function collect_parallel_loop_private_vars(expr)
    function private_var_collector(node, cbdata, index, top_level, read)
        if cbdata[1] != nothing
            if isa(node, Expr) && node.head == :(=)
                lhs = node.args[1]
                if isa(lhs, Symbol)
                    push!(cbdata[1], lhs)
                elseif isa(lhs, SymbolNode)
                    push!(cbdata[1], lhs.name)
                elseif isa(lhs, GenSym)
                    push!(cbdata[1], lhs)
                end
                return node
            elseif isa(node, Expr) && node.head == :loophead
                private_var_collector(node.args[1], cbdata, index, top_level, read)
                return node
            end
        end
        if isa(node, Expr) && node.head == :parallel_loophead
            cbdata[1] = node.args[4]
            return node
        elseif isa(node, Expr) && node.head == :parallel_loopend
            cbdata[1] = nothing
            return node
        elseif isa(node, Expr) && node.head in [:loopend, :loophead]
            return node
        end
        ASTWALK_RECURSE
    end
    AstWalk(expr, private_var_collector, Any[nothing])
end


text:
include("sampler_test_utils.jl")

describe("Truncate to Support Sampler") do

	describe("construction") do

		subsA = DataGenerators.NormalSampler([102.0, 45.81])
		s = DataGenerators.TruncateToSupportSampler(subsA)

		test("numparams and paramranges") do
			@check DataGenerators.numparams(s) == DataGenerators.numparams(subsA)
			@check DataGenerators.paramranges(s) == DataGenerators.paramranges(subsA)
		end

		test("default params") do
			@check DataGenerators.getparams(s) == DataGenerators.getparams(subsA)
		end

		test("set params") do
			DataGenerators.setparams!(s, [-9.8, 0.002])
			@check DataGenerators.getparams(subsA) == [-9.8, 0.002]
			@check DataGenerators.getparams(s) == [-9.8, 0.002]
		end

	end

	describe("sampling from continuous with infinite support") do

		subsA = DataGenerators.NormalSampler([42.8, 98.77])
		s = DataGenerators.TruncateToSupportSampler(subsA)

		@repeat test("random support") do
			bounds = rand(2)*40.0-20.0
			l, u = minimum(bounds), maximum(bounds)
			x, trace = DataGenerators.sample(s, (l, u))
			@check typeof(x) <: Float64
			@check l <= x <= u
		end

		test("equal support") do
			bound = rand()*40.0-20.0
			for i in 1:100
				x, trace = DataGenerators.sample(s, (bound, bound))
				@check x == bound
			end
		end

		@repeat test("entire float range") do
			x, trace = DataGenerators.sample(s, (-realmax(Float64), realmax(Float64)))
			@check typeof(x) <: Float64
			@check isfinite(x)
			@mcheck_that_sometimes x < 0.0
			@mcheck_that_sometimes x > 0.0
		end

	end

	describe("sampling from continuous with finite support") do

		subsA = DataGenerators.UniformSampler([42.5, 79.67])
		s = DataGenerators.TruncateToSupportSampler(subsA)

		@repeat test("random support") do
			bounds = rand(2)*100.0
			l, u = minimum(bounds), maximum(bounds)
			x, trace = DataGenerators.sample(s, (l, u))
			@check typeof(x) <: Float64
			@check min(max(l,42.5),79.67) <= x <= max(min(u,79.67),42.5) # outer max/min because this is the fallback made by distributions.jl when the truncation and distribution support do not overlap
		end

		test("equal support") do
			bound = 51.8
			for i in 1:100
				x, trace = DataGenerators.sample(s, (float64(bound), float64(bound)))
				@check x == bound
			end
		end

		@repeat test("entire float range") do
			x, trace = DataGenerators.sample(s, (-realmax(Float64), realmax(Float64)))
			@check typeof(x) <: Float64
			@check isfinite(x)
			@check 42.5 <= x <= 79.67
		end

	end

	describe("estimate parameters") do

		test("estimates parameters of subsampler (when not too constrained)") do

			subs1A = DataGenerators.NormalSampler()
			s1 = DataGenerators.TruncateToSupportSampler(subs1A)
			params = [94.2, 34.1]
			DataGenerators.setparams!(s1, params)

			subs2A = DataGenerators.NormalSampler()
			s2 = DataGenerators.TruncateToSupportSampler(subs2A)
			otherparams = [-42.0, 50.1]
			DataGenerators.setparams!(s2, otherparams)

			traces = map(1:100) do i
			 	# note: large support that rarely truncates, otherwise would not be able to estimate subsampler with same params
				x, trace = DataGenerators.sample(s1, (-500.0, 500.0))
				trace
			end

			estimateparams!(s2, traces)

			@check isconsistentnormal(subs2A, params[1:2])

		end

	end

end

ext: jl


text:
## --- test ArrayStats.jl

    # Sorting, counting, matching
    A = rand(1:100.,100); B = sort(A)
    @test A[1:count_unique!(A)] == unique(B)
    @test findclosest(3.3:5.3,1:10) == 3:5
    @test findclosestbelow(3.3:5.3,1:10) == 3:5
    @test findclosestabove(3.3:5.3,1:10) == 4:6
    @test findmatches(40:60,1:100) == 40:60
    @test findnth(fill(true,50), 25) == 25

    # Interpolation
    @test cntr(0:2:100) == 1:2:99

    # Integration
    @test trapezoidalquadrature(1:10,fill(1,10)) == 9
    @test trapz(collect(1:10.),ones(10)) == 9
    @test midpointquadrature(1:10,ones(10)) == 10

    # Distributions
    A = draw_from_distribution(ones(100), 10000)::AbstractArray
    @test length(A) == 10000
    @test isapprox(nanmean(A), 0.5, atol=0.08)
    @test isapprox(nanstd(A), sqrt(1/12), atol=0.08)

    # Strings
    @test contains("JuliaLang is pretty cool!", "Julia")
    @test !contains("JuliaLang is pretty cool!", "julia")

    @test containsi("JuliaLang is pretty cool!", "Julia")
    @test containsi("JuliaLang is pretty cool!", "julia")
    @test !containsi("JuliaLang is pretty cool!", "tomatoes")

## ---

max_stars_repo_name: brenhinkeller/StatGeochemBase

max_issues_repo_name: brenhinkeller/StatGeochemBase

max_forks_repo_name: brenhinkeller/StatGeochemBase


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/NodeNumbering.jl/blob/master/LICENSE

using NodeNumbering

@testset "create adjacency graph" begin
    expected = Dict(
        1 => [3, 8, 9],
        2 => [3, 8, 7, 6],
        3 => [1, 2],
        4 => [8, 9],
        5 => [7, 8],
        6 => [2, 7],
        7 => [5, 2, 6],
        8 => [1, 2, 4, 5],
        9 => [1, 4])

       elements=Dict(
            1 => [9, 1, 8, 4],
            2 => [1, 3, 2, 8],
            3 => [8, 2, 7, 5],
            4 => [2, 6, 7])

        element_types = Dict{Int, Symbol}(
            1 => :Quad4,
            2 => :Quad4,
            3 => :Quad4,
            4 => :Tri3)

    result = NodeNumbering.create_adjacency_graph(elements, element_types)
    result_bool = true
    for (k, v) in expected
        if (length(intersect(v, result[k]))==length(v))==false
            result_bool = false
            break
        end
    end
    #@test result == expected
    @test result_bool==true
end


text:
import PackageCompilerX

@inline function default_project()
    return "@v$(Base.VERSION.major).$(Base.VERSION.minor)"
end

@inline function fast_and_simple(project::AbstractString = _default_project();
                                 sysimage_path::Union{AbstractString, Nothing} = nothing,
                                 notest::AbstractVector{<:Symbol} = Symbol[],
                                 replace_default::Bool = false,
                                 shared_project::Bool = false,
                                 test_stdlibs::Bool = false,
                                 force_replace_default::Bool = false)
    if replace_default && project != default_project()
        msg = "`replace_default` is `true` but `project` is not the default project (`$(default_project())`)"
        if force_replace_default
            @warn(msg)
        else
            throw(ArgumentError(msg))
        end
    end
    _with_different_project(project; shared_project = shared_project) do activated_project
        _add_all_project_recursive_dependencies()
        project_recursive_dependencies = unique(Symbol.(_get_project_recursive_dependencies()))
        precompile_execution = _precompile_execution(project_recursive_dependencies;
                                                     notest = notest,
                                                     test_stdlibs = test_stdlibs)
        PackageCompilerX.create_sysimage(project_recursive_dependencies;
                                         project = dirname(activated_project),
                                         replace_default = replace_default,
                                         sysimage_path = sysimage_path,
                                         precompile_execution_file=precompile_execution)
        return nothing
    end
    return nothing
end


text:
module SLMC
    import ..Actions:GaugeActionParam_standard,
                    GaugeActionParam,
                    GaugeActionParam_autogenerator

    using LinearAlgebra
    mutable struct SLMC_data
        KdagK::Array{Float64,2}
        KdagY::Array{Float64,1}
        factor::Float64

        function SLMC_data(n,NC)
            KdagK = zeros(Float64,n+1,n+1)
            KdagY= zeros(Float64,n+1)
            factor = -1/NC
            return new(KdagK,KdagY,factor)
        end
    end

    function update_slmcdata!(s::SLMC_data,inputdata,outputdata)
        n = length(s.KdagY)
        KdagK = zero(s.KdagK)
        KdagY = zero(s.KdagY)
        S = outputdata

        K = zeros(Float64,n)

        K[1] = 1
        for i=2:n
            K[i] = s.factor*inputdata[i-1]
        end
        println("K = ",K)

        s.KdagY += K[:]*S #KdagK
        s.KdagK += K*K' #KdagY

        #numdata = s.KdagK[1,1]
        #s.KdagY ./= numdata
        #s.KdagK ./= numdata
        return
    end

    function make_beta(s::SLMC_data,num)
        #if num == 1
            #detKdagK = det(s.KdagK[1:num+1,1:num+1])
            #println("det",detKdagK)
            #invKdagK = [s.KdagK[2,2] -s.KdagK[1,2]
            #s.KdagK[2,1] s.KdagK[1,1]
            #]/detKdagK
            #println("invKdagK ",invKdagK)
            #println("KdagY ",s.KdagY[1:num+1] )

        #    slmc_betas = invKdagK*s.KdagY[1:num+1]
        #else
            slmc_betas = s.KdagK[1:num+1,1:num+1] \ s.KdagY[1:num+1]
        #end
        return slmc_betas
    end

    function show_effbeta(s::SLMC_data,effectiveterms)
        #effectiveterms = ["plaq"]
        n = length(s.KdagY)
        be1=zeros(n)
        IsSucs=true

        try
            for i=1:n-1
                print("#Estimation$i learning ")
                print("$(effectiveterms[1:i]) only: ")
                be1[1:i+1] = make_beta(s,i)
                #println("slmc_betas ",be1 )
                for k = 1:length(be1)-1
                    print("$(be1[k+1]) \t")
                end
                println("constant is ",be1[1])
            end
        catch
            println("Failed to calculate on-fly SLMC parameter.")
            println(s.KdagK,"\t",s.KdagY)

            IsSucs=false
        end
        return be1[2:end],be1[1],IsSucs
    end

    function show_effbeta(s::SLMC_data)
        effectiveterms = ["plaq"]
        return show_effbeta(s,effectiveterms)
    end

    function show_effbeta(s::SLMC_data,gparam::GaugeActionParam)

        if typeof(gparam) == GaugeActionParam_autogenerator
            effectiveterms = gparam.couplinglist
            println(gparam.couplinglist )
        elseif typeof(gparam) == GaugeActionParam_standard
            effectiveterms = ["plaq"]
        else
            error("error!")
        end
        return show_effbeta(s,effectiveterms)
    end

end


text:
function readinput(filename)
    (rules, messages) = split(read(filename, String), "\n\n")
    ruledict = Dict()
    for group in eachmatch(r"(\d+): (.*)", rules)
        ruledict[group[1]] = group[2]
    end
    return ruledict, split(messages)
end

function parserule(ruledict, target)
    targetrules = ruledict[target]
    if occursin(r"[ab]", targetrules)
        return strip(targetrules, '\"')
    end
    for rule in eachmatch(r"(\d+)", targetrules)
        targetrules = replace(
            targetrules,
            rule[1] => parserule(ruledict, rule[1]);
            count = 1
        )
    end
    if occursin("|", targetrules)
        return "($targetrules)"
    else
        return targetrules
    end
end

function makeregex(rules, target)
    rule = replace(parserule(rules, target), " " => "")
    Regex("^$rule\$")
end

(rules, messages) = readinput("d19.in")
rule = makeregex(rules, "0")
println("P1: ", count(i->!isnothing(match(rule, i)), messages))

function repeat(ruledict, target)
    rule = ruledict[target]
    template = rule
    for i in 1:2
        rule = replace(rule, target => "($rule)")
    end
    replace(rule, "($template)"=>split(template, "|")[1])
end

rules["8"] = "42 | 42 8" # original
rules["8"] = "42 | 42 (42)+" # no self-reference
rules["11"] = "42 31 | 42 11 31" # original
rules["11"] = repeat(rules, "11") # no self-reference

rule = makeregex(rules, "0")
println("P2: ", count(i->!isnothing(match(rule, i)), messages))

max_stars_repo_name: eduellery/adventofcode

max_issues_repo_name: eduellery/adventofcode

max_forks_repo_name: eduellery/adventofcode


text:
using BCBIStats.COOccur
using StatsBase

#Co-occurrence matrix
occur_file = "/Users/isa/Dropbox/BrownAgain/Projects/BCBI/test/occur_sp.jdl"
labels_file = "/Users/isa/Dropbox/BrownAgain/Projects/BCBI/test/mesh2ind.jdl"

#Contingency
x=[1,0,0,1,0,1,1,0,0,1,0,1,1,1,1,1,1,1,0,0,0,0,1,0,1,1,0,0,1,0,1]
y=[1,1,1,0,0,0,0,1,1,0,0,1,0,1,1,1,1,0,0,1,1,1,1,1,1,1,1,0,0,0,0]

occur = [x y]
label2ind = Dict("x" => 1, "y" => 2)

try
  global occur = read_sp_occur(occur_file)
  global label2ind = read_occur_label_dict(labels_file)
catch
  nothing
end

COO = occur2coo(occur, label2ind)

corr = corrcoef(occur)
pmi = pmi_mat(COO.coo_matrix, sum(occur))

d = counts(x,y, span(occur))
stat = chi2_statistic(x,y,span(occur))
@test isapprox(stat,  0.4057853910795092)
stat = chi2_statistic(d)
@test isapprox(stat, 0.4057853910795092)
stat = chi2_statistic(d, min_freq=6)
@test isnan(stat)

max_stars_repo_name: bcbi/BCBIStats

max_issues_repo_name: bcbi/BCBIStats

max_forks_repo_name: bcbi/BCBIStats


text:
const BACKSPACE          = UInt8('\b')
const TAB                = UInt8('\t')
const NEWLINE            = UInt8('\n')
const FORM               = UInt8('\f')
const RETURN             = UInt8('\r')
const SPACE              = UInt8(' ')
const QUOTE              = UInt8('"')
const PLUS               = UInt8('+')
const COMMA              = UInt8(',')
const MINUS              = UInt8('-')
const PERIOD             = UInt8('.')
const FORWARDSLASH       = UInt8('/')
const BACKSLASH          = UInt8('\\')
const ZERO               = UInt8('0')
const NEG_ONE            = UInt8('0') - 0x01
const TEN                = UInt8('9') + 0x01
const BIG_A              = UInt8('A')
const BIG_E              = UInt8('E')
const BIG_F              = UInt8('F')
const BIG_I              = UInt8('I')
const BIG_N              = UInt8('N')
const BIG_T              = UInt8('T')
const BIG_Y              = UInt8('Y')
const LITTLE_A           = UInt8('a')
const LITTLE_B           = UInt8('b')
const LITTLE_E           = UInt8('e')
const LITTLE_F           = UInt8('f')
const LITTLE_I           = UInt8('i')
const LITTLE_L           = UInt8('l')
const LITTLE_N           = UInt8('n')
const LITTLE_R           = UInt8('r')
const LITTLE_S           = UInt8('s')
const LITTLE_T           = UInt8('t')
const LITTLE_U           = UInt8('u')
const LITTLE_Y           = UInt8('y')
const OPEN_CURLY_BRACE   = UInt8('{')
const CLOSE_CURLY_BRACE  = UInt8('}')
const COLON              = UInt8(':')
const OPEN_SQUARE_BRACE  = UInt8('[')
const CLOSE_SQUARE_BRACE = UInt8(']')

function escapechar(b)
    b == QUOTE     && return QUOTE
    b == BACKSLASH && return BACKSLASH
    b == BACKSPACE && return LITTLE_B
    b == FORM      && return LITTLE_F
    b == NEWLINE   && return LITTLE_N
    b == RETURN    && return LITTLE_R
    b == TAB       && return LITTLE_T
    return 0x00
end

function reverseescapechar(b)
    b == QUOTE     && return QUOTE
    b == BACKSLASH && return BACKSLASH
    b == LITTLE_B  && return BACKSPACE
    b == LITTLE_F  && return FORM
    b == LITTLE_N  && return NEWLINE
    b == LITTLE_R  && return RETURN
    b == LITTLE_T  && return TAB
    return 0x00
end

charvalue(b) = (NEG_ONE < b < TEN)         ? b - ZERO            :
               (LITTLE_A <= b <= LITTLE_F) ? b - LITTLE_A - 0x0a :
               (BIG_A <= b <= BIG_F)       ? b - BIG_A - 0x0a    :
               throw(ArgumentError("JSON invalid unicode hex value"))

function readhexchar(bytes, i)
    n = 0x0000
    foreach(1:4) do _
        @inbounds b = bytes[i]
        i += 1
        n <<= 4
        n += charvalue(b)
    end
    return Char(n), i
end

iscntrl(c::Char) = c <= '\x1f' || '\x7f' <= c <= '\u9f'
function escaped(b)
    if b == FORWARDSLASH
        return [FORWARDSLASH]
    elseif b >= 0x80
        return [b]
    elseif b in map(UInt8, ('"', '\\', '\b', '\f', '\n', '\r', '\t'))
        return [BACKSLASH, escapechar(b)]
    elseif iscntrl(Char(b))
        return UInt8[BACKSLASH, LITTLE_U, string(b, base=16, pad=4)...]
    else
        return [b]
    end
end

const ESCAPECHARS = [escaped(b) for b = 0x00:0xff]
const NEEDESCAPE = [length(x) > 1 for x in ESCAPECHARS]

function needescape(str)
    bytes = codeunits(str)
    @simd for i = 1:length(bytes)
        @inbounds NEEDESCAPE[bytes[i] + 0x01] && return true
    end
    return false
end

function write(io::IO, obj::AbstractString)
    Base.write(io, '"')
    if needescape(obj)
        bytes = codeunits(obj)
        for i = 1:length(bytes)
            @inbounds b = JSON2.ESCAPECHARS[bytes[i] + 0x01]
            Base.write(io, b)
        end
    else
        Base.write(io, obj)
    end
    Base.write(io, '"')
    return
end

# read
const BACKSLASH = UInt8('\\')

function unescape(str)
    bytes = codeunits(str)
    # skip the opening/closing '"'
    i = 1
    len = length(bytes)
    while i <= len
        @inbounds b = bytes[i]
        i += 1
        if b == BACKSLASH
            @inbounds b = bytes[i]
            i += 1
            if b == LITTLE_U
                c, i = readhexchar(bytes, i)
                Base.write(BUF, c)
            else
                d = reverseescapechar(b)
                d == 0x00 && throw(ArgumentError("JSON invalid escape character: '$(Char(b))'"))
                Base.write(BUF, d)
            end
            continue
        elseif b < SPACE
            throw(ArgumentError("JSON encountered unescaped control character: '$(Char(b))'"))
        end
        Base.write(BUF, b)
    end
    return String(take!(BUF))
end

function read(io::IO, T::Type{String})
    @expect '"'
    b = readbyte(io)
    hasescapechars = false
    while b != UInt8('"')
        Base.write(BUF, b)
        if b == UInt8('\\')
            eof(io) && throw(ArgumentError("early EOF"))
            hasescapechars = true
            b = readbyte(io)
            Base.write(BUF, b)
        end
        eof(io) && throw(ArgumentError("early EOF"))
        b = readbyte(io)
    end
    str = String(take!(BUF))
    return hasescapechars ? unescape(str) : str
end

function read(io::IOBuffer, T::Type{String})
    JSON2.@expect '"'
    ptr = pointer(io.data, io.ptr)
    b = readbyte(io)
    len = 0
    hasescapechars = false
    while b != JSON2.QUOTE
        len += 1
        if b == JSON2.BACKSLASH
            eof(io) && throw(ArgumentError("early EOF"))
            hasescapechars = true
            b = readbyte(io)
            len += 1
        end
        eof(io) && throw(ArgumentError("early EOF"))
        b = readbyte(io)
    end
    str = unsafe_string(ptr, len)
    return hasescapechars ? unescape(str) : str
end

function read(io::IOBuffer, T::Type{Symbol})
    JSON2.@expect '"'
    ptr = pointer(io.data, io.ptr)
    b = readbyte(io)
    len = 0
    hasescapechars = false
    while b != JSON2.QUOTE
        len += 1
        if b == JSON2.BACKSLASH
            eof(io) && throw(ArgumentError("early EOF"))
            hasescapechars = true
            b = readbyte(io)
            len += 1
        end
        eof(io) && throw(ArgumentError("early EOF"))
        b = readbyte(io)
    end
    sym = ccall(:jl_symbol_n, Ref{Symbol}, (Ptr{UInt8}, Int), ptr, len)
    return hasescapechars ? Symbol(unescape(unsafe_string(ptr, len))) : sym
end
# slow path
read(io::IO, ::Type{Symbol}) = Symbol(read(io, String))


text:
module BoundaryConditions

export
    BCType, Flux, Gradient, Value, NoPenetration,
    BoundaryCondition, bctype, getbc, setbc!,
    PeriodicBoundaryCondition, NoPenetrationBoundaryCondition, NoFluxBoundaryCondition,
    FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition,
    CoordinateBoundaryConditions,
    FieldBoundaryConditions, UVelocityBoundaryConditions, VVelocityBoundaryConditions,
    WVelocityBoundaryConditions, TracerBoundaryConditions, PressureBoundaryConditions,
    DiffusivityBoundaryConditions,
    BoundaryFunction,
    apply_z_bcs!, apply_y_bcs!,
    fill_halo_regions!, zero_halo_regions!

using CUDAnative

using Oceananigans: Cell, Face
using Oceananigans.Architectures
using Oceananigans.Grids

include("boundary_condition_types.jl")
include("boundary_condition.jl")
include("coordinate_boundary_conditions.jl")
include("field_boundary_conditions.jl")
include("boundary_function.jl")
include("show_boundary_conditions.jl")

include("fill_halo_regions.jl")
include("zero_halo_regions.jl")

include("apply_flux_bcs.jl")
include("apply_value_gradient_bcs.jl")

end


text:
"""
    is_option(x)

Check if `x` is an option type or not.
"""
is_option(x) = false

"""
    compare_options(a, b, xs...)::Bool

Compare option types check if they are the same.
"""
function compare_options(a, b, xs...)::Bool
    compare_options(a, b) || return false
    compare_options(b, xs...)
end

compare_options(a, b) = false

function compare_options(a::A, b::A) where {A}
    for each in fieldnames(A)
        getfield(a, each) == getfield(b, each) || return false
    end
    return true
end


text:
mutable struct sem_mesh_data
    nspec::Int64
    nglob::Int64
    xyz_glob::Array{Float64,2}
    ibool::Array{Int64,4}
    idoubling::Array{Int64,1}
    ispec_is_tiso::Array{Bool,1}

    # xix::Array{Float64,4}
    # xiy::Array{Float64,4}
    # xiz::Array{Float64,4}
    # etax::Array{Float64,4}
    # etay::Array{Float64,4}
    # etaz::Array{Float64,4}
    # gammax::Array{Float64,4}
    # gammay::Array{Float64,4}
    # gammaz::Array{Float64,4}

    sem_mesh_data() = new()
end

mutable struct sem_mesh_location
    stat::Int64
    eid::Int64

    uvw::Vector{Float64}
    misloc::Float64
    lagrange::Array{Float64,3}

    sem_mesh_location() = new()
end

mutable struct sem_boundary_disc
    NSPEC2D_MOHO::Int64
    NSPEC2D_400::Int64
    NSPEC2D_670::Int64
    ibelm_moho_top::Array{Int64,1}
    ibelm_moho_bot::Array{Int64,1}
    ibelm_400_top::Array{Int64,1}
    ibelm_400_bot::Array{Int64,1}
    ibelm_670_top::Array{Int64,1}
    ibelm_670_bot::Array{Int64,1}

    sem_boundary_disc() = new()
end

max_stars_repo_name: ziyixi/seisflow

max_issues_repo_name: ziyixi/seisflow

max_forks_repo_name: ziyixi/seisflow


text:
include("mp-datastructures.jl")
include("mp-functions.jl")
using Plots
pyplot()

"""    probabilities(ngames, nmoves)
returns percentage (%) probabilities of landing on each square"""
function probabilities(ngames, nmoves)
    countsquares = fill(0, 40)
    for game in 1:ngames
        mygame = newgame(1500)
        for move in 1:nmoves
            visited = taketurn!(mygame)
            # mygame.currentplayer = 1 + (mygame.currentplayer % length(mygame.players))
            for i in visited
                countsquares[i] += 1
            end
        end
    end
    probs = 100 * countsquares ./ sum(countsquares)
    return probs
end

@time probs = probabilities(100_000, 80)
x = [i.name for i in origsquares]
fig = plot(probs
    ,xticks=(1:40, x)
    ,seriestype=:bar, legend=false, size=(900,600)
    ,title="Probabilities of Finishing on Square"
    ,xrotation=90, xtickfontsize=9
    ,ytickfontsize=10, gridalpha=.5
)


text:
const PVCAM_LIB_VERSION_MAJOR = 2
const MAX_CAM = 16
const CAM_NAME_LEN = 32
const PARAM_NAME_LEN = 32
const ERROR_MSG_LEN = 255
const CCD_NAME_LEN = 17
const MAX_ALPHA_SER_NUM_LEN = 32
const MAX_PP_NAME_LEN = 32
const MAX_SYSTEM_NAME_LEN = 32
const MAX_VENDOR_NAME_LEN = 32
const MAX_PRODUCT_NAME_LEN = 32
const MAX_CAM_PART_NUM_LEN = 32
const MAX_GAIN_NAME_LEN = 32
const MAX_CAM_SYSTEMS_INFO_LEN = 1024
const PP_MAX_PARAMETERS_PER_FEATURE = 10
const PL_MD_FRAME_SIGNATURE = 5328208
const PL_MD_EXT_TAGS_MAX_SUPPORTED = 255

const TYPE_INT16 = 1
const TYPE_INT32 = 2
const TYPE_FLT64 = 4
const TYPE_UNS8 = 5
const TYPE_UNS16 = 6
const TYPE_UNS32 = 7
const TYPE_UNS64 = 8
const TYPE_ENUM = 9
const TYPE_BOOLEAN = 11
const TYPE_INT8 = 12
const TYPE_CHAR_PTR = 13
const TYPE_VOID_PTR = 14
const TYPE_VOID_PTR_PTR = 15
const TYPE_INT64 = 16
const TYPE_SMART_STREAM_TYPE = 17
const TYPE_SMART_STREAM_TYPE_PTR = 18
const TYPE_FLT32 = 19
const CLASS0 = 0
const CLASS2 = 2
const CLASS3 = 3

const PARAM_DD_INFO_LENGTH = ( ( CLASS0 << 16 ) + ( TYPE_INT16 << 24 ) + 1 )
const PARAM_DD_VERSION = ( ( CLASS0 << 16 ) + ( TYPE_UNS16 << 24 ) + 2 )
const PARAM_DD_RETRIES = ( ( CLASS0 << 16 ) + ( TYPE_UNS16 << 24 ) + 3 )
const PARAM_DD_TIMEOUT = ( ( CLASS0 << 16 ) + ( TYPE_UNS16 << 24 ) + 4 )
const PARAM_DD_INFO = ( ( CLASS0 << 16 ) + ( TYPE_CHAR_PTR << 24 ) + 5 )
const PARAM_CAM_INTERFACE_TYPE = ( ( CLASS0 << 16 ) + ( TYPE_ENUM << 24 ) + 10 )
const PARAM_CAM_INTERFACE_MODE = ( ( CLASS0 << 16 ) + ( TYPE_ENUM << 24 ) + 11 )
const PARAM_ADC_OFFSET = ( ( CLASS2 << 16 ) + ( TYPE_INT16 << 24 ) + 195 )
const PARAM_CHIP_NAME = ( ( CLASS2 << 16 ) + ( TYPE_CHAR_PTR << 24 ) + 129 )
const PARAM_SYSTEM_NAME = ( ( CLASS2 << 16 ) + ( TYPE_CHAR_PTR << 24 ) + 130 )
const PARAM_VENDOR_NAME = ( ( CLASS2 << 16 ) + ( TYPE_CHAR_PTR << 24 ) + 131 )
const PARAM_PRODUCT_NAME = ( ( CLASS2 << 16 ) + ( TYPE_CHAR_PTR << 24 ) + 132 )
const PARAM_CAMERA_PART_NUMBER = ( ( CLASS2 << 16 ) + ( TYPE_CHAR_PTR << 24 ) + 133 )
const PARAM_COOLING_MODE = ( ( CLASS2 << 16 ) + ( TYPE_ENUM << 24 ) + 214 )
const PARAM_PREAMP_DELAY = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 502 )
const PARAM_COLOR_MODE = ( ( CLASS2 << 16 ) + ( TYPE_ENUM << 24 ) + 504 )
const PARAM_MPP_CAPABLE = ( ( CLASS2 << 16 ) + ( TYPE_ENUM << 24 ) + 224 )
const PARAM_PREAMP_OFF_CONTROL = ( ( CLASS2 << 16 ) + ( TYPE_UNS32 << 24 ) + 507 )
const PARAM_PREMASK = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 53 )
const PARAM_PRESCAN = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 55 )
const PARAM_POSTMASK = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 54 )
const PARAM_POSTSCAN = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 56 )
const PARAM_PIX_PAR_DIST = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 500 )
const PARAM_PIX_PAR_SIZE = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 63 )
const PARAM_PIX_SER_DIST = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 501 )
const PARAM_PIX_SER_SIZE = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 62 )
const PARAM_SUMMING_WELL = ( ( CLASS2 << 16 ) + ( TYPE_BOOLEAN << 24 ) + 505 )
const PARAM_FWELL_CAPACITY = ( ( CLASS2 << 16 ) + ( TYPE_UNS32 << 24 ) + 506 )
const PARAM_PAR_SIZE = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 57 )
const PARAM_SER_SIZE = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 58 )
const PARAM_ACCUM_CAPABLE = ( ( CLASS2 << 16 ) + ( TYPE_BOOLEAN << 24 ) + 538 )
const PARAM_FLASH_DWNLD_CAPABLE = ( ( CLASS2 << 16 ) + ( TYPE_BOOLEAN << 24 ) + 539 )
const PARAM_READOUT_TIME = ( ( CLASS2 << 16 ) + ( TYPE_FLT64 << 24 ) + 179 )
const PARAM_CLEARING_TIME = ( ( CLASS2 << 16 ) + ( TYPE_INT64 << 24 ) + 180 )
const PARAM_POST_TRIGGER_DELAY = ( ( CLASS2 << 16 ) + ( TYPE_INT64 << 24 ) + 181 )
const PARAM_PRE_TRIGGER_DELAY = ( ( CLASS2 << 16 ) + ( TYPE_INT64 << 24 ) + 182 )
const PARAM_CLEAR_CYCLES = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 97 )
const PARAM_CLEAR_MODE = ( ( CLASS2 << 16 ) + ( TYPE_ENUM << 24 ) + 523 )
const PARAM_FRAME_CAPABLE = ( ( CLASS2 << 16 ) + ( TYPE_BOOLEAN << 24 ) + 509 )
const PARAM_PMODE = ( ( CLASS2 << 16 ) + ( TYPE_ENUM << 24 ) + 524 )
const PARAM_TEMP = ( ( CLASS2 << 16 ) + ( TYPE_INT16 << 24 ) + 525 )
const PARAM_TEMP_SETPOINT = ( ( CLASS2 << 16 ) + ( TYPE_INT16 << 24 ) + 526 )
const PARAM_CAM_FW_VERSION = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 532 )
const PARAM_HEAD_SER_NUM_ALPHA = ( ( CLASS2 << 16 ) + ( TYPE_CHAR_PTR << 24 ) + 533 )
const PARAM_PCI_FW_VERSION = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 534 )
const PARAM_FAN_SPEED_SETPOINT = ( ( CLASS2 << 16 ) + ( TYPE_ENUM << 24 ) + 710 )
const PARAM_CAM_SYSTEMS_INFO = ( ( CLASS2 << 16 ) + ( TYPE_CHAR_PTR << 24 ) + 536 )
const PARAM_EXPOSURE_MODE = ( ( CLASS2 << 16 ) + ( TYPE_ENUM << 24 ) + 535 )
const PARAM_EXPOSE_OUT_MODE = ( ( CLASS2 << 16 ) + ( TYPE_ENUM << 24 ) + 560 )
const PARAM_BIT_DEPTH = ( ( CLASS2 << 16 ) + ( TYPE_INT16 << 24 ) + 511 )
const PARAM_IMAGE_FORMAT = ( ( CLASS2 << 16 ) + ( TYPE_ENUM << 24 ) + 248 )
const PARAM_IMAGE_COMPRESSION = ( ( CLASS2 << 16 ) + ( TYPE_ENUM << 24 ) + 249 )
const PARAM_SCAN_MODE = ( ( CLASS3 << 16 ) + ( TYPE_ENUM << 24 ) + 250 )
const PARAM_SCAN_DIRECTION = ( ( CLASS3 << 16 ) + ( TYPE_ENUM << 24 ) + 251 )
const PARAM_SCAN_DIRECTION_RESET = ( ( CLASS3 << 16 ) + ( TYPE_BOOLEAN << 24 ) + 252 )
const PARAM_SCAN_LINE_DELAY = ( ( CLASS3 << 16 ) + ( TYPE_UNS16 << 24 ) + 253 )
const PARAM_SCAN_LINE_TIME = ( ( CLASS3 << 16 ) + ( TYPE_INT64 << 24 ) + 254 )
const PARAM_SCAN_WIDTH = ( ( CLASS3 << 16 ) + ( TYPE_UNS16 << 24 ) + 255 )
const PARAM_GAIN_INDEX = ( ( CLASS2 << 16 ) + ( TYPE_INT16 << 24 ) + 512 )
const PARAM_SPDTAB_INDEX = ( ( CLASS2 << 16 ) + ( TYPE_INT16 << 24 ) + 513 )
const PARAM_GAIN_NAME = ( ( CLASS2 << 16 ) + ( TYPE_CHAR_PTR << 24 ) + 514 )
const PARAM_READOUT_PORT = ( ( CLASS2 << 16 ) + ( TYPE_ENUM << 24 ) + 247 )
const PARAM_PIX_TIME = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 516 )
const PARAM_SHTR_CLOSE_DELAY = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 519 )
const PARAM_SHTR_OPEN_DELAY = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 520 )
const PARAM_SHTR_OPEN_MODE = ( ( CLASS2 << 16 ) + ( TYPE_ENUM << 24 ) + 521 )
const PARAM_SHTR_STATUS = ( ( CLASS2 << 16 ) + ( TYPE_ENUM << 24 ) + 522 )
const PARAM_IO_ADDR = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 527 )
const PARAM_IO_TYPE = ( ( CLASS2 << 16 ) + ( TYPE_ENUM << 24 ) + 528 )
const PARAM_IO_DIRECTION = ( ( CLASS2 << 16 ) + ( TYPE_ENUM << 24 ) + 529 )
const PARAM_IO_STATE = ( ( CLASS2 << 16 ) + ( TYPE_FLT64 << 24 ) + 530 )
const PARAM_IO_BITDEPTH = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 531 )
const PARAM_GAIN_MULT_FACTOR = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 537 )
const PARAM_GAIN_MULT_ENABLE = ( ( CLASS2 << 16 ) + ( TYPE_BOOLEAN << 24 ) + 541 )
const PARAM_PP_FEAT_NAME = ( ( CLASS2 << 16 ) + ( TYPE_CHAR_PTR << 24 ) + 542 )
const PARAM_PP_INDEX = ( ( CLASS2 << 16 ) + ( TYPE_INT16 << 24 ) + 543 )
const PARAM_ACTUAL_GAIN = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 544 )
const PARAM_PP_PARAM_INDEX = ( ( CLASS2 << 16 ) + ( TYPE_INT16 << 24 ) + 545 )
const PARAM_PP_PARAM_NAME = ( ( CLASS2 << 16 ) + ( TYPE_CHAR_PTR << 24 ) + 546 )
const PARAM_PP_PARAM = ( ( CLASS2 << 16 ) + ( TYPE_UNS32 << 24 ) + 547 )
const PARAM_READ_NOISE = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 548 )
const PARAM_PP_FEAT_ID = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 549 )
const PARAM_PP_PARAM_ID = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 550 )
const PARAM_SMART_STREAM_MODE_ENABLED = ( ( CLASS2 << 16 ) + ( TYPE_BOOLEAN << 24 ) + 700 )
const PARAM_SMART_STREAM_MODE = ( ( CLASS2 << 16 ) + ( TYPE_UNS16 << 24 ) + 701 )
const PARAM_SMART_STREAM_EXP_PARAMS = ( ( CLASS2 << 16 ) + ( TYPE_VOID_PTR << 24 ) + 702 )
const PARAM_SMART_STREAM_DLY_PARAMS = ( ( CLASS2 << 16 ) + ( TYPE_VOID_PTR << 24 ) + 703 )
const PARAM_EXP_TIME = ( ( CLASS3 << 16 ) + ( TYPE_UNS16 << 24 ) + 1 )
const PARAM_EXP_RES = ( ( CLASS3 << 16 ) + ( TYPE_ENUM << 24 ) + 2 )
const PARAM_EXP_RES_INDEX = ( ( CLASS3 << 16 ) + ( TYPE_UNS16 << 24 ) + 4 )
const PARAM_EXPOSURE_TIME = ( ( CLASS3 << 16 ) + ( TYPE_UNS64 << 24 ) + 8 )
const PARAM_BOF_EOF_ENABLE = ( ( CLASS3 << 16 ) + ( TYPE_ENUM << 24 ) + 5 )
const PARAM_BOF_EOF_COUNT = ( ( CLASS3 << 16 ) + ( TYPE_UNS32 << 24 ) + 6 )
const PARAM_BOF_EOF_CLR = ( ( CLASS3 << 16 ) + ( TYPE_BOOLEAN << 24 ) + 7 )
const PARAM_CIRC_BUFFER = ( ( CLASS3 << 16 ) + ( TYPE_BOOLEAN << 24 ) + 299 )
const PARAM_FRAME_BUFFER_SIZE = ( ( CLASS3 << 16 ) + ( TYPE_UNS64 << 24 ) + 300 )
const PARAM_BINNING_SER = ( ( CLASS3 << 16 ) + ( TYPE_ENUM << 24 ) + 165 )
const PARAM_BINNING_PAR = ( ( CLASS3 << 16 ) + ( TYPE_ENUM << 24 ) + 166 )
const PARAM_METADATA_ENABLED = ( ( CLASS3 << 16 ) + ( TYPE_BOOLEAN << 24 ) + 168 )
const PARAM_ROI_COUNT = ( ( CLASS3 << 16 ) + ( TYPE_UNS16 << 24 ) + 169 )
const PARAM_CENTROIDS_ENABLED = ( ( CLASS3 << 16 ) + ( TYPE_BOOLEAN << 24 ) + 170 )
const PARAM_CENTROIDS_RADIUS = ( ( CLASS3 << 16 ) + ( TYPE_UNS16 << 24 ) + 171 )
const PARAM_CENTROIDS_COUNT = ( ( CLASS3 << 16 ) + ( TYPE_UNS16 << 24 ) + 172 )
const PARAM_CENTROIDS_MODE = ( ( CLASS3 << 16 ) + ( TYPE_ENUM << 24 ) + 173 )
const PARAM_CENTROIDS_BG_COUNT = ( ( CLASS3 << 16 ) + ( TYPE_ENUM << 24 ) + 174 )
const PARAM_CENTROIDS_THRESHOLD = ( ( CLASS3 << 16 ) + ( TYPE_UNS32 << 24 ) + 175 )
const PARAM_TRIGTAB_SIGNAL = ( ( CLASS3 << 16 ) + ( TYPE_ENUM << 24 ) + 180 )
const PARAM_LAST_MUXED_SIGNAL = ( ( CLASS3 << 16 ) + ( TYPE_UNS8 << 24 ) + 181 )
const PARAM_FRAME_DELIVERY_MODE = ( ( CLASS3 << 16 ) + ( TYPE_ENUM << 24 ) + 400 )

# temporary placeholder for script hook callback pointer
pm_script_hook = Ptr{Cvoid}

struct PVCAM_FRAME_INFO_GUID
    f1::UInt32
    f2::UInt16
    f3::UInt16
    f4::NTuple{8, UInt8}
end

struct FRAME_INFO
    FrameInfoGUID::PVCAM_FRAME_INFO_GUID
    hCam::Int16
    FrameNr::Int32
    TimeStampe::Int64
    ReadouTime::Int32
    TimeStampBOF::Int64
end

@cenum PL_OPEN_MODES::UInt32 begin
    OPEN_EXCLUSIVE = 0
end

@cenum PL_COOL_MODES::UInt32 begin
    NORMAL_COOL = 0
    CRYO_COOL = 1
    NO_COOL = 2
end

@cenum PL_MPP_MODES::UInt32 begin
    MPP_UNKNOWN = 0
    MPP_ALWAYS_OFF = 1
    MPP_ALWAYS_ON = 2
    MPP_SELECTABLE = 3
end

@cenum PL_SHTR_MODES::UInt32 begin
    SHTR_FAULT = 0
    SHTR_OPENING = 1
    SHTR_OPEN = 2
    SHTR_CLOSING = 3
    SHTR_CLOSED = 4
    SHTR_UNKNOWN = 5
end

@cenum PL_PMODES::UInt32 begin
    PMODE_NORMAL = 0
    PMODE_FT = 1
    PMODE_MPP = 2
    PMODE_FT_MPP = 3
    PMODE_ALT_NORMAL = 4
    PMODE_ALT_FT = 5
    PMODE_ALT_MPP = 6
    PMODE_ALT_FT_MPP = 7
end

@cenum PL_COLOR_MODES::UInt32 begin
    COLOR_NONE = 0
    COLOR_RESERVED = 1
    COLOR_RGGB = 2
    COLOR_GRBG = 3
    COLOR_GBRG = 4
    COLOR_BGGR = 5
end

@cenum PL_IMAGE_FORMATS::UInt32 begin
    PL_IMAGE_FORMAT_MONO16 = 0
    PL_IMAGE_FORMAT_BAYER16 = 1
    PL_IMAGE_FORMAT_MONO8 = 2
    PL_IMAGE_FORMAT_BAYER8 = 3
    PL_IMAGE_FORMAT_MONO24 = 4
    PL_IMAGE_FORMAT_BAYER24 = 5
    PL_IMAGE_FORMAT_RGB24 = 6
    PL_IMAGE_FORMAT_RGB48 = 7
    PL_IMAGE_FORMAT_RGB72 = 8
end

@cenum PL_IMAGE_COMPRESSIONS::UInt32 begin
    PL_IMAGE_COMPRESSION_NONE = 0
    PL_IMAGE_COMPRESSION_RESERVED8 = 8
    PL_IMAGE_COMPRESSION_BITPACK9 = 9
    PL_IMAGE_COMPRESSION_BITPACK10 = 10
    PL_IMAGE_COMPRESSION_BITPACK11 = 11
    PL_IMAGE_COMPRESSION_BITPACK12 = 12
    PL_IMAGE_COMPRESSION_BITPACK13 = 13
    PL_IMAGE_COMPRESSION_BITPACK14 = 14
    PL_IMAGE_COMPRESSION_BITPACK15 = 15
    PL_IMAGE_COMPRESSION_RESERVED16 = 16
    PL_IMAGE_COMPRESSION_BITPACK17 = 17
    PL_IMAGE_COMPRESSION_BITPACK18 = 18
    PL_IMAGE_COMPRESSION_RESERVED24 = 24
    PL_IMAGE_COMPRESSION_RESERVED32 = 32
end

@cenum PL_PARAM_ATTRIBUTES::UInt32 begin
    ATTR_CURRENT = 0
    ATTR_COUNT = 1
    ATTR_TYPE = 2
    ATTR_MIN = 3
    ATTR_MAX = 4
    ATTR_DEFAULT = 5
    ATTR_INCREMENT = 6
    ATTR_ACCESS = 7
    ATTR_AVAIL = 8
end

@cenum PL_PARAM_ACCESS::UInt32 begin
    ACC_READ_ONLY = 1
    ACC_READ_WRITE = 2
    ACC_EXIST_CHECK_ONLY = 3
    ACC_WRITE_ONLY = 4
end

@cenum PL_IO_TYPE::UInt32 begin
    IO_TYPE_TTL = 0
    IO_TYPE_DAC = 1
end

@cenum PL_IO_DIRECTION::UInt32 begin
    IO_DIR_INPUT = 0
    IO_DIR_OUTPUT = 1
    IO_DIR_INPUT_OUTPUT = 2
end

@cenum PL_READOUT_PORTS::UInt32 begin
    READOUT_PORT_0 = 0
    READOUT_PORT_1 = 1
end

@cenum PL_CLEAR_MODES::UInt32 begin
    CLEAR_NEVER = 0
    CLEAR_PRE_EXPOSURE = 1
    CLEAR_PRE_SEQUENCE = 2
    CLEAR_POST_SEQUENCE = 3
    CLEAR_PRE_POST_SEQUENCE = 4
    CLEAR_PRE_EXPOSURE_POST_SEQ = 5
    MAX_CLEAR_MODE = 6
end

@cenum PL_SHTR_OPEN_MODES::UInt32 begin
    OPEN_NEVER = 0
    OPEN_PRE_EXPOSURE = 1
    OPEN_PRE_SEQUENCE = 2
    OPEN_PRE_TRIGGER = 3
    OPEN_NO_CHANGE = 4
end

@cenum PL_EXPOSURE_MODES::UInt32 begin
    TIMED_MODE = 0
    STROBED_MODE = 1
    BULB_MODE = 2
    TRIGGER_FIRST_MODE = 3
    FLASH_MODE = 4
    VARIABLE_TIMED_MODE = 5
    INT_STROBE_MODE = 6
    MAX_EXPOSE_MODE = 7
    EXT_TRIG_INTERNAL = 1792
    EXT_TRIG_TRIG_FIRST = 2048
    EXT_TRIG_EDGE_RISING = 2304
end

@cenum PL_EXPOSE_OUT_MODES::UInt32 begin
    EXPOSE_OUT_FIRST_ROW = 0
    EXPOSE_OUT_ALL_ROWS = 1
    EXPOSE_OUT_ANY_ROW = 2
    EXPOSE_OUT_ROLLING_SHUTTER = 3
    MAX_EXPOSE_OUT_MODE = 4
end

@cenum PL_FAN_SPEEDS::UInt32 begin
    FAN_SPEED_HIGH = 0
    FAN_SPEED_MEDIUM = 1
    FAN_SPEED_LOW = 2
    FAN_SPEED_OFF = 3
end

@cenum PL_TRIGTAB_SIGNALS::UInt32 begin
    PL_TRIGTAB_SIGNAL_EXPOSE_OUT = 0
end

@cenum PL_FRAME_DELIVERY_MODES::UInt32 begin
    PL_FRAME_DELIVERY_MODE_MAX_FPS = 0
    PL_FRAME_DELIVERY_MODE_CONSTANT_INTERVALS = 1
end

@cenum PL_CAM_INTERFACE_TYPES::UInt32 begin
    PL_CAM_IFC_TYPE_UNKNOWN = 0
    PL_CAM_IFC_TYPE_1394 = 256
    PL_CAM_IFC_TYPE_1394_A = 257
    PL_CAM_IFC_TYPE_1394_B = 258
    PL_CAM_IFC_TYPE_USB = 512
    PL_CAM_IFC_TYPE_USB_1_1 = 513
    PL_CAM_IFC_TYPE_USB_2_0 = 514
    PL_CAM_IFC_TYPE_USB_3_0 = 515
    PL_CAM_IFC_TYPE_USB_3_1 = 516
    PL_CAM_IFC_TYPE_PCI = 1024
    PL_CAM_IFC_TYPE_PCI_LVDS = 1025
    PL_CAM_IFC_TYPE_PCIE = 2048
    PL_CAM_IFC_TYPE_PCIE_LVDS = 2049
    PL_CAM_IFC_TYPE_PCIE_X1 = 2050
    PL_CAM_IFC_TYPE_PCIE_X4 = 2051
    PL_CAM_IFC_TYPE_VIRTUAL = 4096
    PL_CAM_IFC_TYPE_ETHERNET = 8192
end

@cenum PL_CAM_INTERFACE_MODES::UInt32 begin
    PL_CAM_IFC_MODE_UNSUPPORTED = 0
    PL_CAM_IFC_MODE_CONTROL_ONLY = 1
    PL_CAM_IFC_MODE_IMAGING = 2
end

@cenum PL_CENTROIDS_MODES::UInt32 begin
    PL_CENTROIDS_MODE_LOCATE = 0
    PL_CENTROIDS_MODE_TRACK = 1
    PL_CENTROIDS_MODE_BLOB = 2
end

@cenum PL_SCAN_MODES::UInt32 begin
    PL_SCAN_MODE_AUTO = 0
    PL_SCAN_MODE_PROGRAMMABLE_LINE_DELAY = 1
    PL_SCAN_MODE_PROGRAMMABLE_SCAN_WIDTH = 2
end

@cenum PL_SCAN_DIRECTIONS::UInt32 begin
    PL_SCAN_DIRECTION_DOWN = 0
    PL_SCAN_DIRECTION_UP = 1
    PL_SCAN_DIRECTION_DOWN_UP = 2
end

@cenum PP_FEATURE_IDS::UInt32 begin
    PP_FEATURE_RING_FUNCTION = 0
    PP_FEATURE_BIAS = 1
    PP_FEATURE_BERT = 2
    PP_FEATURE_QUANT_VIEW = 3
    PP_FEATURE_BLACK_LOCK = 4
    PP_FEATURE_TOP_LOCK = 5
    PP_FEATURE_VARI_BIT = 6
    PP_FEATURE_RESERVED = 7
    PP_FEATURE_DESPECKLE_BRIGHT_HIGH = 8
    PP_FEATURE_DESPECKLE_DARK_LOW = 9
    PP_FEATURE_DEFECTIVE_PIXEL_CORRECTION = 10
    PP_FEATURE_DYNAMIC_DARK_FRAME_CORRECTION = 11
    PP_FEATURE_HIGH_DYNAMIC_RANGE = 12
    PP_FEATURE_DESPECKLE_BRIGHT_LOW = 13
    PP_FEATURE_DENOISING = 14
    PP_FEATURE_DESPECKLE_DARK_HIGH = 15
    PP_FEATURE_ENHANCED_DYNAMIC_RANGE = 16
    PP_FEATURE_MAX = 17
end

@cenum PP_PARAMETER_IDS::UInt32 begin
    PP_PARAMETER_RF_FUNCTION = 0
    PP_FEATURE_BIAS_ENABLED = 10
    PP_FEATURE_BIAS_LEVEL = 11
    PP_FEATURE_BERT_ENABLED = 20
    PP_FEATURE_BERT_THRESHOLD = 21
    PP_FEATURE_QUANT_VIEW_ENABLED = 30
    PP_FEATURE_QUANT_VIEW_E = 31
    PP_FEATURE_BLACK_LOCK_ENABLED = 40
    PP_FEATURE_BLACK_LOCK_BLACK_CLIP = 41
    PP_FEATURE_TOP_LOCK_ENABLED = 50
    PP_FEATURE_TOP_LOCK_WHITE_CLIP = 51
    PP_FEATURE_VARI_BIT_ENABLED = 60
    PP_FEATURE_VARI_BIT_BIT_DEPTH = 61
    PP_FEATURE_DESPECKLE_BRIGHT_HIGH_ENABLED = 80
    PP_FEATURE_DESPECKLE_BRIGHT_HIGH_THRESHOLD = 81
    PP_FEATURE_DESPECKLE_BRIGHT_HIGH_MIN_ADU_AFFECTED = 82
    PP_FEATURE_DESPECKLE_DARK_LOW_ENABLED = 90
    PP_FEATURE_DESPECKLE_DARK_LOW_THRESHOLD = 91
    PP_FEATURE_DESPECKLE_DARK_LOW_MAX_ADU_AFFECTED = 92
    PP_FEATURE_DEFECTIVE_PIXEL_CORRECTION_ENABLED = 100
    PP_FEATURE_DYNAMIC_DARK_FRAME_CORRECTION_ENABLED = 110
    PP_FEATURE_HIGH_DYNAMIC_RANGE_ENABLED = 120
    PP_FEATURE_DESPECKLE_BRIGHT_LOW_ENABLED = 130
    PP_FEATURE_DESPECKLE_BRIGHT_LOW_THRESHOLD = 131
    PP_FEATURE_DESPECKLE_BRIGHT_LOW_MAX_ADU_AFFECTED = 132
    PP_FEATURE_DENOISING_ENABLED = 140
    PP_FEATURE_DENOISING_NO_OF_ITERATIONS = 141
    PP_FEATURE_DENOISING_GAIN = 142
    PP_FEATURE_DENOISING_OFFSET = 143
    PP_FEATURE_DENOISING_LAMBDA = 144
    PP_FEATURE_DESPECKLE_DARK_HIGH_ENABLED = 150
    PP_FEATURE_DESPECKLE_DARK_HIGH_THRESHOLD = 151
    PP_FEATURE_DESPECKLE_DARK_HIGH_MIN_ADU_AFFECTED = 152
    PP_FEATURE_ENHANCED_DYNAMIC_RANGE_ENABLED = 160
    PP_PARAMETER_ID_MAX = 161
end

struct smart_stream_type
    entries::UInt16
    params::Ptr{UInt32}
end

@cenum PL_SMT_MODES::UInt32 begin
    SMTMODE_ARBITRARY_ALL = 0
    SMTMODE_MAX = 1
end

@cenum PL_IMAGE_STATUSES::UInt32 begin
    READOUT_NOT_ACTIVE = 0
    EXPOSURE_IN_PROGRESS = 1
    READOUT_IN_PROGRESS = 2
    READOUT_COMPLETE = 3
    FRAME_AVAILABLE = 3
    READOUT_FAILED = 4
    ACQUISITION_IN_PROGRESS = 5
    MAX_CAMERA_STATUS = 6
end

@cenum PL_CCS_ABORT_MODES::UInt32 begin
    CCS_NO_CHANGE = 0
    CCS_HALT = 1
    CCS_HALT_CLOSE_SHTR = 2
    CCS_CLEAR = 3
    CCS_CLEAR_CLOSE_SHTR = 4
    CCS_OPEN_SHTR = 5
    CCS_CLEAR_OPEN_SHTR = 6
end

@cenum PL_IRQ_MODES::UInt32 begin
    NO_FRAME_IRQS = 0
    BEGIN_FRAME_IRQS = 1
    END_FRAME_IRQS = 2
    BEGIN_END_FRAME_IRQS = 3
end

@cenum PL_CIRC_MODES::UInt32 begin
    CIRC_NONE = 0
    CIRC_OVERWRITE = 1
    CIRC_NO_OVERWRITE = 2
end

@cenum PL_EXP_RES_MODES::UInt32 begin
    EXP_RES_ONE_MILLISEC = 0
    EXP_RES_ONE_MICROSEC = 1
    EXP_RES_ONE_SEC = 2
end

@cenum PL_SRC_MODES::UInt32 begin
    SCR_PRE_OPEN_SHTR = 0
    SCR_POST_OPEN_SHTR = 1
    SCR_PRE_FLASH = 2
    SCR_POST_FLASH = 3
    SCR_PRE_INTEGRATE = 4
    SCR_POST_INTEGRATE = 5
    SCR_PRE_READOUT = 6
    SCR_POST_READOUT = 7
    SCR_PRE_CLOSE_SHTR = 8
    SCR_POST_CLOSE_SHTR = 9
end

@cenum PL_CALLBACK_EVENT::UInt32 begin
    PL_CALLBACK_BOF = 0
    PL_CALLBACK_EOF = 1
    PL_CALLBACK_CHECK_CAMS = 2
    PL_CALLBACK_CAM_REMOVED = 3
    PL_CALLBACK_CAM_RESUMED = 4
    PL_CALLBACK_MAX = 5
end

struct rgn_type
    s1::UInt16
    s2::UInt16
    sbin::UInt16
    p1::UInt16
    p2::UInt16
    pbin::UInt16
end

struct io_struct
    io_port::UInt16
    io_type::UInt32
    state::Cdouble
    next::Ptr{io_struct}
end

const io_entry = io_struct

struct io_list
    pre_open::io_entry
    post_open::io_entry
    pre_flash::io_entry
    post_flash::io_entry
    pre_integrate::io_entry
    post_integrate::io_entry
    pre_readout::io_entry
    post_readout::io_entry
    pre_close::io_entry
    post_close::io_entry
end

struct active_camera_type
    shutter_close_delay::UInt16
    shutter_open_delay::UInt16
    rows::UInt16
    cols::UInt16
    prescan::UInt16
    postscan::UInt16
    premask::UInt16
    postmask::UInt16
    preflash::UInt16
    clear_count::UInt16
    preamp_delay::UInt16
    mpp_selectable::UInt16
    frame_selectable::UInt16
    do_clear::Int16
    open_shutter::Int16
    mpp_mode::UInt16
    frame_transfer::UInt16
    alt_mode::UInt16
    exp_res::UInt32
    io_hdr::Ptr{io_list}
end

@cenum PL_MD_FRAME_FLAGS::UInt32 begin
    PL_MD_FRAME_FLAG_ROI_TS_SUPPORTED = 1
    PL_MD_FRAME_FLAG_UNUSED_2 = 2
    PL_MD_FRAME_FLAG_UNUSED_3 = 4
    PL_MD_FRAME_FLAG_UNUSED_4 = 16
    PL_MD_FRAME_FLAG_UNUSED_5 = 32
    PL_MD_FRAME_FLAG_UNUSED_6 = 64
    PL_MD_FRAME_FLAG_UNUSED_7 = 128
end

@cenum PL_MD_ROI_FLAGS::UInt32 begin
    PL_MD_ROI_FLAG_INVALID = 1
    PL_MD_ROI_FLAG_HEADER_ONLY = 2
    PL_MD_ROI_FLAG_UNUSED_3 = 4
    PL_MD_ROI_FLAG_UNUSED_4 = 16
    PL_MD_ROI_FLAG_UNUSED_5 = 32
    PL_MD_ROI_FLAG_UNUSED_6 = 64
    PL_MD_ROI_FLAG_UNUSED_7 = 128
end

struct md_frame_header
    signature::UInt32
    version::UInt8
    frameNr::UInt32
    roiCount::UInt16
    timestampBOF::UInt32
    timestampEOF::UInt32
    timestampResNs::UInt32
    exposureTime::UInt32
    exposureTimeResNs::UInt32
    roiTimestampResNs::UInt32
    bitDepth::UInt8
    colorMask::UInt8
    flags::UInt8
    extendedMdSize::UInt16
    imageFormat::UInt8
    imageCompression::UInt8
    _reserved::NTuple{6, UInt8}
end

struct md_frame_roi_header
    roiNr::UInt16
    timestampBOR::UInt32
    timestampEOR::UInt32
    roi::rgn_type
    flags::UInt8
    extendedMdSize::UInt16
    roiDataSize::UInt32
    _reserved::NTuple{3, UInt8}
end

@cenum PL_MD_EXT_TAGS::UInt32 begin
    PL_MD_EXT_TAG_PARTICLE_ID = 0
    PL_MD_EXT_TAG_PARTICLE_M0 = 1
    PL_MD_EXT_TAG_PARTICLE_M2 = 2
    PL_MD_EXT_TAG_MAX = 3
end

struct md_ext_item_info
    tag::PL_MD_EXT_TAGS
    type::UInt16
    size::UInt16
    name::Cstring
end

struct md_ext_item
    tagInfo::Ptr{md_ext_item_info}
    value::Ptr{Cvoid}
end

struct md_ext_item_collection
    list::NTuple{255, md_ext_item}
    map::NTuple{255, Ptr{md_ext_item}}
    count::UInt16
end

struct md_frame_roi
    header::Ptr{md_frame_roi_header}
    data::Ptr{Cvoid}
    dataSize::UInt32
    extMdData::Ptr{Cvoid}
    extMdDataSize::UInt16
end

struct md_frame
    header::Ptr{md_frame_header}
    extMdData::Ptr{Cvoid}
    extMdDataSize::UInt16
    impliedRoi::rgn_type
    roiArray::Ptr{md_frame_roi}
    roiCapacity::UInt16
    roiCount::UInt16
end

const PPVCAM_FRAME_INFO_GUID = Ptr{PVCAM_FRAME_INFO_GUID}
const PFRAME_INFO = Ptr{FRAME_INFO}

const PL_CALLBACK_SIG_LEGACY = Ptr{Cvoid}
const PL_CALLBACK_SIG_EX = Ptr{Cvoid}
const PL_CALLBACK_SIG_EX2 = Ptr{Cvoid}
const PL_CALLBACK_SIG_EX3 = Ptr{Cvoid}


text:
module LinAlgExtensions

using LinearAlgebra
using SparseArrays
using SuiteSparse
using Random

include("utils.jl")
include("randutils.jl")
include("pinv.jl")

end # module


text:
module Microscopy

include("microscope_modes.jl")
include("objectives.jl")

end

max_stars_repo_name: JuliaMicroscopy/Microscopy

max_forks_repo_name: JuliaMicroscopy/Microscopy


text:
module TestMGT

include("square.jl")

using Test
import TreeLas.MGT
import GraphIdx
import GraphIdx: EdgeGraph
import GraphIdx.Io: prettyarr

Base.sort(i::T, j::T) where {T} =
    i < j ? (i, j) : (j, i)

@testset "MGT     : Square               " begin
    edges, n = square_edges()
    m = length(edges)
    lambda = GraphIdx.Const(0.5)
    y = Float64[0 1; 0 0.5]
    x = MGT.gaplas(y, edges, lambda)
    @test size(x) == size(y)
    @test sum(x) ≈ sum(y)
    @test x ≈ mean(y) * ones(size(y)...)
end

@testset "MGT     : demo 3x7             " begin
    include("demo3x7.jl")
    n1, n2 = 3, 7
    n = n1 * n2
    m = 32
    root = 1
    s = collect(1:n) |> x -> reshape(x, n1, n2)
    @test prettyarr(s) == """
    3×7 Array{Int64,2}:
     1  4  7  10  13  16  19
     2  5  8  11  14  17  20
     3  6  9  12  15  18  21
    """ |> chomp
    @test size(y) == (n1,n2)

    edges = GraphIdx.Edge{Int}[
        (1,  2),  (4,  5),  (7, 8), (10, 11), (13, 14), (16, 17),
        (19, 20), (2,  3),  (5, 6), (8,   9), (11, 12), (14, 15),
        (17, 18), (20, 21), (1, 4), (4,   7), (7,  10), (10, 13),
        (13, 16), (16, 19), (2, 5), (5,   8), (8,  11), (11, 14),
        (14, 17), (17, 20), (3, 6), (6,   9), (9,  12), (12, 15),
        (15, 18), (18, 21)
    ]

    @test 7*2 + 6*3 == m
    @test length(edges) == m

    γ = [-0.0,  -1.0,  -0.24, -0.04, -0.43, -0.73,
         -0.7,  -0.28, -1.41, -1.43, -0.51, -0.52,
         -1.14, -0.02, -0.11, -0.02, -0.79, -0.33,
         -0.74, -0.65, -1.11, -0.78, -0.51, -0.14,
         -0.44, -0.78, -0.58, -0.80, -1.43, -0.13,
         -1.06, -0.38]

    @testset "compute γ manually" begin
        gam = [-abs(y[u] - y[v]) for (u, v) in edges]
        @test γ ≈ gam
    end

    tree_mask = GraphIdx.kruskal_mst(n, edges, γ)
    tree_cost = γ[tree_mask]
    tree_edges = edges[tree_mask]
    idx = GraphIdx.IncidenceIndex(n, edges)
    local parent::Vector{Int}, selected =
        GraphIdx.prim_mst_edges(n, edges, γ, root)

    @test GraphIdx.Tree.find_root(parent) == root
    @test parent[root] == root
    @test parent ==
        [1, 5, 6, 1, 4, 5, 8, 9, 6, 7, 8, 9, 10, 15, 18, 13, 16, 17, 20, 17, 18]
    @test -sum(tree_cost) ≈ 15.79

    prim_edges = [sort(i, j) for (i, j) in enumerate(parent) if i != j]
    @test length(prim_edges) == n-1
    s_prim, s_kruskal = Set(prim_edges), Set(tree_edges)
    prim_mask = [e ∈ s_prim for e in edges]
    prim_cost = γ[prim_mask]
    @test length(s_prim) == n-1
    @test length(s_kruskal) == n-1
    @test -sum(prim_cost) ≈ 15.79

    sdiff = symdiff(s_prim, s_kruskal)
    sdiff_mask = [e ∈ sdiff for e in edges]
    sdiff_cost = Set(γ[sdiff_mask])
    @test length(sdiff_cost) == 1
    @test sdiff_cost == Set([-0.51])

    @test length(selected) == n - 1
    @test Set(edges[selected]) == s_prim

    @testset "extract non-tree" begin
        z = copy(y)
        lambda = GraphIdx.Ones{Float64}()
        tlam = fill(NaN, length(y))
        alpha = zeros(length(edges))
        MGT.extract_non_tree!(EdgeGraph(n, edges), parent, z, alpha, tlam, lambda)
        @test z ≈ y
        @test all(tlam[1:end .!= root] .≈ 1.0)
    end

    @testset "call gaplas" begin
        local lam = GraphIdx.Ones{Float64}()
        x = MGT.gaplas(y, edges, lam, max_iter=1, verbose=false)
        @test rationalize.(x) == [
            639//800  639//800  547//650  547//650  547//650  547//650  547//650
            639//800  639//800  547//650  547//650  547//650  547//650  547//650
            639//800  639//800  639//800  639//800  547//650  547//650  547//650
        ]
        x = MGT.gaplas(y, edges, lam, max_iter=5, verbose=true)
        @test mean(x) ≈ mean(y)
    end
end

end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
function fast_non_dominated_sort!(pop)
    n = length(pop)

    dom_list = [ Int[] for i in 1:n ]
    rank = zeros(Int, n)
    dom_count = zeros(Int, n)

    for i in 1:n
        for j in i+1:n
            comparison = compare(pop[i], pop[j])
            if comparison == 1 #is_better(pop[i], pop[j])
                push!(dom_list[i], j)
                dom_count[j] += 1
            elseif comparison == 2 #is_better(pop[j], pop[i])
                push!(dom_list[j], i)
                dom_count[i] += 1
            end
        end
        if dom_count[i] == 0
            rank[i] = 1
        end
    end

    k = UInt16(2)
    while any(==(k-one(UInt16)), (rank[p] for p in 1:n)) #ugly workaround for #15276
        for p in 1:n
            if rank[p] == k-one(UInt16)
                for q in dom_list[p]
                    dom_count[q] -= one(UInt16)
                    if dom_count[q] == zero(UInt16)
                        rank[q] = k
                    end
                end
            end
        end
        k += one(UInt16)
    end

    for i in 1:length(pop)
        pop[i].rank = rank[i]
    end

    sort!(pop, by = s -> s.rank, alg = Base.Sort.QuickSort)
    pop
end

function update_crowding_distance!(pop)

    for sol in pop
        sol.crowding = 0.0
    end

    if length(pop[1].f) == 2
        sort!(pop, by = x -> x.f[1])
        pop[1].f[1] == pop[end].f[1] && return #Don't waste time if all indivs are the same
        pop[1].crowding = pop[end].crowding = Inf

        width_y1 = (pop[end].f[1] - pop[1].f[1])
        width_y2 = (pop[1].f[2] - pop[end].f[2])
        for i = 2:length(pop)-1
            pop[i].crowding = (pop[i+1].f[1] - pop[i-1].f[1]) / width_y1 + (pop[i-1].f[2] - pop[i+1].f[2]) / width_y2
        end

        return pop
    end

    for j = 1:length(first(pop).f) # Foreach objective
        let j = j
            sort!(pop, by = x -> x.f[j]) #sort by the objective value
        end
        pop[1].crowding = pop[end].crowding = Inf #Assign infinite value to extremas
        if pop[1].f[j] != pop[end].f[j]
            for i = 2:length(pop)-1
                pop[i].crowding += (pop[i+1].f[j] - pop[i-1].f[j]) / (pop[end].f[j] - pop[1].f[j])
            end
        end
    end
    pop
end

function truncate_population!(population, N)
    fast_non_dominated_sort!(population)

    let f::Int = 1
        ind = 0
        indnext = findlast(x -> x.rank == f, population)
        while 0 < indnext <= N
            ind = indnext
            f += 1
            indnext = findlast(x -> x.rank == f, population)
        end
        indnext == 0 && (indnext = length(population))
        update_crowding_distance!(view(population, ind+1:indnext))
        sort!(view(population, ind+1:indnext), by = x -> x.crowding, rev = true, alg = PartialQuickSort(N-ind))
    end
    deleteat!(population, N+1:length(population))
end


text:
#####
##### Linear
#####

"""
    fit_echo_linear(
        phas::AbstractArray{<:AbstractFloat, N > 1},
        W::AbstractArray{<:AbstractFloat, N > 1},
        TEs::NTuple{NT > 1, Real};
        phase_offset::Bool = true
    ) -> Tuple{typeof(similar(phas)){N-1}, [typeof(similar(phas)){N-1}]}

Weighted least squares for multi-echo data.

### Arguments
- `phas::AbstractArray{<:AbstractFloat, N > 1}`: unwrapped multi-echo phase
- `W::AbstractArray{<:AbstractFloat, N > 1}`: reciprocal of error variance of voxel
- `TEs::NTuple{NT > 1, Real}`: echo times

### Keywords
- `phase_offset::Bool = true`: model phase offset (`true`)

### Returns
- `typeof(similar(phas)){N-1}`: weighted least-squares estimate for phase
- [`typeof(similar(phas)){N-1}`]: weighted least-squares estimate for phase offset
    if `phase_offset = true`
"""
function fit_echo_linear(
    phas::AbstractArray{<:AbstractFloat, N},
    W::AbstractArray{<:AbstractFloat, N},
    TEs::NTuple{NT, Real};
    phase_offset::Bool = true
) where {N, NT}
    N > 1 || throw(ArgumentError("array must contain echoes in last dimension"))
    NT > 1 || throw(ArgumentError("data must be multi-echo"))

    p = similar(phas, size(phas)[1:N-1])
    if !phase_offset
        return fit_echo_linear!(p, phas, W, TEs)
    else
        return fit_echo_linear!(p, similar(p), phas, W, TEs)
    end
end

"""
    fit_echo_linear!(
        p::AbstractArray{<:AbstractFloat, N},
        phas::AbstractArray{<:AbstractFloat, M > 1},
        W::AbstractArray{<:AbstractFloat, M > 1},
        TEs::NTuple{NT > 1, Real}
    ) -> p

Weighted least squares for multi-echo data (phase offset = 0).

### Arguments
- `p::AbstractArray{<:AbstractFloat, N}`: weighted least-squares estimate for phase
- `phas::AbstractArray{<:AbstractFloat, M > 1}`: unwrapped multi-echo phase
- `W::AbstractArray{<:AbstractFloat, M > 1}`: reciprocal of error variance of voxel
- `TEs::NTuple{NT > 1, Real}`: echo times

### Returns
- `p`: weighted least-squares estimate for phase
"""
function fit_echo_linear!(
    p::AbstractArray{Tp, N},
    phas::AbstractArray{Tphas, M},
    W::AbstractArray{TW, M},
    TEs::NTuple{NT, Real}
) where {Tp<:AbstractFloat, Tphas<:AbstractFloat, TW<:Real, N, M, NT}
    M > 1 || throw(ArgumentError("array must contain echoes in last dimension"))
    NT > 1 || throw(ArgumentError("data must be multi-echo"))

    size(phas, M) == NT || throw(DimensionMismatch())
    size(W) == size(phas) || throw(DimensionMismatch())
    length(p) == length(phas) ÷ NT || throw(DimensionMismatch())

    vphas = reshape(phas, :, NT)
    vW = reshape(W, :, NT)
    vp = vec(p)

    T = promote_type(Tp, Tphas, TW)
    tes = convert.(T, TEs)

    _zeroTp = zero(Tp)

    @inbounds @batch for I in eachindex(vp)
        w = vW[I,1] * vW[I,1] * tes[1]
        den = w * tes[1]
        num = w * vphas[I,1]
        for t in 2:NT
            w = vW[I,t] * vW[I,t] * tes[t]
            den = muladd(w, tes[t], den)
            num = muladd(w, vphas[I,t], num)
        end
        vp[I] = iszero(den) ? _zeroTp : num * inv(den)
    end

    return p
end

"""
    fit_echo_linear!(
        p::AbstractArray{<:AbstractFloat, N},
        p0::AbstractArray{<:AbstractFloat, N},
        phas::AbstractArray{<:AbstractFloat, M > 1},
        W::AbstractArray{<:AbstractFloat, M > 1},
        TEs::NTuple{NT > 1, Real}
    ) -> (p, p0)

Weighted least squares for multi-echo data (estimate phase offset).

### Arguments
- `p::AbstractArray{<:AbstractFloat, N}`: weighted least-squares estimate for phase
- `p0::AbstractArray{<:AbstractFloat, N}`: weighted least-squares estimate for phase offset
- `phas::AbstractArray{<:AbstractFloat, M > 1}`: unwrapped multi-echo phase
- `W::AbstractArray{<:AbstractFloat, M > 1}`: reciprocal of error variance of voxel
- `TEs::NTuple{NT > 1, Real}`: echo times

### Returns
- `p`: weighted least-squares estimate for phase
- `p0`: weighted least-squares estimate for phase offset
"""
function fit_echo_linear!(
    p::AbstractArray{Tp, N},
    p0::AbstractArray{Tp0, N},
    phas::AbstractArray{Tphas, M},
    W::AbstractArray{TW, M},
    TEs::NTuple{NT, Real}
) where {Tp<:AbstractFloat, Tp0<:AbstractFloat, Tphas<:AbstractFloat, TW<:Real, N, M, NT}
    M > 1 || throw(ArgumentError("array must contain echoes in last dimension"))
    NT > 1 || throw(ArgumentError("data must be multi-echo"))

    size(phas, M) == NT || throw(DimensionMismatch())
    size(W) == size(phas) || throw(DimensionMismatch())
    length(p) == length(phas) ÷ NT || throw(DimensionMismatch())
    length(p0) == length(phas) ÷ NT || throw(DimensionMismatch())

    vphas = reshape(phas, :, NT)
    vW = reshape(W, :, NT)
    vp0 = vec(p0)
    vp = vec(p)

    T = promote_type(Tp, Tp0, Tphas, TW)
    tes = convert.(T, TEs)

    _zeroTp = zero(Tp)
    _zeroTp0 = zero(Tp0)

    @inbounds @batch for I in eachindex(vp)
        w = vW[I,1]
        x = w * tes[1]
        y = w * vphas[I,1]
        for t in 2:NT
            w += vW[I,t]
            x = muladd(vW[I,t], tes[t], x)
            y = muladd(vW[I,t], vphas[I,t], y)
        end

        if iszero(w)
            vp[I] = _zeroTp
            vp0[I] = _zeroTp0
            continue
        end

        w = inv(w)
        x *= w
        y *= w

        xx = tes[1] - x
        yy = vphas[I,1] - y
        ww = vW[I,1] * xx
        num = ww * yy
        den = ww * xx
        for t in 2:NT
            xx = tes[t] - x
            yy = vphas[I,t] - y
            ww = vW[I,t] * xx
            num = muladd(ww, yy, num)
            den = muladd(ww, xx, den)
        end

        if !iszero(den)
            vp[I] = num * inv(den)
            vp0[I] = y - vp[I] * x
        else
            vp[I] = _zeroTp
            vp0[I] = _zeroTp0
        end
    end

    return p, p0
end


text:
using Test
import ClimateMachine
using ClimateMachine.VariableTemplates: @vars, Vars
using ClimateMachine.Mesh.Grids:
    EveryDirection, HorizontalDirection, VerticalDirection
using ClimateMachine.MPIStateArrays: weightedsum
using ClimateMachine.BalanceLaws

import GaussQuadrature
using MPI
using LinearAlgebra

ClimateMachine.init()

@testset "Exponential and Cutoff filter matrix" begin
    let
        # Values computed with:
        #   https://github.com/tcew/nodal-dg/blob/master/Codes1.1/Codes1D/Filter1D.m
        #! format: off
        W = [0x3fe98f3cd0d725e8  0x3fddfd863c6c9a44  0xbfe111110d0fd334  0x3fddbe357bce0b5c  0xbfc970267f929618
             0x3fb608a150f6f927  0x3fe99528b1a1cd8d  0x3fcd41d41f8bae45  0xbfc987d5fabab8d5  0x3fb5da1cd858af87
             0xbfb333332eb1cd92  0x3fc666666826f178  0x3fe999999798faaa  0x3fc666666826f176  0xbfb333332eb1cd94
             0x3fb5da1cd858af84  0xbfc987d5fabab8d4  0x3fcd41d41f8bae46  0x3fe99528b1a1cd8e  0x3fb608a150f6f924
             0xbfc970267f929618  0x3fddbe357bce0b5c  0xbfe111110d0fd333  0x3fddfd863c6c9a44  0x3fe98f3cd0d725e8]
        #! format: on
        W = reinterpret.(Float64, W)

        N = size(W, 1) - 1

        topology = ClimateMachine.Mesh.Topologies.BrickTopology(
            MPI.COMM_SELF,
            -1.0:2.0:1.0,
        )

        grid = ClimateMachine.Mesh.Grids.DiscontinuousSpectralElementGrid(
            topology;
            polynomialorder = N,
            FloatType = Float64,
            DeviceArray = Array,
        )

        filter = ClimateMachine.Mesh.Filters.ExponentialFilter(grid, 0, 32)
        nf = length(filter.filter_matrices)
        @test all(ntuple(i -> filter.filter_matrices[i] ≈ W, nf))
    end

    let
        # Values computed with:
        #   https://github.com/tcew/nodal-dg/blob/master/Codes1.1/Codes1D/Filter1D.m
        #! format: off
        W = [0x3fd822e5f54ecb62   0x3fedd204a0f08ef8   0xbfc7d3aa58fd6968   0xbfbf74682ac4d276
             0x3fc7db36e726d8c1   0x3fe59d16feee478b   0x3fc6745bfbb91e20   0xbfa30fbb7a645448
             0xbfa30fbb7a645455   0x3fc6745bfbb91e26   0x3fe59d16feee478a   0x3fc7db36e726d8c4
             0xbfbf74682ac4d280   0xbfc7d3aa58fd6962   0x3fedd204a0f08ef7   0x3fd822e5f54ecb62]
        #! format: on
        W = reinterpret.(Float64, W)

        N = size(W, 1) - 1

        topology = ClimateMachine.Mesh.Topologies.BrickTopology(
            MPI.COMM_SELF,
            -1.0:2.0:1.0,
        )
        grid = ClimateMachine.Mesh.Grids.DiscontinuousSpectralElementGrid(
            topology;
            polynomialorder = N,
            FloatType = Float64,
            DeviceArray = Array,
        )

        filter = ClimateMachine.Mesh.Filters.ExponentialFilter(grid, 1, 4)
        nf = length(filter.filter_matrices)
        @test all(ntuple(i -> filter.filter_matrices[i] ≈ W, nf))
    end

    let
        T = Float64
        N = (5, 3)
        Nc = (4, 2)

        topology = ClimateMachine.Mesh.Topologies.BrickTopology(
            MPI.COMM_SELF,
            -1.0:2.0:1.0,
        )
        grid = ClimateMachine.Mesh.Grids.DiscontinuousSpectralElementGrid(
            topology;
            polynomialorder = N,
            FloatType = T,
            DeviceArray = Array,
        )

        ξ = ClimateMachine.Mesh.Grids.referencepoints(grid)
        ξ1 = ξ[1]
        ξ2 = ξ[2]
        a1, b1 = GaussQuadrature.legendre_coefs(T, N[1])
        a2, b2 = GaussQuadrature.legendre_coefs(T, N[2])
        V1 = GaussQuadrature.orthonormal_poly(ξ1, a1, b1)
        V2 = GaussQuadrature.orthonormal_poly(ξ2, a2, b2)

        Σ1 = ones(T, N[1] + 1)
        Σ2 = ones(T, N[2] + 1)
        Σ1[(Nc[1]:N[1]) .+ 1] .= 0
        Σ2[(Nc[2]:N[2]) .+ 1] .= 0

        W1 = V1 * Diagonal(Σ1) / V1
        W2 = V2 * Diagonal(Σ2) / V2

        filter = ClimateMachine.Mesh.Filters.CutoffFilter(grid, Nc)
        @test filter.filter_matrices[1] ≈ W1
        @test filter.filter_matrices[2] ≈ W2
    end

    let
        T = Float64
        N = (5, 3)
        Nc = (4, 2)

        topology = ClimateMachine.Mesh.Topologies.BrickTopology(
            MPI.COMM_SELF,
            -1.0:2.0:1.0,
        )
        grid = ClimateMachine.Mesh.Grids.DiscontinuousSpectralElementGrid(
            topology;
            polynomialorder = N,
            FloatType = T,
            DeviceArray = Array,
        )

        ξ = ClimateMachine.Mesh.Grids.referencepoints(grid)
        ξ1 = ξ[1]
        ξ2 = ξ[2]
        a1, b1 = GaussQuadrature.legendre_coefs(T, N[1])
        a2, b2 = GaussQuadrature.legendre_coefs(T, N[2])
        V1 = GaussQuadrature.orthonormal_poly(ξ1, a1, b1)
        V2 = GaussQuadrature.orthonormal_poly(ξ2, a2, b2)

        Σ1 = ones(T, N[1] + 1)
        Σ2 = ones(T, N[2] + 1)
        Σ1[(Nc[1]:N[1]) .+ 1] .= 0
        Σ2[(Nc[2]:N[2]) .+ 1] .= 0

        W1 = V1 * Diagonal(Σ1) / V1
        W2 = V2 * Diagonal(Σ2) / V2

        filter =
            ClimateMachine.Mesh.Filters.MassPreservingCutoffFilter(grid, Nc)
        @test filter.filter_matrices[1] ≈ W1
        @test filter.filter_matrices[2] ≈ W2
    end

end

struct FilterTestModel{N} <: ClimateMachine.BalanceLaws.BalanceLaw end
ClimateMachine.BalanceLaws.vars_state(::FilterTestModel, ::Auxiliary, FT) =
    @vars()
ClimateMachine.BalanceLaws.init_state_auxiliary!(::FilterTestModel, _...) =
    nothing

# Legendre Polynomials
l0(r) = 1
l1(r) = r
l2(r) = (3 * r^2 - 1) / 2
l3(r) = (5 * r^3 - 3r) / 2

low(x, y, z) = l0(x) * l0(y) + 4 * l1(x) * l1(y) + 5 * l1(z) + 6 * l1(z) * l1(x)

high(x, y, z) = l2(x) * l3(y) + l3(x) + l2(y) + l3(z) * l1(y)

filtered(::EveryDirection, dim, x, y, z) = high(x, y, z)
filtered(::VerticalDirection, dim, x, y, z) =
    (dim == 2) ? l2(x) * l3(y) + l2(y) : l3(z) * l1(y)
filtered(::HorizontalDirection, dim, x, y, z) =
    (dim == 2) ? l2(x) * l3(y) + l3(x) : l2(x) * l3(y) + l3(x) + l2(y)

ClimateMachine.BalanceLaws.vars_state(::FilterTestModel{4}, ::Prognostic, FT) =
    @vars(q1::FT, q2::FT, q3::FT, q4::FT)
function ClimateMachine.BalanceLaws.init_state_prognostic!(
    ::FilterTestModel{4},
    state::Vars,
    aux::Vars,
    localgeo,
    filter_direction,
    dim,
)
    (x, y, z) = localgeo.coord
    state.q1 = low(x, y, z) + high(x, y, z)
    state.q2 = low(x, y, z) + high(x, y, z)
    state.q3 = low(x, y, z) + high(x, y, z)
    state.q4 = low(x, y, z) + high(x, y, z)

    if !isnothing(filter_direction)
        state.q1 -= filtered(filter_direction, dim, x, y, z)
        state.q3 -= filtered(filter_direction, dim, x, y, z)
    end
end

@testset "Exponential and Cutoff filter application" begin
    N = 3
    Ne = (1, 1, 1)

    @testset for FT in (Float64, Float32)
        @testset for dim in 2:3
            @testset for direction in (
                EveryDirection,
                HorizontalDirection,
                VerticalDirection,
            )
                brickrange = ntuple(
                    j -> range(FT(-1); length = Ne[j] + 1, stop = 1),
                    dim,
                )
                topl = ClimateMachine.Mesh.Topologies.BrickTopology(
                    MPI.COMM_WORLD,
                    brickrange,
                    periodicity = ntuple(j -> true, dim),
                )

                grid =
                    ClimateMachine.Mesh.Grids.DiscontinuousSpectralElementGrid(
                        topl,
                        FloatType = FT,
                        DeviceArray = ClimateMachine.array_type(),
                        polynomialorder = N,
                    )

                filter = ClimateMachine.Mesh.Filters.CutoffFilter(grid, 2)

                model = FilterTestModel{4}()
                dg = ClimateMachine.DGMethods.DGModel(
                    model,
                    grid,
                    nothing,
                    nothing,
                    nothing;
                    state_gradient_flux = nothing,
                )

                @testset for target in ((1, 3), (:q1, :q3))
                    Q = ClimateMachine.DGMethods.init_ode_state(
                        dg,
                        nothing,
                        dim,
                    )
                    ClimateMachine.Mesh.Filters.apply!(
                        Q,
                        target,
                        grid,
                        filter,
                        direction = direction(),
                    )
                    P = ClimateMachine.DGMethods.init_ode_state(
                        dg,
                        direction(),
                        dim,
                    )
                    @test Array(Q.data) ≈ Array(P.data)
                end
            end
        end
    end
end

@testset "Mass Preserving Cutoff filter application" begin
    N = 3
    Ne = (1, 1, 1)

    @testset for FT in (Float64, Float32)
        @testset for dim in 2:3
            @testset for direction in (
                EveryDirection,
                HorizontalDirection,
                VerticalDirection,
            )
                brickrange = ntuple(
                    j -> range(FT(-1); length = Ne[j] + 1, stop = 1),
                    dim,
                )
                topl = ClimateMachine.Mesh.Topologies.BrickTopology(
                    MPI.COMM_WORLD,
                    brickrange,
                    periodicity = ntuple(j -> true, dim),
                )

                grid =
                    ClimateMachine.Mesh.Grids.DiscontinuousSpectralElementGrid(
                        topl,
                        FloatType = FT,
                        DeviceArray = ClimateMachine.array_type(),
                        polynomialorder = N,
                    )

                filter = ClimateMachine.Mesh.Filters.MassPreservingCutoffFilter(
                    grid,
                    2,
                )

                model = FilterTestModel{4}()
                dg = ClimateMachine.DGMethods.DGModel(
                    model,
                    grid,
                    nothing,
                    nothing,
                    nothing;
                    state_gradient_flux = nothing,
                )

                @testset for target in ((1, 3), (:q1, :q3))
                    Q = ClimateMachine.DGMethods.init_ode_state(
                        dg,
                        nothing,
                        dim,
                    )
                    ClimateMachine.Mesh.Filters.apply!(
                        Q,
                        target,
                        grid,
                        filter,
                        direction = direction(),
                    )
                    P = ClimateMachine.DGMethods.init_ode_state(
                        dg,
                        direction(),
                        dim,
                    )
                    @test Array(Q.data) ≈ Array(P.data)
                end
            end
        end
    end
end

ClimateMachine.BalanceLaws.vars_state(
    ::FilterTestModel{1},
    ::Prognostic,
    FT,
) where {N} = @vars(q::FT)
function ClimateMachine.BalanceLaws.init_state_prognostic!(
    ::FilterTestModel{1},
    state::Vars,
    aux::Vars,
    localgeo,
)
    (x, y, z) = localgeo.coord
    state.q = abs(x) - 0.1
end

@testset "TMAR filter application" begin

    N = 4
    Ne = (2, 2, 2)

    @testset for FT in (Float64, Float32)
        @testset for dim in 2:3
            brickrange =
                ntuple(j -> range(FT(-1); length = Ne[j] + 1, stop = 1), dim)
            topl = ClimateMachine.Mesh.Topologies.BrickTopology(
                MPI.COMM_WORLD,
                brickrange,
                periodicity = ntuple(j -> true, dim),
            )

            grid = ClimateMachine.Mesh.Grids.DiscontinuousSpectralElementGrid(
                topl,
                FloatType = FT,
                DeviceArray = ClimateMachine.array_type(),
                polynomialorder = N,
            )

            model = FilterTestModel{1}()
            dg = ClimateMachine.DGMethods.DGModel(
                model,
                grid,
                nothing,
                nothing,
                nothing;
                state_gradient_flux = nothing,
            )

            @testset for target in ((1,), (:q,), :)
                Q = ClimateMachine.DGMethods.init_ode_state(dg)

                initialsumQ = weightedsum(Q)
                @test minimum(Q.realdata) < 0

                ClimateMachine.Mesh.Filters.apply!(
                    Q,
                    target,
                    grid,
                    ClimateMachine.Mesh.Filters.TMARFilter(),
                )

                sumQ = weightedsum(Q)

                @test minimum(Q.realdata) >= 0
                @test isapprox(initialsumQ, sumQ; rtol = 10 * eps(FT))
            end
        end
    end
end

function cubedshellwarp(a, b, c, R = max(abs(a), abs(b), abs(c)))

    function f(sR, ξ, η)
        X, Y = tan(π * ξ / 4), tan(π * η / 4)
        x1 = sR / sqrt(X^2 + Y^2 + 1)
        x2, x3 = X * x1, Y * x1
        x1, x2, x3
    end

    fdim = argmax(abs.((a, b, c)))
    if fdim == 1 && a < 0
        # (-R, *, *) : Face I from Ronchi, Iacono, Paolucci (1996)
        x1, x2, x3 = f(-R, b / a, c / a)
    elseif fdim == 2 && b < 0
        # ( *,-R, *) : Face II from Ronchi, Iacono, Paolucci (1996)
        x2, x1, x3 = f(-R, a / b, c / b)
    elseif fdim == 1 && a > 0
        # ( R, *, *) : Face III from Ronchi, Iacono, Paolucci (1996)
        x1, x2, x3 = f(R, b / a, c / a)
    elseif fdim == 2 && b > 0
        # ( *, R, *) : Face IV from Ronchi, Iacono, Paolucci (1996)
        x2, x1, x3 = f(R, a / b, c / b)
    elseif fdim == 3 && c > 0
        # ( *, *, R) : Face V from Ronchi, Iacono, Paolucci (1996)
        x3, x2, x1 = f(R, b / c, a / c)
    elseif fdim == 3 && c < 0
        # ( *, *,-R) : Face VI from Ronchi, Iacono, Paolucci (1996)
        x3, x2, x1 = f(-R, b / c, a / c)
    else
        error("invalid case for cubedshellwarp: $a, $b, $c")
    end

    return x1, x2, x3
end

@testset "Mass Preserving Cutoff Filter Conservation Test" begin
    N = 3
    Ne = (1, 1, 1)
    dim = 3
    # dim, direction
    @testset for FT in (Float64,)
        Rrange = [FT(1.0), FT(1.2)]

        topl = ClimateMachine.Mesh.Grids.StackedCubedSphereTopology(
            MPI.COMM_WORLD,
            1,
            Rrange,
            boundary = (5, 6),
        )

        grid = ClimateMachine.Mesh.Grids.DiscontinuousSpectralElementGrid(
            topl,
            FloatType = FT,
            DeviceArray = ClimateMachine.array_type(),
            polynomialorder = (N, N),
            meshwarp = cubedshellwarp,
        )

        mp_filter =
            ClimateMachine.Mesh.Filters.MassPreservingCutoffFilter(grid, 2)

        reg_filter = ClimateMachine.Mesh.Filters.CutoffFilter(grid, 2)

        model = FilterTestModel{4}()
        dg = ClimateMachine.DGMethods.DGModel(
            model,
            grid,
            nothing,
            nothing,
            nothing;
            state_gradient_flux = nothing,
        )

        # test mp filter
        filter = mp_filter
        Q = ClimateMachine.DGMethods.init_ode_state(dg, nothing, dim)
        sum_before_1 = weightedsum(Q, 1)
        sum_before_2 = weightedsum(Q, 2)
        sum_before_3 = weightedsum(Q, 3)
        target = 1:3
        ClimateMachine.Mesh.Filters.apply!(Q, target, grid, filter)
        sum_after_1 = weightedsum(Q, 1)
        sum_after_2 = weightedsum(Q, 2)
        sum_after_3 = weightedsum(Q, 3)

        @test sum_before_1 ≈ sum_after_1
        @test sum_before_2 ≈ sum_after_2
        @test sum_before_3 ≈ sum_after_3

        # test regular filter
        filter = reg_filter
        Q = ClimateMachine.DGMethods.init_ode_state(dg, nothing, dim)
        sum_before_1 = weightedsum(Q, 1)
        sum_before_2 = weightedsum(Q, 2)
        sum_before_3 = weightedsum(Q, 3)
        target = 1:3
        ClimateMachine.Mesh.Filters.apply!(Q, target, grid, filter)
        sum_after_1 = weightedsum(Q, 1)
        sum_after_2 = weightedsum(Q, 2)
        sum_after_3 = weightedsum(Q, 3)

        @test !(sum_before_1 ≈ sum_after_1)
        @test !(sum_before_2 ≈ sum_after_2)
        @test !(sum_before_3 ≈ sum_after_3)

    end
end


text:
include("../utils.jl")

@testset "Symmetric Positive Semidefinite Matrices of Fixed Rank" begin
    @testset "Real Matrices" begin
        M = SymmetricPositiveSemidefiniteFixedRank(4, 2)
        @test repr(M) == "SymmetricPositiveSemidefiniteFixedRank(4, 2, ℝ)"
        @test manifold_dimension(M) == 7
        q = [1.0 0.0; 0.0 1.0; 0.0 0.0; 0.0 0.0]
        @test is_point(M, q)
        Y = [1.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]
        @test_throws DomainError is_point(M, Y, true)
        @test is_vector(M, q, Y)
        q2 = [2.0 1.0; 0.0 0.0; 0.0 1.0; 0.0 0.0]
        q3 = [0.0 0.0; 1.0 0.0; 0.0 1.0; 0.0 0.0]
        X = log(M, q, q2)
        X3 = vector_transport_to(M, q, X, q3, ProjectionTransport())
        X3t = project(M, q3, X)
        @test is_vector(M, q3, X3)
        @test isapprox(M, q3, X3, X3t)

        types = [Matrix{Float64}]
        TEST_FLOAT32 && push!(types, Matrix{Float32})
        TEST_STATIC_SIZED && push!(types, MMatrix{4,2,Float64,8})
        for T in types
            pts = [convert(T, q), convert(T, q2), convert(T, q3)]
            @testset "Type $T" begin
                test_manifold(
                    M,
                    pts,
                    exp_log_atol_multiplier=5,
                    is_tangent_atol_multiplier=5,
                    test_project_tangent=true,
                    test_inplace=true,
                )
            end
        end
    end
    @testset "Complex Matrices" begin
        M = SymmetricPositiveSemidefiniteFixedRank(4, 2, ℂ)
        @test repr(M) == "SymmetricPositiveSemidefiniteFixedRank(4, 2, ℂ)"
        @test manifold_dimension(M) == 12
    end
end


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function subroutine(🐢::Turtle)
    for i in 1:10
        Forward(🐢, 40)
        Turn(🐢, 51)
        Circle(🐢, 1)
    end
end

function test_turtles(fname)
    Drawing(1200, 1200, fname)
    origin()
    background("black")

    # let's have two turtles
    raphael = Turtle(0, 0, true, 0, (1.0, 0.25, 0.25))
    michaelangelo = Turtle(0, 0, true, 0, (1.0, 0.25, 0.25))

    Pencolor(raphael, 1.0, 0.4, 0.2);       Pencolor(michaelangelo, 0.2, 0.9, 1.0)
    Penwidth(raphael, 0.75);                Penwidth(michaelangelo, 1.5)
    Reposition(raphael, -400, -300);        Reposition(michaelangelo, 0, 200)

    Orientation(raphael, pi);               Orientation(michaelangelo, -pi)

    for 🐢 in [raphael, michaelangelo]
        setlinecap("round")
        setopacity(0.5)
        Penwidth(🐢, .5)
        n = 400
        Penup(🐢)
        Pendown(🐢)
        for i in 1:200
            Forward(🐢, n)
            Turn(🐢, 29)
            🐢 == raphael && HueShift(🐢, 1)
            n += 1.25
            subroutine(🐢)
        end
        Reposition(🐢, O)
        Rectangle(🐢, 10, 10)
    end

    fontsize(30)
    Reposition(raphael, Point(400, -100));   Reposition(michaelangelo, 400, -50)
    Message(raphael, "Raphael");      Message(michaelangelo, "Michaelangelo")

    Push(raphael)
    Pop(raphael)

    Pencolor(raphael, "white")
    Randomize_saturation(michaelangelo)

    # test warnings
    @info("stack was already empty warning expected...")
    Pop(michaelangelo) # should get "stack empty" warning!
    @test finish() == true
end

fname = "turtles-all-the-way-down.png"
test_turtles(fname)
println("...finished test: output in $(fname)")