text:
using FiniteDifferences, Test
using FiniteDifferences: jvp, j′vp
using ChainRules
using ChainRulesCore: AbstractDifferential

const _fdm = central_fdm(5, 1)

"""
    test_scalar(f, x; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), test_wirtinger=x isa Complex, kwargs...)

Given a function `f` with scalar input an scalar output, perform finite differencing checks,
at input point `x` to confirm that there are correct ChainRules provided.

# Arguments
- `f`: Function for which the `frule` and `rrule` should be tested.
- `x`: input at which to evaluate `f` (should generally be set to an arbitary point in the domain).

- `test_wirtinger`: test whether the wirtinger derivative is correct, too

All keyword arguments except for `fdm` and `test_wirtinger` are passed to `isapprox`.
"""
function test_scalar(f, x; rtol=1e-9, atol=1e-9, fdm=_fdm, test_wirtinger=x isa Complex, kwargs...)
    ensure_not_running_on_functor(f, "test_scalar")

    @testset "$f at $x, $(nameof(rule))" for rule in (rrule, frule)
        res = rule(f, x)
        @test res !== nothing  # Check the rule was defined
        fx,  prop_rule = res
        @test fx == f(x)  # Check we still get the normal value, right

        if rule == rrule
            ∂self, ∂x = prop_rule(1)
            @test ∂self === NO_FIELDS
        else # rule == frule
            # Got to input extra first aguement for internals
            # But it is only a dummy since this is not a functor
            ∂x, = prop_rule(NamedTuple(), 1)
        end

        # Check that we get the derivative right:
        if !test_wirtinger
            @test isapprox(
                ∂x, fdm(f, x);
                rtol=rtol, atol=atol, kwargs...
            )
        else
            # For complex arguments, also check if the wirtinger derivative is correct
            ∂Re = fdm(ϵ -> f(x + ϵ), 0)
            ∂Im = fdm(ϵ -> f(x + im*ϵ), 0)
            ∂ = 0.5(∂Re - im*∂Im)
            ∂̅ = 0.5(∂Re + im*∂Im)
            @test isapprox(
                wirtinger_primal(∂x), ∂;
                rtol=rtol, atol=atol, kwargs...
            )
            @test isapprox(
                wirtinger_conjugate(∂x), ∂̅;
                rtol=rtol, atol=atol, kwargs...
            )
        end
    end
end

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

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

"""
    frule_test(f, (x, ẋ)...; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), kwargs...)

# Arguments
- `f`: Function for which the `frule` should be tested.
- `x`: input at which to evaluate `f` (should generally be set to an arbitary point in the domain).
- `ẋ`: differential w.r.t. `x` (should generally be set randomly).

All keyword arguments except for `fdm` are passed to `isapprox`.
"""
function frule_test(f, (x, ẋ); rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    return frule_test(f, ((x, ẋ),); rtol=rtol, atol=atol, fdm=fdm, kwargs...)
end

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

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

"""
    rrule_test(f, ȳ, (x, x̄)...; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), kwargs...)

# Arguments
- `f`: Function to which rule should be applied.
- `ȳ`: adjoint w.r.t. output of `f` (should generally be set randomly).
- `x`: input at which to evaluate `f` (should generally be set to an arbitary point in the domain).
- `x̄`: currently accumulated adjoint (should generally be set randomly).

All keyword arguments except for `fdm` are passed to `isapprox`.
"""
function rrule_test(f, ȳ, (x, x̄)::Tuple{Any, Any}; rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    ensure_not_running_on_functor(f, "rrule_test")

    # Check correctness of evaluation.
    fx, pullback = ChainRules.rrule(f, x)
    @test fx ≈ f(x)
    (∂self, x̄_ad) = pullback(ȳ)
    @test ∂self === NO_FIELDS  # No internal fields
    # Correctness testing via finite differencing.
    x̄_fd = j′vp(fdm, f, ȳ, x)
    @test isapprox(x̄_ad, x̄_fd; rtol=rtol, atol=atol, kwargs...)

    # Assuming x̄_ad to be correct, check that other ChainRules mechanisms are correct.
    test_accumulation(x̄, x̄_ad)
    test_accumulation(Zero(), x̄_ad)
end

function _make_fdm_call(fdm, f, ȳ, xs, ignores)
    sig = Expr(:tuple)
    call = Expr(:call, f)
    newxs = Any[]
    arginds = Int[]
    i = 1
    for (x, ignore) in zip(xs, ignores)
        if ignore
            push!(call.args, x)
        else
            push!(call.args, Symbol(:x, i))
            push!(sig.args, Symbol(:x, i))
            push!(newxs, x)
            push!(arginds, i)
        end
        i += 1
    end
    fdexpr = :(j′vp($fdm, $sig -> $call, $ȳ, $(newxs...)))
    fd = eval(fdexpr)
    fd isa Tuple || (fd = (fd,))
    args = Any[nothing for _ in 1:length(xs)]
    for (dx, ind) in zip(fd, arginds)
        args[ind] = dx
    end
    return (args...,)
end

function rrule_test(f, ȳ, xx̄s::Tuple{Any, Any}...; rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    ensure_not_running_on_functor(f, "rrule_test")

    # Check correctness of evaluation.
    xs, x̄s = collect(zip(xx̄s...))
    y, pullback = rrule(f, xs...)
    @test f(xs...) == y

    @assert !(isa(ȳ, Thunk))
    ∂s = pullback(ȳ)
    ∂self = ∂s[1]
    x̄s_ad = ∂s[2:end]
    @test ∂self === NO_FIELDS

    # Correctness testing via finite differencing.
    x̄s_fd = _make_fdm_call(fdm, f, ȳ, xs, x̄s .== nothing)
    for (x̄_ad, x̄_fd) in zip(x̄s_ad, x̄s_fd)
        if x̄_fd === nothing
            # The way we've structured the above, this tests that the rule is a DNERule
            @test x̄_ad isa DNE
        else
            @test isapprox(x̄_ad, x̄_fd; rtol=rtol, atol=atol, kwargs...)
        end
    end

    # Assuming the above to be correct, check that other ChainRules mechanisms are correct.
    for (x̄, x̄_ad) in zip(x̄s, x̄s_ad)
        x̄ === nothing && continue
        test_accumulation(x̄, x̄_ad)
        test_accumulation(Zero(), x̄_ad)
    end
end

function Base.isapprox(ad::Wirtinger, fd; kwargs...)
    error("Finite differencing with Wirtinger rules not implemented")
end
function Base.isapprox(d_ad::Casted, d_fd; kwargs...)
    return all(isapprox.(extern(d_ad), d_fd; kwargs...))
end
function Base.isapprox(d_ad::DNE, d_fd; kwargs...)
    error("Tried to differentiate w.r.t. a DNE")
end
function Base.isapprox(d_ad::AbstractDifferential, d_fd; kwargs...)
    return isapprox(extern(d_ad), d_fd; kwargs...)
end

function test_accumulation(x̄, ∂x)
    @test all(extern(x̄ + ∂x) .≈ extern(x̄) .+ extern(∂x))
    test_accumulate(x̄, ∂x)
    test_accumulate!(x̄, ∂x)
    test_store!(x̄, ∂x)
end

function test_accumulate(x̄::Zero, ∂x)
    @test extern(accumulate(x̄, ∂x)) ≈ extern(∂x)
end

function test_accumulate(x̄::Number, ∂x)
    @test extern(accumulate(x̄, ∂x)) ≈ extern(x̄) + extern(∂x)
end

function test_accumulate(x̄::AbstractArray, ∂x)
    x̄_old = copy(x̄)
    @test all(extern(accumulate(x̄, ∂x)) .≈ (extern(x̄) .+ extern(∂x)))
    @test x̄ == x̄_old  # make sure didn't mutate x̄
end

test_accumulate!(x̄::Zero, ∂x) = nothing

function test_accumulate!(x̄::Number, ∂x)
    # This case won't have been inplace as `Number` is immutable
    @test accumulate!(x̄, ∂x) ≈ accumulate(x̄, ∂x)
end

function test_accumulate!(x̄::AbstractArray, ∂x)
    x̄_copy = copy(x̄)

    accumulate!(x̄_copy, ∂x)  # this should have actually been in-place
    @test extern(x̄_copy) ≈ (extern(x̄) .+ extern(∂x))
end

test_store!(x̄::Zero, ∂x) = nothing
test_store!(x̄::Number, ∂x) = nothing

function test_store!(x̄::AbstractArray, ∂x)
    x̄_store = copy(x̄)
    store!(x̄_store, ∂x)
    @test x̄_store ≈ extern(∂x)

    # store! is the same as `accumulate!` to a zero array
    x̄_acc = false.*x̄
    accumulate!(x̄_acc, ∂x)
    @test x̄_acc ≈ x̄_store
end


text:
using YaoBlocks
using LuxurySparse
using LinearAlgebra
using ..SymEngine
using ..SymEngine: BasicType, BasicOp, BasicTrigFunction

op_types = [:Mul, :Add, :Pow]
const BiVarOp = Union{[SymEngine.BasicType{Val{i}} for i in op_types]...}

export @vars

simag = SymFunction("Im")
sreal = SymFunction("Re")
sabs = SymFunction("abs")

Base.promote_rule(::Type{Bool}, ::Type{Basic}) = Basic
Base.conj(x::Basic) = Basic(conj(SymEngine.BasicType(x)))
Base.conj(x::BasicType) = real(x) - im * imag(x)
Base.conj(x::BiVarOp) = juliafunc(x)(conj.(get_args(x.x))...)
Base.conj(x::BasicTrigFunction) = juliafunc(x)(conj.(get_args(x.x)...)...)
# WARNING: symbols and constants are assumed real!
Base.imag(x::BasicType{Val{:Constant}}) = Basic(0)
Base.imag(x::BasicType{Val{:Symbol}}) = Basic(0)
Base.abs(x::Basic) = Basic(abs(SymEngine.BasicType(x)))
Base.abs(x::BasicType{Val{:Constant}}) = x
Base.abs(x::BasicType{Val{:Symbol}}) = x
function Base.imag(x::BasicType{Val{:Add}})
    args = get_args(x.x)
    mapreduce(imag, +, args)
end

function Base.real(x::BasicType{Val{:Add}})
    args = get_args(x.x)
    mapreduce(real, +, args)
end

function Base.abs(x::BasicType{Val{:Add}})
    args = get_args(x.x)
    mapreduce(abs, +, args)
end

function Base.imag(x::BasicType{Val{:Mul}})
    args = (get_args(x.x)...,)
    get_mul_imag(args)
end

function Base.real(x::BasicType{Val{:Pow}})
    a, b = get_args(x.x)
    if imag(a) == 0 && imag(b) == 0
        return x.x
    else
        if imag(a) == 0
            return a^real(b) * cos(log(a) * imag(b))
        else
            return sreal(x.x)
        end
    end
end

function Base.imag(x::BasicType{Val{:Pow}})
    a, b = get_args(x.x)
    if imag(a) == 0 && imag(b) == 0
        return Basic(0)
    else
        if imag(a) == 0
            return a^real(b) * sin(log(a) * imag(b))
        else
            return simag(x.x)
        end
    end
end

function Base.abs(x::BasicType{Val{:Pow}})
    a, b = get_args(x.x)
    abs(a)^real(b)
end

function Base.real(x::BasicType{Val{:Mul}})
    args = (get_args(x.x)...,)
    get_mul_real(args)
end

function get_mul_imag(args::NTuple{N,Any}) where {N}
    imag(args[1]) * get_mul_real(args[2:end]) + real(args[1]) * get_mul_imag(args[2:end])
end
get_mul_imag(args::Tuple{Basic}) = imag(args[1])

function get_mul_real(args::NTuple{N,Any}) where {N}
    real(args[1]) * get_mul_real(args[2:end]) - imag(args[1]) * get_mul_imag(args[2:end])
end
get_mul_real(args::Tuple{Basic}) = real(args[1])

function Base.real(x::BasicTrigFunction)
    a, = get_args(x.x)
    if imag(a) == 0
        return x.x
    else
        return sreal(x.x)
    end
end

function Base.abs(x::BasicTrigFunction)
    a, = get_args(x.x)
    if imag(a) == 0
        return x.x
    else
        return sabs(x.x)
    end
end

function Base.imag(x::BasicTrigFunction)
    a, = get_args(x.x)
    if imag(a) == 0
        return Basic(0)
    else
        return simag(x.x)
    end
end

@generated function juliafunc(x::BasicType{Val{T}}) where {T}
    SymEngine.map_fn(T, SymEngine.fn_map)
end

const SymReal = Union{Basic,SymEngine.BasicRealNumber}
YaoBlocks.RotationGate(block::GT, theta::T) where {D,T<:SymReal,GT<:AbstractBlock{D}} =
    RotationGate{D,T,GT}(block, theta)

YaoBlocks.phase(θ::SymReal) = PhaseGate(θ)
YaoBlocks.shift(θ::SymReal) = ShiftGate(θ)

YaoBlocks.mat(::Type{Basic}, gate::GT) where GT<:ConstantGate = _pretty_basic.(mat(gate))
YaoBlocks.mat(::Type{Basic}, gate::ConstGate.TGate) = Diagonal(Basic[1, exp(Basic(im)*Basic(π)/4)])
YaoBlocks.mat(::Type{Basic}, gate::ConstGate.TdagGate) = Diagonal(Basic[1, exp(-Basic(im)*Basic(π)/4)])
YaoBlocks.mat(::Type{Basic}, ::HGate) = 1 / sqrt(Basic(2)) * Basic[1 1; 1 -1]
YaoBlocks.mat(::Type{Basic}, gate::ShiftGate) = Diagonal([1, exp(im * gate.theta)])
YaoBlocks.mat(::Type{Basic}, gate::PhaseGate) = exp(im * gate.theta) * IMatrix{2}()
function YaoBlocks.mat(::Type{Basic}, R::RotationGate{D}) where {D}
    I = IMatrix{D^nqudits(R)}()
    return I * cos(R.theta / 2) - im * sin(R.theta / 2) * mat(Basic, R.block)
end
for GT in [:XGate, :YGate, :ZGate]
    @eval YaoBlocks.mat(::Type{Basic}, R::RotationGate{D,T,<:$GT}) where {D,T} =
        invoke(mat, Tuple{Type{Basic},RotationGate}, Basic, R)
end

for T in [:(RotationGate{D,<:SymReal} where D), :(PhaseGate{<:SymReal}), :(ShiftGate{<:SymReal})]
    @eval YaoBlocks.mat(gate::$T) = mat(Basic, gate)
end

YaoBlocks.PSwap(n::Int, locs::Tuple{Int,Int}, θ::SymReal) =
    YaoBlocks.PutBlock(n, rot(ConstGate.SWAPGate(), θ), locs)

YaoBlocks.pswap(n::Int, i::Int, j::Int, α::SymReal) = PSwap(n, (i, j), α)
YaoBlocks.pswap(i::Int, j::Int, α::SymReal) = n -> pswap(n, i, j, α)

export subs
SymEngine.subs(c::AbstractBlock, args...; kwargs...) = subs(Basic, c, args...; kwargs...)
function SymEngine.subs(::Type{T}, c::AbstractBlock, args...; kwargs...) where {T}
    c = setiparams(c, map(x -> T(subs(x, args...; kwargs...)), getiparams(c))...)
    chsubblocks(c, [subs(T, blk, args..., kwargs...) for blk in subblocks(c)])
end

# dumpload
YaoBlocks.tokenize_param(param::Basic) = QuoteNode(Symbol(param))
YaoBlocks.parse_param(x::QuoteNode) = :(Basic($x))

ext: jl


text:
module TOKCC

points(weight, repetitions) = (weight/8) ^ 1.35 * repetitions

score(weight, repetitions, bodyweight) = points(weight, repetitions) / (bodyweight ^ (2/3))

end

max_stars_repo_name: wjholden/TOKCC

max_issues_repo_name: wjholden/TOKCC

max_forks_repo_name: wjholden/TOKCC


text:
using Test

include("../src/types.jl")
include("../src/intersecttest.jl")

@testset "address error from focus point 1.5, 1.5, 1.5 with single triangle mesh" begin

lx = [-100663297., 100663297., 100663297., -67108864.]
ly = [-1., 1., 0., 0.]
lc = [150994945., -150994946., -150994944., 100663296.]

@test ccrY0(lx, ly, lc, Dim(4)) == []
@test ccr(lx, ly, lc, Dim(4)) == []

end

@testset "address problem with first simplex for 3D sphere mesh" begin

  lx = [0.5765708843703083, 0.21828465814030382, -1.5049122768452694, 0.7100567343346572]
  ly = [-0.9714406736017811, -0.7988509382303219, 0.9398460751671794, 0.8304455366649236]
  lc = [0.0, -0.0, -0.0, 1.0]

  # FIXME: need to use approx and instance for rtoldefault
  @test_broken ccrXY(lx, ly, lc, Dim(4)) ≈ [Intersect2D(-0.8312942, -0.4933905, -0.4902460, -0.7849979)]
  @test_broken ccr(lx, ly, lc, Dim(4)) ≈ [Intersect2D(-0.8312942, -0.4933905, -0.4902460, -0.7849979)]
end


text:
export
    WindowCreator,
    create_window,
    SyncWindow,
    AsyncWindow,
    CoherentWindow

"""
    abstract type WindowCreator end

A type that identifies how to create windows of instances for filtering algorithms.

T represents the time type.

Must implement `create_window`.
"""

abstract type WindowCreator{T} end

struct SyncWindow <: WindowCreator{Int} end

"""
    create_window(::SyncWindow, runtime::Runtime, variables::Vector{<:Variable}, time::Int)::Vector{Instance}

Creates a window by instantiating all variables at all intermediate time steps from the earliest parent to the given time.
The `variables` argument is ignored.
"""
function create_window(::SyncWindow, runtime::Runtime, variables_to_sample::Vector{<:Variable}, time::Int)::Vector{Instance}
    prevtime = time
    net = get_network(runtime)
    ord = topsort(get_transition_graph(net))
    for v in variables_to_sample
        pars = get_transition_parents(net, v)
        for p in pars
            t = get_time(latest_instance_before(runtime, p, time, p != v))
            prevtime = min(prevtime, t)
        end
    end

    insts = Instance[]
    for n in ord
        ph = Placeholder{output_type(n)}(n.name)
        push!(insts, PlaceholderInstance(ph, prevtime))
    end
    for t in prevtime+1:time
        for n in ord
            push!(insts, ensure_instance!(runtime, n, t))
        end
    end
    return insts
end

struct AsyncWindow{T <: Number} <: WindowCreator{T} end

"""
    create_window(::AsyncWindow, runtime::Runtime, variables::Vector{<:Variable}, time::Int)::Vector{Instance}

Creates a window by instantiating only the given variables at the given time.
"""

function create_window(::AsyncWindow{T}, runtime::Runtime, variables::Vector{<:Variable}, time::T)::Vector{Instance} where T
    insts = Instance[]
    done = Set{Variable}()
    for v in variables
        for p in get_transition_parents(get_network(runtime), v)
            if !(p in done)
                parinst = latest_instance_before(runtime, p, time, true) # changed false to true!!!
                partime = get_time(parinst)
                ph = Placeholder{output_type(p)}(p.name)
                phinst = PlaceholderInstance(ph, partime)
                push!(insts, phinst)
                push!(done, p)
            end
        end
        push!(insts, ensure_instance!(runtime, v, time))
        push!(done, v)
    end
    return insts
end

"""
    struct CoherentWindow <: WindowCreator end

A variant of AsyncWindow that ensures that parent values are never stale for any variable that
gets updated in a filter step. In other words, if any parent of a direct parent has been updated more recently than a variable
to be updated, the direct parent is added to the variables to be updated. This condition then recurses for the direct parents.
"""
struct CoherentWindow{T <: Number} <: WindowCreator{T} end

function create_window(::CoherentWindow{T}, runtime::Runtime, variables::Vector{<:Variable}, time::T)::Vector{Instance} where T
    # Note: This method does not allow placeholder parents of dynamic variables
    net = get_network(runtime)
    parents = get_transition_graph(net)
    fullvars = Set{Variable}()
    order = topsort(parents)
    times = Dict([(n, get_time(current_instance(runtime, n))) for n in order])

    function ensure(v)
        if !(v in fullvars)
            push!(fullvars, v)
            for anc in ancestors(parents, v, Set{Node}())
                for grandanc in get(parents, anc, [])
                    if times[grandanc] > times[anc]
                        ensure(anc)
                        break
                    end
                end
            end
            times[v] = time # Need to do this to ensure instantiation through a chain of dependencies
        end
    end

    for var in variables
        ensure(var)
    end
    # We must make sure variables get sampled in the correct order to maintain coherence
    orderedvars = Variable[]
    for v in order
        if v in fullvars
            push!(orderedvars, v)
        end
    end

    create_window(AsyncWindow{T}(), runtime, orderedvars, time)
end


text:
using Makie, GeometryTypes, Colors

scene = Scene()
# this is just a work around that will be gone immediatly after I figure out how to best
# insert default cameras into the scene
r = range(0, stop=3, length=4)
a = axis(r, r)
center!(scene, 0.2)
cam = scene[:screen].cameras[:orthographic_pixel]
pos = lift_node(getindex.(scene, (:mouseposition, :time))...) do mpos, t
    map(range(0, stop=2pi, length=60)) do i
        circle = Point2f0(sin(i), cos(i))
        mouse = to_world(Point2f0(mpos), cam)
        secondary = (sin((i * 10f0) + t) * 0.09) * normalize(circle)
        (secondary .+ circle) .+ mouse
    end
end

p1 = lines(pos)
p2 = scatter(
    pos, markersize = 0.1f0,
    marker = :star5,
    color = lift_node(x-> x .+ 0.2, p1[:color])
)

p1[:color] = RGBA(1, 0, 0, 0.1)
p2[:marker] = 'π'
p2[:markersize] = 0.2
p2[:marker] = 'o'

for i = range(0.01, stop=0.4, length=100)
    p2[:markersize] = i
    yield()
    sleep(0.01)
end


text:
# This file includes code that was formerly a part of Julia. License is MIT: http://julialang.org/license

## Error messages for Unicode / UTF support

const UTF_ERR_SHORT             = "invalid UTF-8 sequence starting at index <<1>> (0x<<2>> missing one or more continuation bytes)"
const UTF_ERR_CONT              = "invalid UTF-8 sequence starting at index <<1>> (0x<<2>> is not a continuation byte)"
const UTF_ERR_LONG              = "invalid UTF-8 sequence, overlong encoding starting at index <<1>> (0x<<2>>)"
const UTF_ERR_NOT_LEAD          = "not a leading Unicode surrogate code unit at index <<1>> (0x<<2>>)"
const UTF_ERR_NOT_TRAIL         = "not a trailing Unicode surrogate code unit at index <<1>> (0x<<2>>)"
const UTF_ERR_NOT_SURROGATE     = "not a valid Unicode surrogate code unit at index <<1>> (0x<<2>>)"
const UTF_ERR_MISSING_SURROGATE = "missing trailing Unicode surrogate code unit after index <<1>> (0x<<2>>)"
const UTF_ERR_INVALID           = "invalid Unicode character starting at index <<1>> (0x<<2>> > 0x10ffff)"
const UTF_ERR_SURROGATE         = "surrogate encoding not allowed in UTF-8 or UTF-32, at index <<1>> (0x<<2>>)"
const UTF_ERR_NULL_16_TERMINATE = "UTF16String data must be NULL-terminated"
const UTF_ERR_NULL_32_TERMINATE = "UTF32String data must be NULL-terminated"
const UTF_ERR_ODD_BYTES_16      = "UTF16String can't have odd number of bytes <<1>>"
const UTF_ERR_ODD_BYTES_32      = "UTF32String must have multiple of 4 bytes <<1>>"
const UTF_ERR_INVALID_CHAR      = "invalid Unicode character (0x<<2>> > 0x10ffff)"
const UTF_ERR_INVALID_8         = "invalid UTF-8 data"
const UTF_ERR_INVALID_16        = "invalid UTF-16 data"
const UTF_ERR_INVALID_INDEX     = "invalid character index"
const UTF_ERR_MAP_CHAR          = "map(f,s::AbstractString) requires f to return Char; try map(f,collect(s)) or a comprehension instead"

## Functions to check validity of UTF-8, UTF-16, and UTF-32 encoded strings,
#  and also to return information necessary to convert to other encodings

is_surrogate_lead(c::Unsigned) = ((c & ~0x003ff) == 0xd800)
is_surrogate_trail(c::Unsigned) = ((c & ~0x003ff) == 0xdc00)
is_surrogate_codeunit(c::Unsigned) = ((c & ~0x007ff) == 0xd800)
is_valid_continuation(c) = ((c & 0xc0) == 0x80)

## Return flags for check_string function

const UTF_LONG = 1              ##< Long encodings are present
const UTF_LATIN1 = 2            ##< characters in range 0x80-0xFF present
const UTF_UNICODE2 = 4          ##< characters in range 0x100-0x7ff present
const UTF_UNICODE3 = 8          ##< characters in range 0x800-0xd7ff, 0xe000-0xffff
const UTF_UNICODE4 = 16         ##< non-BMP characters present
const UTF_SURROGATE = 32        ##< surrogate pairs present

# Get a UTF-8 continuation byte, give error if invalid, return updated character value
@inline function get_continuation(ch::UInt32, byt::UInt8, pos)
    if !is_valid_continuation(byt)
        throw(UnicodeError(UTF_ERR_CONT, pos, byt))
    end
    (ch << 6) | (byt & 0x3f)
end

"""
Validates and calculates number of characters in a UTF-8,UTF-16 or UTF-32 encoded vector/string

Warning: this function does not check the bounds of the start or end positions
Use `checkstring` to make sure the bounds are checked

Input Arguments:

* `dat`    UTF-8 (`Vector{UInt8}`), UTF-16 (`Vector{UInt16}`) or UTF-32 (`Vector{UInt32}`, `AbstractString`) encoded string

Optional Input Arguments:

* `pos`    start position (defaults to 1)
* `endpos` end position   (defaults to `lastindex(dat)`)

Keyword Arguments:

* `accept_long_null`  = `true`  # Modified UTF-8 (`\\0` represented as `b\"\\xc0\\x80\"`)
* `accept_surrogates` = `true`  # `CESU-8`
* `accept_long_char`  = `false` # Accept arbitrary long encodings

Returns:

* (total characters, flags, 4-byte, 3-byte, 2-byte)

Throws:

* `UnicodeError`
"""
function unsafe_checkstring end

function unsafe_checkstring(dat::AbstractVector{UInt8},
                      pos = 1,
                      endpos = length(dat)
                      ;
                      accept_long_null  = true,
                      accept_surrogates = true,
                      accept_long_char  = false)
    local byt::UInt8, ch::UInt32, surr::UInt32
    flags::UInt = 0
    totalchar = num2byte = num3byte = num4byte = 0
    @inbounds while pos <= endpos
        ch = dat[pos]
        pos += 1
        totalchar += 1
        if ch > 0x7f
            # Check UTF-8 encoding
            if ch < 0xe0
                # 2-byte UTF-8 sequence (i.e. characters 0x80-0x7ff)
                (pos > endpos) && throw(UnicodeError(UTF_ERR_SHORT, pos, ch))
                byt = dat[pos]
                pos += 1
                ch = get_continuation(ch & 0x3f, byt, pos)
                if ch > 0x7f
                    num2byte += 1
                    flags |= (ch > 0xff) ? UTF_UNICODE2 : UTF_LATIN1
                elseif accept_long_char
                    flags |= UTF_LONG
                elseif (ch == 0) && accept_long_null
                    flags |= UTF_LONG
                else
                    throw(UnicodeError(UTF_ERR_LONG, pos, ch))
                end
             elseif ch < 0xf0
                # 3-byte UTF-8 sequence (i.e. characters 0x800-0xffff)
                (pos + 1 > endpos) && throw(UnicodeError(UTF_ERR_SHORT, pos, ch))
                byt = dat[pos]
                pos += 1
                ch = get_continuation(ch & 0x0f, byt, pos)
                byt = dat[pos]
                pos += 1
                ch = get_continuation(ch, byt, pos)
                # check for surrogate pairs, make sure correct
                if is_surrogate_codeunit(ch)
                    !is_surrogate_lead(ch) && throw(UnicodeError(UTF_ERR_NOT_LEAD, pos-2, ch))
                    # next character *must* be a trailing surrogate character
                    (pos + 2 > endpos) && throw(UnicodeError(UTF_ERR_MISSING_SURROGATE, pos-2, ch))
                    byt = dat[pos]
                    pos += 1
                    (byt != 0xed) && throw(UnicodeError(UTF_ERR_NOT_TRAIL, pos, byt))
                    byt = dat[pos]
                    pos += 1
                    surr = get_continuation(0x0000d, byt, pos)
                    byt = dat[pos]
                    pos += 1
                    surr = get_continuation(surr, byt, pos)
                    !is_surrogate_trail(surr) && throw(UnicodeError(UTF_ERR_NOT_TRAIL, pos-2, surr))
                    !accept_surrogates && throw(UnicodeError(UTF_ERR_SURROGATE, pos-2, surr))
                    flags |= UTF_SURROGATE
                    num4byte += 1
                elseif ch > 0x07ff
                    num3byte += 1
                elseif accept_long_char
                    flags |= UTF_LONG
                    num2byte += 1
                else
                    throw(UnicodeError(UTF_ERR_LONG, pos-2, ch))
                end
            elseif ch < 0xf5
                # 4-byte UTF-8 sequence (i.e. characters > 0xffff)
                (pos + 2 > endpos) && throw(UnicodeError(UTF_ERR_SHORT, pos, ch))
                byt = dat[pos]
                pos += 1
                ch = get_continuation(ch & 0x07, byt, pos)
                byt = dat[pos]
                pos += 1
                ch = get_continuation(ch, byt, pos)
                byt = dat[pos]
                pos += 1
                ch = get_continuation(ch, byt, pos)
                if ch > 0x10ffff
                    throw(UnicodeError(UTF_ERR_INVALID, pos-3, ch))
                elseif ch > 0xffff
                    num4byte += 1
                elseif is_surrogate_codeunit(ch)
                    throw(UnicodeError(UTF_ERR_SURROGATE, pos-3, ch))
                elseif accept_long_char
                    # This is an overly long encoded character
                    flags |= UTF_LONG
                    if ch > 0x7ff
                        num3byte += 1
                    elseif ch > 0x7f
                        num2byte += 1
                    end
                else
                    throw(UnicodeError(UTF_ERR_LONG, pos-2, ch))
                end
            else
                throw(UnicodeError(UTF_ERR_INVALID, pos, ch))
            end
        end
    end
    num3byte != 0 && (flags |= UTF_UNICODE3)
    num4byte != 0 && (flags |= UTF_UNICODE4)
    return totalchar, flags, num4byte, num3byte, num2byte
end

AbstractString1632{Tel<:Union{UInt16,UInt32}} = Union{AbstractVector{Tel}, AbstractString}

function unsafe_checkstring(
                      dat::AbstractString1632,
                      pos = 1,
                      endpos = lastindex(dat)
                      ;
                      accept_long_null  = true,
                      accept_surrogates = true,
                      accept_long_char  = false)
    local ch::UInt32
    flags::UInt = 0
    totalchar = num2byte = num3byte = num4byte = 0
    @inbounds while pos <= endpos
        ch = dat[pos]
        pos = nextind(dat, pos)
        totalchar += 1
        if ch > 0x7f
            if ch < 0x100
                num2byte += 1
                flags |= UTF_LATIN1
            elseif ch < 0x800
                num2byte += 1
                flags |= UTF_UNICODE2
            elseif ch > 0x0ffff
                (ch > 0x10ffff) && throw(UnicodeError(UTF_ERR_INVALID, pos, ch))
                num4byte += 1
            elseif !is_surrogate_codeunit(ch)
                num3byte += 1
            elseif is_surrogate_lead(ch)
                pos > endpos && throw(UnicodeError(UTF_ERR_MISSING_SURROGATE, pos, ch))
                # next character *must* be a trailing surrogate character
                ch = dat[pos]
                pos = nextind(dat, pos)
                !is_surrogate_trail(ch) && throw(UnicodeError(UTF_ERR_NOT_TRAIL, pos, ch))
                num4byte += 1
                if !(typeof(dat) <: AbstractVector{UInt16})
                    !accept_surrogates && throw(UnicodeError(UTF_ERR_SURROGATE, pos, ch))
                    flags |= UTF_SURROGATE
                end
            else
                throw(UnicodeError(UTF_ERR_NOT_LEAD, pos, ch))
            end
        end
    end
    num3byte != 0 && (flags |= UTF_UNICODE3)
    num4byte != 0 && (flags |= UTF_UNICODE4)
    return totalchar, flags, num4byte, num3byte, num2byte
end

"""
Validates and calculates number of characters in a UTF-8,UTF-16 or UTF-32 encoded vector/string

This function checks the bounds of the start and end positions
Use `unsafe_checkstring` to avoid that overhead if the bounds have already been checked

Input Arguments:

* `dat`    UTF-8 (`Vector{UInt8}`), UTF-16 (`Vector{UInt16}`) or UTF-32 (`Vector{UInt32}`, `AbstractString`) encoded string

Optional Input Arguments:

* `startpos` start position (defaults to 1)
* `endpos`   end position   (defaults to `lastindex(dat)`)

Keyword Arguments:

* `accept_long_null`  = `true`  # Modified UTF-8 (`\\0` represented as `b\"\\xc0\\x80\"`)
* `accept_surrogates` = `true`  # `CESU-8`
* `accept_long_char`  = `false` # Accept arbitrary long encodings

Returns:

* (total characters, flags, 4-byte, 3-byte, 2-byte)

Throws:

* `UnicodeError`
"""
function checkstring end

# No need to check bounds if using defaults
checkstring(dat::AbstractString; kwargs...) = unsafe_checkstring(dat, 1, lastindex(dat); kwargs...)
checkstring(dat; kwargs...) = unsafe_checkstring(dat, 1, length(dat); kwargs...)

# Make sure that beginning and end positions are bounds checked
function checkstring(dat, startpos, endpos = lastindex(dat); kwargs...)
    checkbounds(dat,startpos)
    checkbounds(dat,endpos)
    endpos < startpos && throw(ArgumentError("End position ($endpos) is less than start position ($startpos)"))
    unsafe_checkstring(dat, startpos, endpos; kwargs...)
end

isvalid(str::T) where {T<:Union{ASCIIString,UTF8String,UTF16String,UTF32String}} = isvalid(T, str.data)
isvalid(::Type{T}, str::T) where {T<:Union{ASCIIString,UTF8String,UTF16String,UTF32String}} = isvalid(T, str.data)

byte_string_classify(data::Vector{UInt8}) =
    ccall(:u8_isvalid, Int32, (Ptr{UInt8}, Int), data, length(data))
byte_string_classify(s::ByteString) = byte_string_classify(s.data)
    # 0: neither valid ASCII nor UTF-8
    # 1: valid ASCII
    # 2: valid UTF-8

isvalid(::Type{ASCIIString}, s::Union{Vector{UInt8},ByteString}) = byte_string_classify(s) == 1
isvalid(::Type{UTF8String}, s::Union{Vector{UInt8},ByteString}) = byte_string_classify(s) != 0

bytestring() = ASCIIString("")
function bytestring(s::AbstractString...)
    str = Base.print_to_string(s...)
    data = Vector{UInt8}(codeunits(str))
    isvalid(ASCIIString, data) ? ASCIIString(data) : UTF8String(data)
end
bytestring(s::Vector{UInt8}) = bytestring(String(s))
bytestring(p::Union{Ptr{UInt8},Ptr{Int8},Cstring}) = unsafe_string(p)
bytestring(p::Union{Ptr{UInt8},Ptr{Int8}},len::Integer) = unsafe_string(p,len)


text:
# Neighbors in 1D

function has_neighbor(m::Mesh1D, cell, face)
    m.isperiodic && return true
    if cell.coord == 1 && face == :l
        return false
    elseif cell.coord == length(m.elements)-1 && face == :r
        return false
    else
        return true
    end
end

function neighbor(m::Mesh1D, cell, face)
    @assert has_neighbor(m, cell, face)
    if face == :l
        if m.isperiodic && cell.coord == 1
            return m[length(m.elements)-1]
        else
            return m[cell.coord - 1]
        end
    else
        @assert face == :r
        if m.isperiodic && cell.coord == length(m.elements)-1
            return m[1]
        else
            return m[cell.coord + 1]
        end
    end
end

coord(face) = face == :l ? -1 : 1
oppcoord(face) = -coord(face)

oppface(cell::Cell1D,face) = face == :l ? :r : :l
oppface(m::Mesh1D,cell,face) = oppface(cell,face)

# Neighbors in 2D

has_neighbor(m::Meshes.Mesh, edge::Edge; ns = computeNeighbors(m)) = has_neighbor(m,edge,ns)
has_neighbor(m::Meshes.Mesh, cell::Cell2D, edgeid::Int64; ns = computeNeighbors(m)) = has_neighbor(m,cell,edgeid,ns)
function has_neighbor(m::Meshes.Mesh, edge::Edge, ns)
    ns[edge.cid, edgeid(m,edge)] != -1
end

function has_neighbor(m::Meshes.Mesh, cell::Cell2D, edgeid::Int64, ns)
    ns[cell.cid, edgeid] != -1
end

has_neighbor(m::Meshes.Mesh, cell::Cell2D, edge::Edge, ns) = has_neighbor(m,edge, ns)
has_neighbor(m::Meshes.Mesh, cell::Cell2D, edge::Edge; ns = computeNeighbors(m)) =
    has_neighbor(m,edge, ns)

function neighbor(m::Meshes.Mesh, cell::Cell2D, edgeid; ns = computeNeighbors(m))
    m[ns[cell.cid, edgeid]]
end

neighbor(m::Meshes.Mesh, cell::Cell2D, edge::Edge; ns = computeNeighbors(m)) =
    neighbor(m,cell,edge,ns)

function neighbor(m::Meshes.Mesh, cell::Cell2D, edge::Edge, ns)
    m[ns[cell.cid, edgeid(cell,edge)]]
end

function neighbor(m::Meshes.Mesh, edge::Edge, ns = computeNeighbors(m))
    m[ns[cid(edge), edgeid(m[cid(edge)],edge)]]
end

function oppface(m::Mesh, cell::Cell2D, face, ns)
    n = neighbor(m, cell,face, ns)
    n[edgeid(n,reverse(face))]
end
oppface(m::Mesh, cell::Cell2D, face;  ns = computeNeighbors(m)) = oppface(m,cell,face,ns)

canon(a,b) = (a<b) ? (a,b) : (b,a)

function computeNeighbors(m)
    fs = ((Int64,Int64)=>(Int64,Int64))[]
    ns = Array(Int64,length(m.faces),3)
    fill!(ns,-1)
    for i in 1:length(m.faces)
        f = m.faces[i]
        for (nn,a,b) in ((1,:v1,:v2),(2,:v2,:v3),(3,:v3,:v1))
            c = canon(f.(a),f.(b))
            if haskey(fs,c)
                (j,num) = fs[c]
                ns[i,nn] = j
                ns[j,num] = i
            else
                fs[c] = (i,nn)
            end
        end
    end
    ns
end

# Tests if the vertex p is contained in cell x
function facingface(x::Cell2D,p::Vertex2)
    v0 = x.p2 - x.p1
    v1 = x.p3 - x.p1
    v2 = p - x.p1

    dot00 = dot(v0, v0)
    dot01 = dot(v0, v1)
    dot02 = dot(v0, v2)
    dot11 = dot(v1, v1)
    dot12 = dot(v1, v2)

    invDenom = 1 / (dot00 * dot11 - dot01 * dot01)
    λ1 = u = (dot11 * dot02 - dot01 * dot12) * invDenom
    λ2 = v = (dot00 * dot12 - dot01 * dot02) * invDenom
    λ3 = 1-u-v

    pt = clamp(λ1,0.0,1.0)*x.p1 + clamp(λ2,0.0,1.0)*x.p2 + clamp(λ3,0.0,1.0)*x.p3

    dist = norm(p-pt)

    if (u >= 0) && (v >= 0) && (u + v < 1)
        return (0,0.0)
    elseif (u <= 0) && (v <= 0)
        return ((abs(u) > abs(v) ? 3 : 2),dist)
    elseif (u <= 0)
        return (3,dist)
    elseif (v <= 0)
        return (2,dist)
    else
        return (1,dist)
    end
end

contains(x::Cell2D,p::Vertex2) = facingface(x,p)[1] == 0


text:
import Conda
@static if Sys.iswindows()
    import WinReg
end

include("setup.jl")

const depfile = joinpath(dirname(@__FILE__), "deps.jl")

try
    if isfile(depfile)
        @eval module DepFile; include($depfile); end
    else
        @eval module DepFile; Rhome=libR=""; end
    end

    if !haskey(ENV,"R_HOME") && isdir(DepFile.Rhome) && validate_libR(DepFile.libR)
        Rhome, libR = DepFile.Rhome, DepFile.libR
        @info "Using previously configured R at $Rhome with libR in $libR."
    else
        Rhome = get(ENV, "R_HOME", "")
        if Rhome == "*"
            # install with Conda
            @info "Installing R via Conda.  To use a different R installation,"*
                " set the \"R_HOME\" environment variable and re-run "*
                "Pkg.build(\"RCall\")."
            Conda.add_channel("r")
            Conda.add("r-base")
            Rhome = joinpath(Conda.LIBDIR, "R")

            libR = locate_libR(Rhome)
        else
            if isempty(Rhome)
                try Rhome = readchomp(`R RHOME`); catch; end
            end
            @static if Sys.iswindows()
                if isempty(Rhome)
                    try Rhome = WinReg.querykey(WinReg.HKEY_LOCAL_MACHINE,
                                                "Software\\R-Core\\R", "InstallPath"); catch; end
                end
                if isempty(Rhome)
                    try Rhome = WinReg.querykey(WinReg.HKEY_CURRENT_USER,
                                                "Software\\R-Core\\R", "InstallPath"); catch; end
                end
            else
                if !isdir(Rhome)
                    error("R_HOME is not a directory.")
                end
            end

            isempty(Rhome) && error("R cannot be found. Set the \"R_HOME\" environment variable to re-run Pkg.build(\"RCall\").")
            libR = locate_libR(Rhome)
        end

        @info "Using R at $Rhome and libR at $libR."
        if DepFile.Rhome != Rhome || DepFile.libR != libR
            open(depfile, "w") do f
                println(f, "const Rhome = \"", escape_string(Rhome), '"')
                println(f, "const libR = \"", escape_string(libR), '"')
            end
        end
    end
catch
    if isfile(depfile)
        rm(depfile, force=true) # delete on error
    end
    rethrow()
end


text:
using MultilevelCoordinateSearch
using Base.Test

function camel(x)
    # 6-hump camel function. Typically evaluated over [-3,3] × [-2,2].
    x1, x2 = x[1], x[2]
    x1s = x1*x1
    x2s = x2*x2
    return (4 - 2.1*x1s + x1s*x1s/3)*x1s + x1*x2 + (-4 + 4*x2s)*x2s
end

u = [-3.0,-2.0]
v = -u

@testset "Init: box edges extend to region edges" begin
    xstar, fstar, box = MultilevelCoordinateSearch.init_boxes(camel, [0.0,0.0], ([-3.0,0.0,3.0], [-2.0,0.0,2.0]), u, v)
    @test xstar == [0.0, 0.0]
    @test fstar == camel(xstar)
    @test length(box.parent.children) == 4 && length(box.parent.parent.children) == 4
    r = MultilevelCoordinateSearch.get_root(box)
    @test r.fvalues[2] < min(r.fvalues[1], r.fvalues[3]) && [r.children[i].level for i = 1:4] == [3, 2, 2, 3]
    b = box.parent
    @test b.fvalues[2] < min(b.fvalues[1], b.fvalues[3]) && [b.children[i].level for i = 1:4] == [4, 3, 3, 4]
    @test MultilevelCoordinateSearch.splitindexes(b.children[1]) == (1, 2, 1)
    @test MultilevelCoordinateSearch.splitindexes(b.children[2]) == (1, 2, 2)
    @test MultilevelCoordinateSearch.splitindexes(b.children[3]) == (2, 3, 2)
    @test MultilevelCoordinateSearch.splitindexes(b.children[4]) == (2, 3, 3)
    pqs = MultilevelCoordinateSearch.init_sweep_queues(box, 4)
    @test isempty(pqs[1])
    @test length(pqs[2]) == 1
    @test length(pqs[3]) == 4
    @test length(pqs[4]) == 2
end

@testset "Init: box edges interior to region edges" begin
    xstar, fstar, box = MultilevelCoordinateSearch.init_boxes(camel, [0.0,0.0], ([-2.0,0.0,2.0], [-1.0,0.0,1.0]), u, v)
    @test xstar == [0.0, 0.0]
    @test fstar == camel(xstar)
    @test length(box.parent.children) == 6 && length(box.parent.parent.children) == 6
    r = MultilevelCoordinateSearch.get_root(box)
    @test r.fvalues[2] < min(r.fvalues[1], r.fvalues[3]) && [r.children[i].level for i = 1:6] == [2, 3, 2, 2, 3, 2]
    b = box.parent
    @test b.fvalues[2] == min(b.fvalues[1], b.fvalues[3]) && [b.children[i].level for i = 1:6] == [3, 4, 3, 4, 3, 3]
    @test MultilevelCoordinateSearch.splitindexes(b.children[1]) == (0, 1, 1)
    @test MultilevelCoordinateSearch.splitindexes(b.children[2]) == (1, 2, 1)
    @test MultilevelCoordinateSearch.splitindexes(b.children[3]) == (1, 2, 2)
    @test MultilevelCoordinateSearch.splitindexes(b.children[4]) == (2, 3, 2)
    @test MultilevelCoordinateSearch.splitindexes(b.children[5]) == (2, 3, 3)
    @test MultilevelCoordinateSearch.splitindexes(b.children[6]) == (3, 4, 3)
    pqs = MultilevelCoordinateSearch.init_sweep_queues(box, 4)
    @test isempty(pqs[1])
    @test length(pqs[2]) == 3
    @test length(pqs[3]) == 6
    @test length(pqs[4]) == 2
end

@testset "Init: box edges half to region edges" begin
    xstar, fstar, box = MultilevelCoordinateSearch.init_boxes(camel, [0.0,0.0], ([-3.0,0.0,2.0], [-1.0,0.0,2.0]), u, v)
    @test xstar == [0.0, 0.0]
    @test fstar == camel(xstar)
    @test length(box.parent.children) == 5 && length(box.parent.parent.children) == 5
    r = MultilevelCoordinateSearch.get_root(box)
    @test r.fvalues[2] < min(r.fvalues[1], r.fvalues[3]) && [r.children[i].level for i = 1:5] == [3, 2, 2, 3, 2]
    b = box.parent
    @test (b.fvalues[2] == b.fvalues[1] && b.fvalues[2] < b.fvalues[3]) && [b.children[i].level for i = 1:5] == [3, 4, 3, 3, 4]
    @test MultilevelCoordinateSearch.splitindexes(r.children[1]) == (1, 2, 1)
    @test MultilevelCoordinateSearch.splitindexes(r.children[2]) == (1, 2, 2)
    @test MultilevelCoordinateSearch.splitindexes(r.children[3]) == (2, 3, 2)
    @test MultilevelCoordinateSearch.splitindexes(r.children[4]) == (2, 3, 3)
    @test MultilevelCoordinateSearch.splitindexes(r.children[5]) == (3, 4, 3)
    @test MultilevelCoordinateSearch.splitindexes(b.children[1]) == (0, 1, 1)
    @test MultilevelCoordinateSearch.splitindexes(b.children[2]) == (1, 2, 1)
    @test MultilevelCoordinateSearch.splitindexes(b.children[3]) == (1, 2, 2)
    @test MultilevelCoordinateSearch.splitindexes(b.children[4]) == (2, 3, 2)
    @test MultilevelCoordinateSearch.splitindexes(b.children[5]) == (2, 3, 3)
    pqs = MultilevelCoordinateSearch.init_sweep_queues(box, 4)
    @test isempty(pqs[1])
    @test length(pqs[2]) == 2
    @test length(pqs[3]) == 5
    @test length(pqs[4]) == 2
end


text:
using KissMCMC, BITEModel, VAWTools; const BM=BITEModel
const I2 = BM.ITMIX2
#include("itmix2-runs.jl")

pmcmc = BM.MCMCNum(;pmcmc_defaults...,
                   niter=400,
                   nthin = 1,
                   nburnin=10,
                   nchains=4)

# number of points to fit 1D quantities on
n_1d_theta = 3

good_exps = Dict(gl=>[] for gl in BM.ITMIXv2_glaciers)
bad_gl = Dict(gl=>[] for gl in BM.ITMIXv2_glaciers)
bb = map(x->BM.ITMIXGlacier(x,2), [:Starbuck, :Columbia])
bb = map(x->BM.ITMIXGlacier(x,2), [:Aqqutikitsoq])

for priority=1:4
    println("\n\n++++++++++++++++ PRIORITY $priority ++++++++++++++++++++++++++\n\n")
    for gid in I2.priorities[priority]
        exp = 1 # to make it available
        println("\n\n============================\nRunning Glacier $(BM.getname(gid))")
        #    try
        gl,gb,pp,pm,pn,pl = BM.init_forward(gid; pl_kws...);
        for exp=1:I2.nexps
            println("\n\n---------------\nRunning experiment $exp of $(BM.getname(gid))")
            gln = BM.ITMIX2.make_ITMIX2_glacier(gl, exp)
            gbn = BM.Bands(gb, gl=gln)
            th0d = BM.theta0_dict_defaults(n_1d_theta, pm)
            (theta0, logposterior, logposterior1d, logprior, logprior_debug,
             loglikelihood, loglikelihood1d, fwdm_fn, fwdm1d_fn,
             pmcmc_defaults, fit_target) =
                 BM.init_inverse(gbn, pp, pm, pn, n_1d_theta = n_1d_theta,
                                                               theta0_dict=th0d)
            # some tests:
            fwdm_fn(theta0.th0);
            print("    Time to run forward model:")
            @time hs1d, ivs1d, hs2d, ivs2d, pm_, taus1d, taus1d_l, taus2d =
                fwdm_fn(theta0.th0);
            print("    Posterior value:")
            @show logposterior(theta0.th0)[1]
            print("    Time to calculate posterior:")
            @time logposterior(theta0.th0)[1]

            # sample posterior
            thetas, theta_expect, accept_ratio, blobs,
            Rhat, sample_size, nthin, reshape_revert,
            mc_hs1d, mc_ivs1d,
            mh2d, sh2d, miv2d, siv2d, min_h2d, max_h2d, min_iv2d, max_iv2d = mcmc(logposterior, theta0, pmcmc; verbose=true);
            print_results(thetas, accept_ratio, names=BM.get_varnames(theta0))
            push!(good_exps[gid], exp)
            I2.save_itmix2_run(mh2d, gln, pl)
        end
        # catch e
            #     println("\nFailed with $e\n\n")
        #     push!(bad_gl[gid], (exp, e))
        # end
    end
end


text:
module MSTM

include("constants.jl")
include("special_functions.jl")
include("data.jl")
include("mie.jl")
include("wrapper.jl")

end


text:
using Distributions

function random_spatial_network(N::Integer; # number of nodes
    d::Integer = 2, # dimension
    Pk = Poisson(10), # Expected degree distribution
    Pf = Normal(0.0, 0.03) # Distance kernel
    )
    ks = rand(Pk, N)
    xs = rand(N,d)
    dists = Array{Float64,2}(undef, N, N)
    for i in 1:N, j in i:N
        dists[i,j] = norm(xs[i,:]-xs[j,:])
        dists[j,i] = dists[i,j]
    end
    f = x -> pdf(Pf, x)
    return random_spatial_network(dists, ks, f, 1), xs
end

function random_spatial_network(distances::AbstractMatrix, ks, f, ρ)
    N = size(distances, 1)
    k = mean(ks)
    ps = (ks * ks') .* f.(distances) / (ρ * k)
    return edge_probability_graph(ps)
end

function edge_probability_graph(ps::AbstractMatrix)
    N = size(ps, 1)
    edges = [(i,j) for i in 1:N-1, j in (i+1):N if rand() < ps[i,j]]
    return SimpleGraphFromIterator(edges)
end

max_stars_repo_name: csimal/metaplex-networks

max_issues_repo_name: csimal/metaplex-networks

max_forks_repo_name: csimal/metaplex-networks


text:
export ex1, fex1, fex2, bayes_drift, test1, levelK, finger_pm, visualize_posterior

#%  Introduction
#%  ~~~~~~~~~~~~
#%
#%
#%  The procedure is as follows.
#%  Consider the diffusion process :math:`(x_t\colon 0 \le t \le T)` given by
#%
#%  	:math:`dx_t = b(x_t) dt + dw_t`
#%
#%
#%  where the drift ``b`` is expressed as linear combination
#%
#%  	:math:`f(x) = \sum_{i =1}^n c_i \phi_i(x)`
#%
#%  (see :ref:`modschauder`) and
#%  prior distribution on the coefficients
#%
#%  	:math:`c_i \sim N(0,\xi_i)`
#%
#%  Then the posterior distribution of :math:`b` given observations :math:`x_t` is given by
#%
#%  	:math:`c_i | x_s \sim N(W^{-1}\mu, W^{-1})`
#%  	:math:`W = \Sigma + (\operatorname{diag}(\xi))^{-1}`,
#%
#%  with the nxn-matrix
#%
#%  	:math:`\Sigma_{ij} = \int_0^T \phi_i(x_t)\phi_j(x_t) dt`
#%
#%  and the n-vector
#%
#%  	:math:`\mu_i = \int_0^T \phi_i(x_t) d x_t`.
#%
#%
#%  Using the recursion detailed in :ref:`modschauder`, one rather computes
#%
#%  	:math:`\Sigma^\prime_{ij} = \int_0^T \psi_i(x_t)\psi_j(x_t) dt`
#%
#%  and the n-vector
#%
#%  	:math:`\mu^\prime_i = \int_0^T \psi_i(x_t) d x_t`
#%
#%  and uses ``pickup_mu!(mu)`` and ``pickup_Sigma!(Sigma)`` to obtain :math:`\mu` and :math:`\Sigma`.
#%

#%
#%  Optional additional basis functions
#%  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#%
#%  One can extend the basis by additional functions, implemented are variants. ``B1`` includes a constant, ``B2`` two linear functions
#%
#%  ``B1``
#%  	:math:`\phi_1 \dots \phi_n, c`
#%
#%  ``B2``
#%  	:math:`\phi_1 \dots \phi_n, \max(1-x, 0), \max(x, 0)`
#%
#%  To compute ``mu``, use
#%
#%  .. code-block:: matlab
#%
#%  	mu = pickup_mu!(fe_mu(y,L, 0))
#%  	mu = fe_muB1(mu, y);
#%
#%  or
#%
#%  .. code-block:: matlab
#%
#%  	mu = pickup_mu!(fe_mu(y,L, 0))
#%  	mu = fe_muB2(mu, y);
#%
#%

#%  Reference
#%  ~~~~~~~~~

#%
#%  Functions taking ``y` without parameter [a,b] expect ``y`` to be shifted into the intervall ``[0,1]``.
#%

#%  .. function:: pickup_mu!(mu)
#%
#%  	computes mu from mu'
#%

function pickup_mu!(mu)
	L = level(mu)
	n = 2^L-1
	for l in 1:L-1
		for i in 2^l:2^l:n
	 		 mu[i] *= 2
			 mu[i] += mu[i-2^(l-1)] + mu[i+2^(l-1)]
		end

	end
	mu
end

#%  .. function:: drop_mu!(mu)
#%
#%  	Computes mu' from mu.
#%

function drop_mu!(mu)
	L = level(mu)
	n = 2^L-1
	for l in L-1:-1:1
		for i in 2^l:2^l:n

			mu[i] -= mu[i-2^(l-1)] + mu[i+2^(l-1)]
			mu[i] /= 2
		end

	end
	mu
end

#%  .. function:: pickup_Sigma!(Sigma)
#%
#%  	Transforms Sigma' into Sigma.
#%

function pickup_Sigma!(Sigma)
	L = level(Sigma)
	n = 2^L-1
	for l in 1:L-1
		for i in 2^l:2^l:n
		 Sigma[i, :] *= 2
		 Sigma[:, i] *= 2
		 Sigma[i, :] += Sigma[i-2^(l-1), :] + Sigma[i+2^(l-1), :]
		 Sigma[:, i] += Sigma[:, i-2^(l-1)] + Sigma[:, i+2^(l-1)]
		end

	end
	Sigma
end

#%  .. function:: drop_Sigma!(Sigma)
#%
#%  	Transforms Sigma into Sigma'.
#%

function drop_Sigma!(Sigma)
	L = level(Sigma)
	n = 2^L-1

	for l in L-1:-1:1
		for i in 2^l:2^l:n
		 Sigma[i, :] -= Sigma[i-2^(l-1), :] + Sigma[i+2^(l-1), :]
		 Sigma[:, i] -= Sigma[:, i-2^(l-1)] + Sigma[:, i+2^(l-1)]
		 Sigma[i, :] /= 2
		 Sigma[:, i] /= 2
		end

	end
	Sigma
end

#use after drop_Sigma
function drop_SigmaB2!(Sigma)
	(L,K) = levelK(Sigma)
	assert(K==2)
	n = 2^L-1
	for i in 1:n
		 Sigma[i, :] -= Sigma[end, :]*i/(n+1) + Sigma[end-1, :]*(n+1-i)/(n+1)
		 Sigma[:, i] -= Sigma[:, end]*i/(n+1)  + Sigma[:, end-1]*(n+1-i)/(n+1)

	end
	Sigma
end

#use after drop_Sigma
function drop_SigmaB1!(Sigma)
	(L,K) = levelK(Sigma)
	assert(K==1)
	n = 2^L-1
	for i in 1:n
		 Sigma[i, :] -= Sigma[end, :]
		 Sigma[:, i] -= Sigma[:, end]

	end
	Sigma
end

compose(f, g) = x -> f(g(x))

drop_mu = compose(drop_mu!,copy)
drop_Sigma = compose(drop_Sigma!,copy)

#%  .. function:: fe_mu(y, L, K)
#%
#%  	Computes mu' from the observations `y` using ``2^L-1`` basis elements
#%  	and returns a vector with ``K`` trailing zeros (in case one ones to customize
#%  	the basis.

function fe_mu_jl(y, L, K)
	n = 2^(L-1) #number of even element/elements in lowest level!
	mu = zeros(2n-1 + K)

	dy = [diff(y),0]
	for i in 1:n - 1
		mu[2i-1] = dot(hat(y*n .- (i - 1)),dy)
		mu[2i] = dot(hat(y*n .- (i -.5)),dy)
	end
	mu[2n-1] = dot(hat(y*n.- (n-1)),dy)
	mu
end

#%  .. function:: fe_muB1(mu, y)
#%
#%  	Append :math:`\mu_{n+1} = \int_0^T \phi_{n+1} d x_t` with :math:`\phi_{n+1} = 1`.
#%

function fe_muB1(mu, y)
	[mu, y[end] - y[1]]
end

#%  .. function:: fe_muB2(mu, y)
#%
#%  	Append :math:`\mu_{n+1} = \int_0^T \phi_{n+1} d x_t` with :math:`\phi_{n+1} =  \max(1-x, 0)`
#%  	and :math:`\mu_{n+2} = \int_0^T \phi_{n+2} d x_t` with :math:`\phi_{n+2} =  \max(x, 0)`
#%

function fe_muB2(mu,y)
	dy = [diff(y),0]
	[mu, dot(max(1 .- y, 0), dy), dot(max(y, 0), dy)]
end

#%  .. function:: fe_Sigma(y, dt, L)
#%
#%  	Computes the matrix Sigma' from the observations `y` uniformly spaced at distance ``dt``
#%	using ``2^L-1`` basis elements.
#%

#int(chol(finger_matrix(pickup_Sigma!(fe_Sigma(rand(10000), 0, 1, 5)))) .!= 0)
#memory intensive, using dot
#dont be confused, even wavelets correspond to odd indices as julia indices start at 1
function fe_Sigma_dot(y, dt::Float64, L)
	n = 2^(L-1)
	S = zeros(2n-1, 2n-1)
	yn = y*n
	evnext = hat(yn)
	for i in 1 : n-1
		ev = evnext
		od = hat(yn .- (i - 0.5))
		evnext = hat(yn .-i)
		S[2i-1, 2i-1] = dot(ev,ev)*dt
		S[2i, 2i] = dot(od,od)*dt
		S[2i-1, 2i] = S[2i, 2i-1] = dot(od, ev)*dt
		S[2i, 2i+1] = S[2i+1, 2i] = dot(evnext,od)*dt
	end
	S[2n-1, 2n-1] = dot(evnext, evnext)*dt
	S
end

#wavelets with a constant
function fe_SigmaB1_dot(y, dt::Float64, L)
	K = 1
	n = 2^(L-1)
	S = zeros(2n-1 + K, 2n-1+K)
	yn = y*n
	evnext = hat(yn)
	for i in 1:n - 1
		ev = evnext
		od = hat(yn .- (i - 0.5))
		evnext = hat(yn .- i)

		S[2i-1, 2i-1] = dot(ev,ev)*dt
		S[2i, 2i] = dot(od,od)*dt
		S[2i-1, 2i] = S[2i, 2i-1] = dot(od, ev)*dt
		S[2i, 2i+1] = S[2i+1, 2i] = dot(evnext,od)*dt
		S[end, 2i-1] = S[2i-1, end] = sum(ev)*dt
		S[2i, end] = S[end, 2i] = sum(od)*dt
	end
	S[2n-1, 2n-1] = dot(evnext, evnext)*dt
	S[end, 2n-1] = S[2n-1, end] = sum(evnext)*dt
	S[end, end]  = length(y)*dt
	S
end

function fe_SigmaB2_dot(y, dt::Float64, L)
	K = 2
	n = 2^(L-1)
	S = zeros(2n-1 + K, 2n-1+K)
	yn = y*n
	evnext = hat(yn)
	k1 = max(1-y, 0)
	k2 = max(y, 0)

	for i in 1:n - 1
		ev = evnext
		od = hat(yn .- (i - 0.5))
		evnext = hat(yn .- i)

		S[2i-1, 2i-1] = dot(ev,ev)*dt
		S[2i, 2i] = dot(od,od)*dt
		S[2i-1, 2i] = S[2i, 2i-1] = dot(od, ev)*dt
		S[2i, 2i+1] = S[2i+1, 2i] = dot(evnext,od)*dt
		S[end, 2i-1] = S[2i-1, end] = dot(ev,k2)*dt
		S[2i, end] = S[end, 2i] = dot(od, k2)*dt
		S[end-1, 2i-1] = S[2i-1, end-1] = dot(ev,k1)*dt
		S[2i, end-1] = S[end-1, 2i] = dot(od, k1)*dt

	end
	S[2n-1, 2n-1] = dot(evnext, evnext)*dt
	S[end, 2n-1] = S[2n-1, end] = dot(k2,evnext)*dt
	S[end, end]  = dot(k2,k2)*dt
	S[end-1, 2n-1] = S[2n-1, end-1] = dot(k1,evnext)*dt
	S[end-1, end-1] = dot(k1,k1)*dt
	S[end, end-1] = S[end-1, end] = dot(k1,k2)*dt
	S
end

#variant accessing y once, fastest if implemented in c

function fe_Sigma_at(y, dt::Float64, L)
	n = 2^(L-1)
	N = length(y)

	S = zeros(2n-1, 2n-1)

	for t = 1:N
		yn = y[t]*n
		i = clamp(ceil(yn), 1, n)
		j = clamp(ceil(yn .- 0.5), 1, n-1)
		S[2i-1, 2i-1] += ((hat(yn .- (i - 1)).^2))*dt
		S[2j, 2j] += ((hat(yn .- (j - 0.5)).^2))*dt
		S[2i-1, 2j] = S[2j, 2i-1] += (hat(yn .- (i - 1)).*hat(yn .- (j - 0.5)))*dt
	end
	S
end

fe_mu = fe_mu_jl
fe_Sigma = fe_Sigma_dot
fe_SigmaB1 = fe_SigmaB1_dot
fe_SigmaB2 = fe_SigmaB2_dot

#can libsigma be used?
if find_library(["libsigma"], [Pkg.dir("SDE","deps")]) != ""
	include("sigma.jl")

end

# fe_Sigma([0, 0.25, 0.5, 0.625,0], 2, 0,1, 2)
#3x3 Float64 Array:
# 2.0  0.0  0.0
# 0.0  2.5  0.5
# 0.0  0.5  0.5

#expect [xiA, ..., xiB, xi1, ..., xiL] where xiA ... xiB are variances for the trailing elements
#and xi is the variance on level i with 2^(i-1) elements
function expand_levelwise_coeff(xi,K)
	L = length(xi)-K
	rev = 2^L
	xifull = zeros(2^L-1 + K)
	for i in 1:L
		xifull[rev-(2^i-1):rev-2^(i-1)] = xi[K+i]
	end
	xifull[end-K+1:end] = xi[1: K]
	p = finger_pm(L,K)
	pm = permutationmatrix(p)

	ipermute!(xifull, p)
end

#%  .. function:: bayes_drift(x, dt, a, b, L, xirem, beta, B)
#%
#%  	Performs estimation of drift on observations ``x`` in [a,b] spaced at distance ``dt``
#%	using the Schauder basis of level ``L`` and level wise coefficients decaying at rate ``beta``.
#%	A Brownian motion like prior is obtained for beta= 0.5. The ``K`` remaining optional
#%  	basiselements have variance ``xirem``.
#%
#%  	The result is returned as ``[designp coeff se]`` where ``coeff`` are coefficients of finite elements with maximum at the designpoints ``designp`` and standard error ``se``.
#%
#%
#%
#%	Observations outside [a,b] may influence the result through ``phi_{n+1}, ..., phi_{n+K}``

function bayes_drift(x, dt, a, b, L, xirem, beta, B)

	n = 2^L -1

	# shift process
	y = (x-a)/(b-a)

	# number of additional elements
	K = 0

	if (B == "B1"); K = 1;
	elseif (B == "B2"); K = 2;
	end

	xil = [xirem, 1.*2.^(-2.*[L-1:-1:0]    -beta*[1:L])]

	println("preparations"); tic()
	xi = expand_levelwise_coeff(xil, K)
 	toq()
	println("Xi ", xi')

	println("mu"); tic()
#	mu = pickedup_mu(y, L)
	mu = pickup_mu!(fe_mu(y,L, 0))
	toq()

	if(B=="B1"); mu = fe_muB1(mu, y);
	elseif(B=="B2"); mu = fe_muB2(mu, y);
	end

	println("mu ", mu)

	println("Sigma pickup"); tic()

	if B=="B1"
		feS = fe_SigmaB1(y,dt, L)
	elseif B == "B2"
		feS = fe_SigmaB2(y,dt, L)
	else
		feS = fe_Sigma(y,dt, L)
	end

	println("Sigma f.e.", feS)
	toc()
	tic()
	sigma = pickup_Sigma!(feS)
	toc()
	println("Sigma schauder ", sigma)

	println("chol, coeff")
	tic()
	M = chol(sigma + diagm(xi.^(-1)))
	toc()

	tic()
	coeff = M\(M'\mu)
	toc()

	ci2 =  inv(M)*inv(M')

 	drop_Sigma!(ci2) #check
	drop!(coeff)
	coeff .*= (b-a)

	if B=="B1"
		drop_SigmaB1!(ci2)
		se = diag(ci2).^(.5)
		coeff = dropB1(coeff)
		designp = [(1:n+1)*(b-a)/(n+1)+a]
	elseif B=="B2"
		drop_SigmaB2!(ci2)
		se = diag(ci2).^(.5)
		se = se[[n+1, 1:n, n+2]]
		coeff = dropB2(coeff)
		designp = [(0:n+1)*(b-a)/(n+1)+a]
	else
		se = diag(ci2).^(.5)
		se = 0.5(se + [se[2],se[1:end-1]])
		designp = [(1:n+K)*(b-a)/(n+1)+a]
	end

	println("coeff f.e.", coeff)

	[designp coeff se*(b-a)]
end


text:
@testset "Amplifiers" begin
  @testset "SimulatedAmplifier" begin
    include("SimulatedAmplifierTest.jl")
  end

  if "hubert" in ARGS
    @testset "HubertAmplifier" begin
      include("HubertAmplifierTest.jl")
    end
  end
end


text:
println("➡ Running Tests for Soliton problem ...")
include("Sech2_Soliton.jl")

#println("----------------------------------")

#println("➡ Running Tests for SPM problem ...")
#include("SPM.jl")


text:
# # Quasi-Linear forced-dissipative barotropic QG beta-plane turbulence
#
#md # This example can be run online via [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/barotropicqgql_betaforced.ipynb).
#md # Also, it can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/barotropicqgql_betaforced.ipynb).
#
# A simulation of forced-dissipative barotropic quasi-geostrophic turbulence on
# a beta plane under the *quasi-linear approximation*. The dynamics include
# linear drag and stochastic excitation.

using FourierFlows, Plots, Statistics, Printf, Random

using FourierFlows: parsevalsum
using FFTW: irfft
using Random: seed!
using Statistics: mean

import GeophysicalFlows.BarotropicQGQL
import GeophysicalFlows.BarotropicQGQL: energy, enstrophy

# ## Choosing a device: CPU or GPU

dev = CPU()     # Device (CPU/GPU)
nothing # hide

# ## Numerical parameters and time-stepping parameters

      n = 128            # 2D resolution = n^2
stepper = "FilteredRK4"  # timestepper
     dt = 0.05           # timestep
 nsteps = 8000           # total number of time-steps
 nsubs  = 10             # number of time-steps for intermediate logging/plotting (nsteps must be multiple of nsubs)
nothing # hide

# ## Physical parameters

L = 2π        # domain size
β = 10.0      # planetary PV gradient
μ = 0.01      # bottom drag
nothing # hide

# ## Forcing

# We force the vorticity equation with stochastic excitation that is delta-correlated
# in time and while spatially homogeneously and isotropically correlated. The forcing
# has a spectrum with power in a ring in wavenumber space of radious $k_f$ and
# width $\delta k_f$, and it injects energy per unit area and per unit time equal
# to $\varepsilon$.

forcing_wavenumber = 14.0    # the central forcing wavenumber for a spectrum that is a ring in wavenumber space
forcing_bandwidth  = 1.5     # the width of the forcing spectrum
ε = 0.001                    # energy input rate by the forcing

gr  = TwoDGrid(n, L)

K = @. sqrt(gr.Krsq)                      # a 2D array with the total wavenumber
k = [gr.kr[i] for i=1:gr.nkr, j=1:gr.nl]  # a 2D array with the zonal wavenumber

forcing_spectrum = @. exp(-(K - forcing_wavenumber)^2 / (2 * forcing_bandwidth^2))
forcing_spectrum[K .<  2 * 2π/L] .= 0    # no power at low wavenumbers
forcing_spectrum[K .> 20 * 2π/L] .= 0    # no power at high wavenumbers
forcing_spectrum[k .< 2π/L ] .= 0        # make sure forcing does not have power at k=0
ε0 = parsevalsum(forcing_spectrum .* gr.invKrsq / 2, gr) / (gr.Lx * gr.Ly)
@. forcing_spectrum *= ε/ε0               # normalize forcing to inject energy at rate ε

seed!(1234) # reset of the random number generator for reproducibility
nothing # hide

# Next we construct function `calcF!` that computes a forcing realization every timestep
function calcF!(Fh, sol, t, clock, vars, params, grid)
  ξ = ArrayType(dev)(exp.(2π * im * rand(eltype(grid), size(sol))) / sqrt(clock.dt))
  @. Fh = ξ * sqrt.(forcing_spectrum)
  Fh[abs.(grid.Krsq) .== 0] .= 0
  return nothing
end
nothing # hide

# ## Problem setup
# We initialize a `Problem` by providing a set of keyword arguments. Not providing
# a viscosity coefficient ν leads to the module's default value: ν=0. In this
# example numerical instability due to accumulation of enstrophy in high wavenumbers
# is taken care with the `FilteredTimestepper` we picked.
prob = BarotropicQGQL.Problem(dev; nx=n, Lx=L, β=β, μ=μ, dt=dt, stepper=stepper,
                              calcF=calcF!, stochastic=true)
nothing # hide

# and define some shortcuts.
sol, cl, vs, pr, gr = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
x, y = gr.x, gr.y
nothing # hide

# First let's see how a forcing realization looks like.
calcF!(vs.Fh, sol, 0.0, cl, vs, pr, gr)

heatmap(x, y, irfft(vs.Fh, gr.nx),
     aspectratio = 1,
               c = :balance,
            clim = (-8, 8),
           xlims = (-gr.Lx/2, gr.Lx/2),
           ylims = (-gr.Ly/2, gr.Ly/2),
          xticks = -3:3,
          yticks = -3:3,
          xlabel = "x",
          ylabel = "y",
           title = "a forcing realization",
      framestyle = :box)

# ## Setting initial conditions

# Our initial condition is simply fluid at rest.
BarotropicQGQL.set_zeta!(prob, zeros(gr.nx, gr.ny))
nothing # hide

# ## Diagnostics

# Create Diagnostics -- `energy` and `enstrophy` are functions imported at the top.
E = Diagnostic(energy, prob; nsteps=nsteps)
Z = Diagnostic(enstrophy, prob; nsteps=nsteps)
nothing # hide

# We can also define our custom diagnostics via functions.
function zetaMean(prob)
  sol = prob.sol
  sol[1, :]
end

zMean = Diagnostic(zetaMean, prob; nsteps=nsteps, freq=10)  # the zonal-mean vorticity
nothing # hide

# We combile all diags in a list.
diags = [E, Z, zMean] # A list of Diagnostics types passed to "stepforward!" will  be updated every timestep.
nothing # hide

# ## Output

# We choose folder for outputing `.jld2` files and snapshots (`.png` files).
filepath = "."
plotpath = "./plots_forcedbetaQLturb"
plotname = "snapshots"
filename = joinpath(filepath, "forcedbetaQLturb.jld2")
nothing # hide

# Do some basic file management,
if isfile(filename); rm(filename); end
if !isdir(plotpath); mkdir(plotpath); end
nothing # hide

# and then create Output.
get_sol(prob) = sol # extracts the Fourier-transformed solution
get_u(prob) = irfft(im*g.l.*g.invKrsq.*sol, g.nx)
out = Output(prob, filename, (:sol, get_sol), (:u, get_u))
nothing # hide

# ## Visualizing the simulation

# We define a function that plots the vorticity and streamfunction fields, the
# corresponding zonal-mean vorticity and zonal-mean zonal velocity and timeseries
# of energy and enstrophy.

function plot_output(prob)
  ζ̄, ζ′= prob.vars.Zeta, prob.vars.zeta
  ζ = @. ζ̄ + ζ′
  ψ̄, ψ′= prob.vars.Psi,  prob.vars.psi
  ψ = @. ψ̄ + ψ′
  ζ̄ₘ = mean(ζ̄, dims=1)'
  ūₘ = mean(prob.vars.U, dims=1)'

  pζ = heatmap(x, y, ζ,
       aspectratio = 1,
            legend = false,
                 c = :balance,
              clim = (-8, 8),
             xlims = (-gr.Lx/2, gr.Lx/2),
             ylims = (-gr.Ly/2, gr.Ly/2),
            xticks = -3:3,
            yticks = -3:3,
            xlabel = "x",
            ylabel = "y",
             title = "vorticity ζ=∂v/∂x-∂u/∂y",
        framestyle = :box)

  pψ = contourf(x, y, ψ,
            levels = -0.32:0.04:0.32,
       aspectratio = 1,
         linewidth = 1,
            legend = false,
              clim = (-0.22, 0.22),
                 c = :viridis,
             xlims = (-gr.Lx/2, gr.Lx/2),
             ylims = (-gr.Ly/2, gr.Ly/2),
            xticks = -3:3,
            yticks = -3:3,
            xlabel = "x",
            ylabel = "y",
             title = "streamfunction ψ",
        framestyle = :box)

  pζm = plot(ζ̄ₘ, y,
            legend = false,
         linewidth = 2,
             alpha = 0.7,
            yticks = -3:3,
             xlims = (-3, 3),
            xlabel = "zonal mean ζ",
            ylabel = "y")
  plot!(pζm, 0*y, y, linestyle=:dash, linecolor=:black)

  pum = plot(ūₘ, y,
            legend = false,
         linewidth = 2,
             alpha = 0.7,
            yticks = -3:3,
             xlims = (-0.5, 0.5),
            xlabel = "zonal mean u",
            ylabel = "y")
  plot!(pum, 0*y, y, linestyle=:dash, linecolor=:black)

  pE = plot(1,
             label = "energy",
         linewidth = 2,
             alpha = 0.7,
             xlims = (-0.1, 4.1),
             ylims = (0, 0.05),
            xlabel = "μt")

  pZ = plot(1,
             label = "enstrophy",
         linecolor = :red,
            legend = :bottomright,
         linewidth = 2,
             alpha = 0.7,
             xlims = (-0.1, 4.1),
             ylims = (0, 5),
            xlabel = "μt")

  l = @layout grid(2, 3)
  p = plot(pζ, pζm, pE, pψ, pum, pZ, layout=l, size = (1000, 600), dpi=150)

  return p
end

# ## Time-stepping the `Problem` forward

# We time-step the `Problem` forward in time.

p = plot_output(prob)

startwalltime = time()

anim = @animate for j=0:Int(nsteps/nsubs)

  cfl = cl.dt*maximum([maximum(vs.v)/gr.dy, maximum(vs.u+vs.U)/gr.dx])

  log = @sprintf("step: %04d, t: %d, cfl: %.2f, E: %.4f, Q: %.4f, walltime: %.2f min",
    cl.step, cl.t, cfl, E.data[E.i], Z.data[Z.i],
    (time()-startwalltime)/60)

  if j%(1000/nsubs)==0; println(log) end

  p[1][1][:z] = @. vs.zeta + vs.Zeta
  p[1][:title] = "vorticity, μt="*@sprintf("%.2f", μ*cl.t)
  p[4][1][:z] = @. vs.psi + vs.Psi
  p[2][1][:x] = mean(vs.Zeta, dims=1)'
  p[5][1][:x] = mean(vs.U, dims=1)'
  push!(p[3][1], μ*E.t[E.i], E.data[E.i])
  push!(p[6][1], μ*Z.t[Z.i], Z.data[Z.i])

  stepforward!(prob, diags, nsubs)
  BarotropicQGQL.updatevars!(prob)

end

mp4(anim, "barotropicqgql_betaforced.mp4", fps=18)

# ## Save

# Finally save the last snapshot.
savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), cl.step)
savefig(savename)


text:
@testset "macro methods" begin
    @testset "tensorexpressions" begin
        using TensorOperations: isassignment, isdefinition, getlhsrhs, isindex, makeindex,
            istensor, maketensor, isgeneraltensor, makegeneraltensor, istensorexpr, isscalarexpr,
            makescalar, hastraceindices, getindices, getallindices, geteltype

        @test isassignment(:(a[-1,-2,-3] = b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))
        @test isassignment(:(a[-1,-2,-3] += b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))
        @test isassignment(:(a[-1,-2,-3] -= b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))
        @test !isassignment(:(a[-1,-2,-3] := b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))
        @test !isassignment(:(a[-1,-2,-3] ≔ b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))
        @test !isassignment(:(b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))

        @test !isdefinition(:(a[-1,-2,-3] = b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))
        @test !isdefinition(:(a[-1,-2,-3] += b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))
        @test !isdefinition(:(a[-1,-2,-3] -= b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))
        @test isdefinition(:(a[-1,-2,-3] := b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))
        @test isdefinition(:(a[-1,-2,-3] ≔ b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))
        @test !isdefinition(:(b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))

        lhs = :(a[-1,-2,-3])
        rhs = :(b[-1,-2,1]*c[1,-3]+d[-2,-3,-1])
        @test getlhsrhs(:(a[-1,-2,-3] = b[-1,-2,1]*c[1,-3]+d[-2,-3,-1])) == (lhs, rhs)
        @test getlhsrhs(:(a[-1,-2,-3] += b[-1,-2,1]*c[1,-3]+d[-2,-3,-1])) == (lhs, rhs)
        @test getlhsrhs(:(a[-1,-2,-3] -= b[-1,-2,1]*c[1,-3]+d[-2,-3,-1])) == (lhs, rhs)
        @test getlhsrhs(:(a[-1,-2,-3] := b[-1,-2,1]*c[1,-3]+d[-2,-3,-1])) == (lhs, rhs)
        @test getlhsrhs(:(a[-1,-2,-3] ≔ b[-1,-2,1]*c[1,-3]+d[-2,-3,-1])) == (lhs, rhs)

        @test isindex(:a)
        @test makeindex(:a) == :a
        @test isindex(:(a'))
        @test makeindex(:(a')) == :(a′)
        @test isindex(:(a'''))
        @test makeindex(:(a''')) == :(a′′′)
        @test isindex(:(β))
        @test makeindex(:(β)) == :(β)
        @test isindex(:(β'))
        @test makeindex(:(β')) == :(β′)
        @test isindex(:(3))
        @test makeindex(:(3)) == 3
        @test isindex(:(-5))
        @test makeindex(:(-5)) == -5
        @test !isindex(:('a'))
        @test !isindex(:(a+b))
        @test !isindex(:("x"))
        @test !isindex(:(5.1))

        @test istensor(:(a[1,2,3]))
        @test maketensor(:(a[1,2,3])) == (esc(:a), Any[1,2,3], Any[])
        @test istensor(:(a[5][a b c]))
        @test maketensor(:(a[5][a b c])) == (esc(:(a[5])), Any[:a,:b,:c], Any[])
        @test istensor(:(cos(y)[a b c; 1 2 3]))
        @test maketensor(:(cos(y)[a b c; 1 2 3])) == (esc(:(cos(y))), Any[:a,:b,:c], Any[1,2,3])
        @test istensor(:(x[_; 1 2 3]))
        @test maketensor(:(x[_; 1 2 3])) == (esc(:x), Any[], Any[1,2,3])
        @test !istensor(:(2*a[1,2,3]))
        @test !istensor(:(a[1 2 3; 4 5 6; 7 8 9]))
        @test !istensor(:(conj(a[5][a b c])))
        @test !istensor(:(cos(y)[a b c; 1 2 3]*b[4,5]))
        @test !istensor(:(3+5))

        @test isgeneraltensor(:(conj(a[1,2,3])))
        @test makegeneraltensor(:(conj(a[1,2,3]))) == (esc(:a), [1,2,3], [], makescalar(:(conj(1))), true)
        @test isgeneraltensor(:(x*a[5][a b c]))
        @test makegeneraltensor(:(x*a[5][a b c])) == (esc(:(a[5])), [:a,:b,:c], [], makescalar(:(x*1)), false)
        @test isgeneraltensor(:(3*conj(a*cos(y)[a b c; 1 2 3])))
        @test makegeneraltensor(:(3*conj(a*cos(y)[a b c; 1 2 3]))) == (esc(:(cos(y))),  Any[:a,:b,:c],  Any[1,2,3], makescalar(:(3*conj(a*1))), true)
        @test !isgeneraltensor(:(1/a[1,2,3]))
        @test !isgeneraltensor(:(a[1 2 3; 4 5 6]\x))
        @test !isgeneraltensor(:(cos(y)[a b c; 1 2 3]*b[4,5]))
        @test !isgeneraltensor(:(3+5))

        @test hastraceindices(:(a[x,y,z,1,x]))
        @test hastraceindices(:(a[x y z;1 x]))
        @test !hastraceindices(:(a[-1,y,z,1,x]))
        @test !hastraceindices(:(a[-1 y z;1 x]))

        @test isscalarexpr(:a)
        @test isscalarexpr(2.)
        @test isscalarexpr(:(2. + im*3))
        @test isscalarexpr(:(3*a+c))
        @test isscalarexpr(:(sin(x)+exp(-y)))
        @test isscalarexpr(:(scalar(a[x,y]*b[y,x])))
        @test !isscalarexpr(:(a[x,y]*b[y,x]))
        @test !isscalarexpr(:(3*scalar(a[x,y]*b[y,x]) + conj(c[z])))

        @test geteltype(:(a[1,2,3]*b[3,4,5]+c[1,2,4,5])) == :(promote_type(promote_type(eltype($(Expr(:escape, :a))), eltype($(Expr(:escape, :b)))), eltype($(Expr(:escape, :c)))))
    end

    @testset "parsecost" begin
        using TensorOperations: parsecost, Power
        @test parsecost(:(3/5)) === 3/5
        @test parsecost(:(5+2)) === 5+2
        @test parsecost(:(Int128(2*8))) === Int128(2*8)
        @test parsecost(:(float(2*8))) === float(2*8)
        @test parsecost(:(big(2*8))) isa BigInt
        @test parsecost(:(x)) == Power{:x}(1,1)
        @test parsecost(:(3*x^3)) == Power{:x}(3,3)
        @test parsecost(:(3*x^3+2)) == Power{:x}(3,3)+2
        @test parsecost(:((3*x^3-2)/5)) == (Power{:x}(3,3)-2)/5
    end

end

@testset "methods for tensoropt" begin
    @testset "set methods" for T in (UInt32,UInt64,UInt128,BitSet,BitVector)
        maxint = 32
        A = Set(rand(1:maxint, 20))
        B = Set(rand(1:maxint, 20))
        U = union(A,B)
        I = intersect(A,B)
        S = setdiff(A,B)

        A2 = TensorOperations.storeset(T, A, maxint)
        B2 = TensorOperations.storeset(T, B, maxint)
        U2 = TensorOperations.storeset(T, U, maxint)
        I2 = TensorOperations.storeset(T, I, maxint)
        S2 = TensorOperations.storeset(T, S, maxint)

        @test TensorOperations._union(A2, B2) == U2
        @test TensorOperations._intersect(A2, B2) == I2
        @test TensorOperations._setdiff(A2, B2) == S2
        @test TensorOperations._isemptyset(TensorOperations.storeset(T, [], maxint))
        @test !TensorOperations._isemptyset(TensorOperations.storeset(T, [1], maxint))
    end

    @testset "poly" begin
        using TensorOperations: Power, degree, Poly
        x = Power{:x}(1,1)
        x⁰ = Power{:x}(1,0)
        x⁵= Power{:x}(1,5)
        @test x^5 == x⁵ == x*x*x*x*x
        @test x^0 == x⁰ == one(x)
        @test zero(x) == 0*x
        @test (5*x)/5 == x
        @test degree(x) == 1
        @test x + x^3 + x^5 isa Poly
        p = x - x^3 + x^5
        coeffs = [0,1,0,-1,0,1]
        @test p == Poly{:x}(coeffs)
        @test 2.5*p == Poly{:x}(2.5*coeffs) == p*2.5
        @test p/2 == Poly{:x}(coeffs/2) == 2\p
        @test p*p == x^10+x^6+x^2-2*x^8+2*x^6-2*x^4 == (x^10+x^6+x^2+1)-2*(x^8-x^6+x^4)-1
        p = x + x^3 + x^5
        for i = 0:6
            @test p[i] == isodd(i)
        end
        @test p + p/2 == 1.5*p
        @test p*one(p) == p == p*one(x)
        @test p*zero(p) == 0 == p*zero(x)
        @test p <= p*x
        @test x*p == p*x
    end

    @testset "cost methods" begin

    end
end

@testset "cache datastructure" begin
    c = TensorOperations.LRU{Int,Int}(; maxsize = 1024) # 1 kilobyte cache
    numints = div(1024,sizeof(Int))
    sizehint!(c, numints)
    for i = 1:10
        c[i] = i
    end
    @test Set(collect(c)) == Set([i=>i for i = 1:10])
    @test haskey(c, 5)
    @test c[5] == 5
    c[5] = 6
    @test c[5] == 6
    @test !haskey(c, 11)
    @test length(c) == 10
    @test !haskey(c, 11)
    @test get!(()->11, c, 10) == 10
    @test c[10] == 10
    @test get!(()->11, c, 11) == 11
    @test haskey(c, 11)
    @test c[11] == 11
    @test !haskey(c, 12)
    @test get!(c, 10, 12) == 10
    @test c[10] == 10
    @test get!(c, 12, 12) == 12
    @test haskey(c, 12)
    @test c[12] == 12
    delete!(c, 12)
    @test !haskey(c, 12)
    @test length(c) == 11
    @test !isempty(c)
    empty!(c)
    @test isempty(c)
    @test length(c) == 0
    for i = 1:(numints+1)
        c[i] = i
    end
    @test !haskey(c, 1)
    @test haskey(c, 2)
    TensorOperations.setsize!(c; maxsize = 2048)
    for i = 1:(numints+1)
        c[i] = i
    end
    @test haskey(c, 1)
end


text:
using ImplicitDomainQuadrature
import Base.==, Base.≈

function allapprox(v1, v2)
    return all(v1 .≈ v2)
end

function allapprox(v1, v2, atol)
    return length(v1) == length(v2) &&
           all([isapprox(v1[i], v2[i], atol = atol) for i = 1:length(v1)])
end

function allequal(v1, v2)
    return all(v1 .== v2)
end

function Base.isequal(c1::CutCell.CellMap, c2::CutCell.CellMap)
    return allequal(c1.yL, c2.yL) && allequal(c1.yR, c2.yR)
end

function ==(c1::CutCell.CellMap, c2::CutCell.CellMap)
    return isequal(c1, c2)
end

function ≈(q1::QuadratureRule, q2::QuadratureRule)
    flag = allapprox(q1.points, q2.points) && allapprox(q1.weights, q2.weights)
end

function plane_distance_function(coords, normal, x0)
    return (coords .- x0)' * normal
end

function plane_distance_function(coords, normal, x0)
    return (coords .- x0)' * normal
end

function circle_distance_function(coords, center, radius)
    difference = (coords .- center) .^ 2
    distance = radius .- sqrt.(vec(mapslices(sum, difference, dims = 1)))
    return distance
end

function corner_distance_function(x::V, xc) where {V<:AbstractVector}
    v = xc - x
    if all(x .<= xc)
        minimum(v)
    elseif all(x .> xc)
        return -sqrt(v' * v)
    elseif x[2] > xc[2]
        return v[2]
    else
        return v[1]
    end
end

function corner_distance_function(points::M, xc) where {M<:AbstractMatrix}
    return vec(
        mapslices(x -> corner_distance_function(x, xc), points, dims = 1),
    )
end

function normal_from_angle(theta)
    return [cosd(theta), sind(theta)]
end

function add_cell_error_squared!(
    err,
    interpolater,
    exactsolution,
    cellmap,
    quad,
)
    detjac = CutCell.determinant_jacobian(cellmap)
    for (p, w) in quad
        numsol = interpolater(p)
        exsol = exactsolution(cellmap(p))
        err .+= (numsol - exsol) .^ 2 * detjac * w
    end
end

function add_interface_error_squared!(
    err,
    interpolater,
    exactsolution,
    cellmap,
    quad,
    facescale,
)
    @assert length(quad) == length(facescale)
    for (i, (p, w)) in enumerate(quad)
        numsol = interpolater(p)
        exsol = exactsolution(cellmap(p))
        err .+= (numsol - exsol) .^ 2 * facescale[i] * w
    end
end

function add_cell_norm_squared!(vals, func, cellmap, quad)
    detjac = CutCell.determinant_jacobian(cellmap)
    for (p, w) in quad
        v = func(cellmap(p))
        vals .+= v .^ 2 * detjac * w
    end
end

function add_interface_norm_squared!(vals, func, cellmap, quad, facescale)
    @assert length(quad) == length(facescale)
    for (i, (p, w)) in enumerate(quad)
        v = func(cellmap(p))
        vals .+= v .^ 2 * facescale[i] * w
    end
end

function integral_norm_on_cut_mesh(func, cellquads, cutmesh, ndofs)
    vals = zeros(ndofs)
    ncells = CutCell.number_of_cells(cutmesh)
    for cellid = 1:ncells
        s = CutCell.cell_sign(cutmesh, cellid)
        cellmap = CutCell.cell_map(cutmesh, cellid)
        @assert s == -1 || s == 0 || s == 1
        if s == 1 || s == 0
            pquad = cellquads[1, cellid]
            add_cell_norm_squared!(vals, func, cellmap, pquad)
        end
        if s == -1 || s == 0
            nquad = cellquads[-1, cellid]
            add_cell_norm_squared!(vals, func, cellmap, nquad)
        end
    end
    return sqrt.(vals)
end

function integral_norm_on_uniform_mesh(func, quad, cellmaps, ndofs)
    vals = zeros(ndofs)
    for cellmap in cellmaps
        add_cell_norm_squared!(vals, func, cellmap, quad)
    end
    return sqrt.(vals)
end

function integral_norm_on_interface(
    func,
    interfacequads,
    levelsetsign,
    cutmesh,
    ndofs,
)
    vals = zeros(ndofs)
    cellsign = CutCell.cell_sign(cutmesh)
    cellids = findall(cellsign .== 0)
    for cellid in cellids
        cellmap = CutCell.cell_map(cutmesh, cellid)
        quad = interfacequads[levelsetsign, cellid]
        normals = CutCell.interface_normals(interfacequads, cellid)
        facescale = CutCell.scale_area(cellmap, normals)

        add_interface_norm_squared!(vals, func, cellmap, quad, facescale)

    end
    return sqrt.(vals)
end

function interface_L2_error(
    nodalsolutions,
    exactsolution,
    levelsetsign,
    basis,
    interfacequads,
    cutmesh,
)
    ndofs = size(nodalsolutions)[1]
    err = zeros(ndofs)
    interpolater = InterpolatingPolynomial(ndofs, basis)

    cellsign = CutCell.cell_sign(cutmesh)
    cellids = findall(cellsign .== 0)

    for cellid in cellids
        quad = interfacequads[levelsetsign, cellid]
        normals = CutCell.interface_normals(interfacequads, cellid)
        cellmap = CutCell.cell_map(cutmesh, cellid)
        facescale = CutCell.scale_area(cellmap, normals)

        nodeids = CutCell.nodal_connectivity(cutmesh, levelsetsign, cellid)
        elementsolution = nodalsolutions[:, nodeids]
        update!(interpolater, elementsolution)

        add_interface_error_squared!(
            err,
            interpolater,
            exactsolution,
            cellmap,
            quad,
            facescale,
        )
    end
    return sqrt.(err)
end

function mesh_L2_error(nodalsolutions, exactsolution, basis, cellquads, cutmesh)
    ndofs = size(nodalsolutions)[1]
    err = zeros(ndofs)
    interpolater = InterpolatingPolynomial(ndofs, basis)
    ncells = CutCell.number_of_cells(cutmesh)
    for cellid = 1:ncells
        s = CutCell.cell_sign(cutmesh, cellid)
        @assert s == -1 || s == 0 || s == 1
        if s == 1 || s == 0
            cellmap = CutCell.cell_map(cutmesh, +1, cellid)
            nodeids = CutCell.nodal_connectivity(cutmesh, 1, cellid)
            elementsolution = nodalsolutions[:, nodeids]
            update!(interpolater, elementsolution)
            quad = cellquads[1, cellid]
            add_cell_error_squared!(
                err,
                interpolater,
                exactsolution,
                cellmap,
                quad,
            )
        end
        if s == -1 || s == 0
            cellmap = CutCell.cell_map(cutmesh, -1, cellid)
            nodeids = CutCell.nodal_connectivity(cutmesh, -1, cellid)
            elementsolution = nodalsolutions[:, nodeids]
            update!(interpolater, elementsolution)
            quad = cellquads[-1, cellid]
            add_cell_error_squared!(
                err,
                interpolater,
                exactsolution,
                cellmap,
                quad,
            )
        end
    end
    return sqrt.(err)
end

function uniform_mesh_L2_error(nodalsolutions, exactsolution, basis, quad, mesh)
    ndofs, nnodes = size(nodalsolutions)
    err = zeros(ndofs)
    interpolater = InterpolatingPolynomial(ndofs, basis)
    ncells = CutCell.number_of_cells(mesh)
    nodalconnectivity = CutCell.nodal_connectivity(mesh)

    for cellid = 1:ncells
        cellmap = CutCell.cell_map(mesh, cellid)
        nodeids = nodalconnectivity[:, cellid]
        elementsolution = nodalsolutions[:, nodeids]
        update!(interpolater, elementsolution)
        add_cell_error_squared!(err, interpolater, exactsolution, cellmap, quad)
    end
    return sqrt.(err)
end

function required_quadrature_order(polyorder)
    ceil(Int, 0.5 * (2polyorder + 1))
end

function mean(v)
    return sum(v) / length(v)
end

function convergence_rate(dx, err)
    return diff(log.(err)) ./ diff(log.(dx))
end


text:
using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
using ..ResultListUtils, ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils, ..ScanForDuplicates

resList = allocateResultList(100,100);
ind= (1,2,3)
getResLinIndex([ind[1] ind[2] ind[3] 1 1 1 ; 0 0 0 0 0 0],1, (200,200,200))

carts = CartesianIndices(zeros(200,200,200))

mappedA = map(ind-> getResLinIndex(ind[1],ind[2],ind[3],1, (200,200,200)) ,carts)
mappedB = map(ind-> getResLinIndex(ind[1],ind[2],ind[3],0, (200,200,200)) ,carts)

@test length(unique(mappedA))== 200*200*200
@test length(unique(vec(vcat(mappedA,mappedB ))))== 200*200*200*2

#################  addResult
using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
using ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils,..ResultListUtils

threads=(32,5);
blocks =2;
mainArrDims= (67,177,90);
dataBdim = (43,21,17);
mainArrCPU= falses(mainArrDims);
mainArrCPU[5,5,5]= true;
mainArrCPU[33,10,7]= true;
mainArrGPU = CuArray(mainArrCPU);
metaData = MetaDataUtils.allocateMetadata(mainArrDims,dataBdim);
metaDataDims= size(metaData);
totalFp,totalFn= 500,500
resList,resListIndicies,maxResListIndex= allocateResultLists(totalFp,totalFn)

xMeta,yMeta,zMeta= 1,1,1 #aaume 0 based
x,y,z = 1,2,3
dir= 5
queueNumber = 1
isGold=4
iterNumb=6
metaData[xMeta+1,yMeta+1,zMeta+1,getResOffsetsBeg()+queueNumber]=7
function addResKernel(metaData ,xMeta,yMeta,zMeta, resList,resListIndicies,maxResListIndex,x,y,z, dir,iterNumb,queueNumber,metaDataDims ,mainArrDims,isGold )
    @ifXY 1 1 @addResult(metaData ,xMeta,yMeta,zMeta, resList,resListIndicies,maxResListIndex,x,y,z, dir,iterNumb,queueNumber,metaDataDims,mainArrDims ,isGold  )
    @ifXY 2 1 @addResult(metaData ,xMeta,yMeta,zMeta, resList,resListIndicies,maxResListIndex,x,y,z, dir,iterNumb,queueNumber,metaDataDims ,mainArrDims,isGold  )
    @ifXY 3 1 @addResult(metaData ,xMeta,yMeta,zMeta, resList,resListIndicies,maxResListIndex,x,y,z, dir,iterNumb,queueNumber,metaDataDims,mainArrDims ,isGold  )
    return
end

@cuda threads=threads blocks=blocks addResKernel(metaData ,xMeta,yMeta,zMeta, resList,resListIndicies,maxResListIndex,x,y,z, dir,iterNumb,queueNumber,metaDataDims ,mainArrDims,isGold )
@test length(filter(it-> it>0,Array(resListIndicies)))==6

# @test  Int64.(Array(resList[7,:]))==[1,2,3,4,5,6]
@test  Int64.(Array(resList[8,:]))==[1,2,3,4,5,6]
@test  Int64.(Array(resList[9,:]))==[1,2,3,4,5,6]
@test  Int64.(Array(resList[10,:]))==[1,2,3,4,5,6]
@test  Int64.(Array(resList[11,:]))==[1,2,3,4,5,6]
@test  Int64.(Array(resList[12,:]))==[1,2,3,4,5,6]
@test  Int64.(Array(resList[14,:]))==[0,0,0,0,0,0]

Int64.(Array(resListIndicies[5:15]))

filter(it-> it>0,Array(resListIndicies))
Int64(sum(resList))


text:
"""
    V_test(θ, θ₀, degrees=false) -> is_significant, V, Vcrit

Perform the modified Rayleigh or V test to determine if a sample of angles
`θ` differ significantly from random towards the angle `θ₀` at the `α` level.

Return `true` if they do differ from random significantly, the test
statistic `V` and the critical value `Vcrit`.
"""
function V_test(θ, θ₀, α=0.05, degrees=false)
    degrees && (θ = deg2rad.(θ); θ₀ = deg2rad(θ₀))
    n = length(θ)
    r = cresultant(θ)
    θ̄ = cmean(θ)
    ϑ = r*cos(θ̄ .- θ₀)
    V = sqrt(2n)*ϑ
    Vcrit = V_critical_value(n, α)
    V > Vcrit, V, Vcrit
end

"""Table of critical values of the V (modified Rayleigh) test.  Table 34 of Kanji (2006)"""
const V_TEST_TABLE = [1.3051 1.6524 2.2505 2.4459 2.7938 3.0825
                      1.3009 1.6509 2.2640 2.4695 2.8502 3.2114
                      1.2980 1.6499 2.2734 2.4858 2.8886 3.2970
                      1.2958 1.6492 2.2803 2.4978 2.9164 3.3578
                      1.2942 1.6484 2.2856 2.5070 2.9375 3.4034
                      1.2929 1.6482 2.2899 2.5143 2.9540 3.4387
                      1.2918 1.6479 2.2933 2.5201 2.9672 3.4669
                      1.2909 1.6476 2.2961 2.5250 2.9782 3.4899
                      1.2902 1.6474 2.2985 2.5290 2.9873 3.5091
                      1.2895 1.6472 2.3006 2.5325 2.9950 3.5253
                      1.2890 1.6470 2.3023 2.5355 3.0017 3.5392
                      1.2885 1.6469 2.3039 2.5381 3.0075 3.5512
                      1.2881 1.6467 2.3052 2.5404 3.0126 3.5617
                      1.2877 1.6466 2.3064 2.5424 3.0171 3.5710
                      1.2874 1.6465 2.3075 2.5442 3.0211 3.5792
                      1.2871 1.6464 2.3085 2.5458 3.0247 3.5866
                      1.2868 1.6464 2.3093 2.5473 3.0279 3.5932
                      1.2866 1.6463 2.3101 2.5486 3.0308 3.5992
                      1.2864 1.6462 2.3108 2.5498 3.0335 3.6047
                      1.2862 1.6462 2.3115 2.5509 3.0359 3.6096
                      1.2860 1.6461 2.3121 2.5519 3.0382 3.6142
                      1.2858 1.6461 2.3127 2.5529 3.0402 3.6184
                      1.2856 1.6460 2.3132 2.5538 3.0421 3.6223
                      1.2855 1.6460 2.3136 2.5546 3.0439 3.6258
                      1.2853 1.6459 2.3141 2.5553 3.0455 3.6292
                      1.2852 1.6459 2.3145 2.5560 3.0471 3.6323
                      1.2843 1.6456 2.3175 2.5610 3.0580 3.6545
                      1.2837 1.6455 2.3193 2.5640 3.0646 3.6677
                      1.2834 1.6454 2.3205 2.5660 3.0689 3.6764
                      1.2831 1.6453 2.3213 2.5674 3.0720 3.6826
                      1.2826 1.6452 2.3228 2.5699 3.0775 3.6936
                      1.2818 1.6449 2.3256 2.5747 3.0877 3.7140
                      1.2817 1.6449 2.3260 2.5752 3.0890 3.7165]
"Values of α (columns) for V_TEST_TABLE"
const V_TEST_ALPHAS = [0.1000, 0.0500, 0.0100, 0.0050, 0.0010, 0.0001]
"Values of n (rows) for V_TEST_TABLE"
const V_TEST_NS = [5:30..., 40, 50, 60, 70, 100, 500, 1000]

"""
    V_critical_value(n, α) -> Vcrit

Return the critical value for the V test `Vcrit` for a sample size of `n` at the
`α` confidence level.
"""
function V_critical_value(n, α)
    n >= V_TEST_NS[1] || throw(ArgumentError("V test tables only valid for `n >= 5` (have $n)"))
    indn = argmin(abs.(n .- V_TEST_NS))
    indα = argmin(abs.(α .- V_TEST_ALPHAS))
    V_TEST_TABLE[indn,indα]
end

"""
    watson_U2n(θ, cdf, α=0.05, degrees=false; axial=false) -> test, U², U²crit

Perform Watson's U²n test to determine if a sample of angles `θ`
fit the supplied cdf of a distribution at the `α` level of
significance.  *`test` is `true` if the sample differs significantly*.
Return also the value of the test statistic `U²` and the critical value
`U²crit`.

    watson_U2n(θ, α=0.05, degrees=false; axial=false) -> test, U², U²crit

Perform the test against the maximum-likelihood von Mises distribution
for the sample data.  This is fit using `fit_vonMises()`.
"""
function watson_U2n(θ, cdf::Function, α=0.05, degrees=false; axial=false)
    axial && (θ = 2θ)
    degrees && (θ = deg2rad.(θ))
    θ = sort(mod.(θ, 2π))
    n = length(θ)
    V = cdf.(θ) .- cdf(0)
    V̄ = sum(V)/n
    U² = sum(V.^2) - sum((2*(1:n) .- 1).*V./n) + n*(1/3 - (V̄ - 1/2)^2)
    U²crit = watson_U2n_crit(n, α)
    U² > U²crit, U², U²crit
end

function watson_U2n(θ, α::Number=0.05, degrees=false; kwargs...)
    μ, κ = fit_vonMises(θ, degrees; kwargs...)
    watson_U2n(θ, x->von_mises_cdf(x, μ, κ, degrees), α, degrees; kwargs...)
end

"Table of critical values of Watson's U²n test.  Table 35 of Kanji (2006)."
const WATSON_U2N_TABLE = [0.143 0.000 0.161 0.164 0.165
                          0.145 0.173 0.194 0.213 0.224
                          0.146 0.176 0.202 0.233 0.252
                          0.148 0.177 0.205 0.238 0.262
                          0.149 0.179 0.208 0.243 0.269
                          0.149 0.180 0.210 0.247 0.274
                          0.150 0.181 0.211 0.250 0.278
                          0.150 0.182 0.212 0.252 0.281
                          0.150 0.182 0.213 0.254 0.283
                          0.150 0.183 0.215 0.256 0.287
                          0.151 0.184 0.216 0.258 0.290
                          0.151 0.184 0.216 0.259 0.291
                          0.151 0.184 0.217 0.259 0.292
                          0.151 0.185 0.217 0.261 0.293
                          0.152 0.185 0.219 0.263 0.296
                          0.152 0.186 0.219 0.264 0.298
                          0.152 0.186 0.220 0.265 0.299
                          0.152 0.186 0.221 0.266 0.301
                          0.152 0.187 0.221 0.267 0.302]
"Values of α (columns) for WATSON_U2N_TABLE"
const WATSON_U2N_ALPHAS = [0.100, 0.050, 0.025, 0.010, 0.005]
"Values of n (rows) for WATSON_U2N_table"
const WATSON_U2N_N = [2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 30, 40, 50, 100, 200]

"""
    watson_U2n_crit(n, α) -> U²crit

Return the critical value at the α level for a sample of `n` angles.
"""
function watson_U2n_crit(n, α)
    2 <= n || throw(ArgumentError("`n` (supplied $n) must be 2 or more"))
    any(isapprox.(α, WATSON_U2N_ALPHAS, atol=0.0001)) == 1 ||
        throw(ArgumentError("Significance level for Watson's U²n test must " *
                            "be one of $(WATSON_U2N_ALPHAS); asked for $α"))
    indn = argmin(abs.(n .- WATSON_U2N_N))
    indα = argmin(abs.(α .- WATSON_U2N_ALPHAS))
    WATSON_U2N_TABLE[indn,indα]
end

"""
    watsonU2(θ, ϕ, α=0.05, degrees=false)

Perform Watson's U² test to determine if two samples `θ` and `ϕ` differ
significantly from each other in terms of mean direction or variance,
assuming both are sampled continuously.  Kanji (2006) recommends avoiding
this test if data are grouped more coarsely than 5°.
"""
function watson_U2(θ, ϕ, α=0.05, degrees=false)
    verbose = false
    degrees && (θ = deg2rad.(θ); ϕ = deg2rad.(ϕ))
    n, m = length(θ), length(ϕ)
    N = n + m
    θ = Float64.(sort(mod.(θ, 2π)))
    ϕ = Float64.(sort(mod.(ϕ, 2π)))
    Σdk, Σdk² = 0.0, 0.0
    i = j = 0
    verbose && println("  i    θi    i/n   j     ϕj   j/n     dk    dk²")
    while true
        i == n && j == m && break
        if i == n
            j += 1
        elseif j == m
            i += 1
        elseif θ[i+1] < ϕ[j+1]
            i += 1
        elseif θ[i+1] > ϕ[j+1]
            j += 1
        else
            i += 1
            j += 1
        end
        dk = i/n - j/m
        dk² = dk^2
        Σdk += dk
        Σdk² += dk²
        if verbose
            θi, ϕj = (i > 0 ? θ[i] : 0), (j>0 ? ϕ[j] : 0)
            @printf("%3d %6.2f %5.2f %3d %6.2f %5.2f %6.2f %6.2f\n",
                i, rad2deg(θi), θi/n, j, rad2deg(ϕj), ϕj/m, dk, dk²)
        end
    end
    verbose && @show Σdk, Σdk²
    i == n && j == m || error("Unexpected mismatch between `i`, `j`, `n` and `m`")
    U² = n*m/N^2*(Σdk² - Σdk^2/N)
    verbose && @show U²
    # U²crit = watson_U2_crit(n, m, α)
    U²crit = watson_U2_crit(n, m, α)
    U² > U²crit, U², U²crit
end

"Values of the larger sample size for WATSON_U2_TABLE"
const WATSON_U2_TABLE_NS = [10, 12, 14, 16, 20, 25, 30, 40, 50, 100]
"Values of α for WATSON_U2_TABLE"
const WATSON_U2_TABLE_ALPHAS = [0.1, 0.05, 0.01, 0.001]
function make_Watson_U2_table()
    t = Dict{Int, Dict{Int,NTuple{length(WATSON_U2_TABLE_ALPHAS),Float64}}}()
    for n1 in WATSON_U2_TABLE_NS
        t[n1] = Dict{Int,NTuple{length(WATSON_U2_TABLE_ALPHAS),Float64}}()
    end
    t[10][6]  = (0.153, 0.182, 0.243, 0.266)
    t[10][7]  = (0.153, 0.183, 0.246, 0.270)
    t[10][8]  = (0.153, 0.183, 0.248, 0.273)
    t[10][9]  = (0.153, 0.184, 0.249, 0.275)
    t[10][10] = (0.153, 0.184, 0.250, 0.277)

    t[12][6]  = (0.153, 0.182, 0.244, 0.268)
    t[12][7]  = (0.153, 0.183, 0.247, 0.272)
    t[12][8]  = (0.153, 0.184, 0.249, 0.275)
    t[12][9]  = (0.153, 0.184, 0.251, 0.278)
    t[12][10] = (0.153, 0.184, 0.252, 0.279)
    t[12][11] = (0.153, 0.184, 0.253, 0.280)
    t[12][12] = (0.153, 0.184, 0.253, 0.282)

    t[14][6]  = (0.153, 0.182, 0.244, 0.269)
    t[14][8]  = (0.153, 0.184, 0.250, 0.276)
    t[14][10] = (0.153, 0.184, 0.253, 0.281)
    t[14][12] = (0.153, 0.185, 0.254, 0.283)
    t[14][14] = (0.153, 0.185, 0.255, 0.285)

    t[16][4]  = (0.151, 0.178, 0.231, 0.251)
    t[16][6]  = (0.152, 0.182, 0.244, 0.269)
    t[16][8]  = (0.153, 0.183, 0.250, 0.277)
    t[16][10] = (0.153, 0.184, 0.253, 0.282)
    t[16][12] = (0.153, 0.185, 0.255, 0.284)
    t[16][14] = (0.153, 0.185, 0.256, 0.286)
    t[16][16] = (0.153, 0.185, 0.257, 0.287)

    t[20][4]  = (0.151, 0.178, 0.231, 0.251)
    t[20][6]  = (0.152, 0.182, 0.245, 0.270)
    t[20][8]  = (0.152, 0.183, 0.251, 0.278)
    t[20][12] = (0.153, 0.185, 0.256, 0.285)
    t[20][16] = (0.153, 0.185, 0.258, 0.289)
    t[20][20] = (0.153, 0.185, 0.259, 0.290)

    t[25][5]  = (0.151, 0.180, 0.240, 0.263)
    t[25][10] = (0.152, 0.184, 0.254, 0.283)
    t[25][15] = (0.153, 0.185, 0.258, 0.289)
    t[25][20] = (0.153, 0.186, 0.260, 0.291)
    t[25][26] = (0.153, 0.186, 0.261, 0.293)

    t[30][5]  = (0.150, 0.179, 0.240, 0.263)
    t[30][10] = (0.152, 0.184, 0.254, 0.284)
    t[30][15] = (0.152, 0.185, 0.259, 0.290)
    t[30][30] = (0.153, 0.186, 0.262, 0.294)

    t[40][5]  = (0.150, 0.179, 0.240, 0.263)
    t[40][10] = (0.152, 0.184, 0.254, 0.284)
    t[40][15] = (0.152, 0.185, 0.259, 0.290)
    t[40][20] = (0.152, 0.186, 0.261, 0.293)
    t[40][40] = (0.152, 0.186, 0.263, 0.296)

    t[50][5]  = (0.150, 0.178, 0.239, 0.263)
    t[50][10] = (0.151, 0.183, 0.254, 0.284)
    t[50][15] = (0.152, 0.185, 0.259, 0.290)
    t[50][20] = (0.152, 0.185, 0.261, 0.293)
    t[50][25] = (0.152, 0.186, 0.262, 0.295)
    t[50][50] = (0.152, 0.186, 0.264, 0.298)

    # Actually t[∞,∞]
    t[100][100] = (0.1517, 0.1869, 0.2684, 0.3035)

    t
end

"""Table of critical values for Watson's two-sample U² test.  Table 9 of
Batschelet (1972).

#### Reference
- Batschelet, E (1972).  Recent statistical methods for orientation data.
  Animal Orientation and Navigation, NASA SP262, 61-91.
"""
const WATSON_U2_TABLE = make_Watson_U2_table()

function watson_U2_crit(n, m, α)
    n1 = max(n, m)
    n2 = min(n, m)
    n1ind = argmin(abs.(n1 .- WATSON_U2_TABLE_NS))
    n1 = WATSON_U2_TABLE_NS[n1ind]
    any(isapprox.(α, WATSON_U2_TABLE_ALPHAS, atol=0.0001)) ||
        throw(ArgumentError("Significance level for Watson's U² test must " *
                            "be one of $(WATSON_U2_TABLE_ALPHAS); asked for $α"))
    αind = argmin(abs.(α .- WATSON_U2_TABLE_ALPHAS))
    t = WATSON_U2_TABLE[n1]
    for k in sort(collect(keys(t)))
        n2 <= k && return t[k][αind]
    end
    return t[100][αind]
end


text:
# # Flock model

using Agents, Random, LinearAlgebra

mutable struct Bird <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    vel::NTuple{2,Float64}
    speed::Float64
    cohere_factor::Float64
    separation::Float64
    seperate_factor::Float64
    match_factor::Float64
    visual_distance::Float64
end

# The function `initialize_model` generates birds and returns a model object using default values.
function initialize_model(;n_birds=100, speed=1.0, cohere_factor=.25, separation=4.0,
    seperate_factor=.25, match_factor=.01, visual_distance=5.0, dims=(100,100), seed=0)
    Random.seed!(seed)

    space2d = ContinuousSpace(2; periodic=true, extend=dims)
    model = ABM(Bird, space2d, scheduler=random_activation)
    for _ in 1:n_birds
        vel = Tuple(rand(2)*2 .- 1)
        add_agent!(
            model, vel, speed, cohere_factor,separation, seperate_factor,
            match_factor,visual_distance
        )
    end
    index!(model)
    return model
end

# ## Defining the agent_step!
# `agent_step!` is the primary function called for each step and computes velocity
# according to the three rules defined above.
function agent_step!(bird, model)
    # Obtain the ids of neibhors within the bird's visual distance
    ids = space_neighbors(bird, model, bird.visual_distance)
    # Compute velocity based on rules defined above
    bird.vel = (bird.vel .+ cohere(bird, model, ids) .+ seperate(bird, model, ids)
        .+ match(bird, model, ids))./2
    bird.vel = bird.vel ./ norm(bird.vel)
    # Move bird according to new velocity and speed
    move_agent!(bird, model, bird.speed)
end

distance(a1, a2) = sqrt(sum((a1.pos .- a2.pos).^2))

get_heading(a1, a2) = a1.pos .- a2.pos

# cohere computes the average position of neighboring birds, weighted by importance
function cohere(bird, model, ids)
    N = max(length(ids), 1)
    coherence = (0.0,0.0)
    for id in ids
        coherence = coherence .+ get_heading(model[id], bird)
    end
    return coherence ./ N .* bird.cohere_factor
end

# seperate repells the bird away from neighboring birds
function seperate(bird, model, ids)
    seperation_vec = (0.0,0.0)
    N = max(length(ids), 1)
    for id in ids
        if distance(bird, model[id]) < bird.separation
            seperation_vec = seperation_vec .- get_heading(model[id], bird)
        end
    end
    return seperation_vec./N.*bird.seperate_factor
end

# match computes the average trajectory of neighboring birds, weighted by importance
function match(bird, model, ids)
    match_vector = (0.0,0.0)
    N = max(length(ids), 1)
    for id in ids
        match_vector = match_vector .+ model[id].vel
    end
    return match_vector ./ N .* bird.match_factor
end


text: a ? b : c

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode