text:
@testset "UnlabeledSlidingWindow" begin
    @test_throws UndefVarError UnlabeledSlidingWindow
    @test MLDataPattern.UnlabeledSlidingWindow <: AbstractVector
    @test MLDataPattern.UnlabeledSlidingWindow <: LearnBase.DataView
    @test !(MLDataPattern.UnlabeledSlidingWindow <: LearnBase.AbstractObsIterator)
    @test !(MLDataPattern.UnlabeledSlidingWindow <: LearnBase.AbstractBatchIterator)

    @testset "constructor" begin
        @test_throws DimensionMismatch slidingwindow((rand(2,10),rand(9)), 1)
        @test_throws DimensionMismatch slidingwindow((rand(2,10),rand(4,9,10),rand(9)), 1)
        @test_throws MethodError slidingwindow(EmptyType(), 1)
        @test_throws MethodError slidingwindow(EmptyType(), 1, ObsDim.Last())
        @test_throws MethodError slidingwindow(EmptyType(), 1, ObsDim.Undefined())
        @test_throws MethodError slidingwindow(EmptyType(), 1, ObsDim.Last())
        @test_throws MethodError slidingwindow((EmptyType(),EmptyType()))
        @test_throws MethodError slidingwindow(CustomType(), 1, ObsDim.Last())
        @test_throws MethodError slidingwindow(EmptyType(), 1, obsdim=1)
        @test_throws MethodError slidingwindow((EmptyType(),EmptyType()), 1)

        for var in (vars..., tuples..., Xs, ys)
            @test_throws MethodError slidingwindow(var...)
            @test_throws MethodError slidingwindow(var..., 1)
            @test_throws MethodError slidingwindow(var..., obsdim=:last)
            @test_throws MethodError slidingwindow(var..., 1, obsdim=:last)
            @test_throws ArgumentError slidingwindow(var, 151)
            @test_throws ArgumentError slidingwindow(var, 0)
            @test_throws ArgumentError slidingwindow(var, 1, 0)
            A = @inferred slidingwindow(var, 5, 6)
            @test A == slidingwindow(var, 5, stride = 6)
            @test @inferred(length(A)) == A.count == 25
            @test A.size == 5
            @test A.stride == 6
            @test A.obsdim == LearnBase.default_obsdim(var)
            @test @inferred(nobs(A)) == length(A)
            @test @inferred(parent(A)) === var
            A = @inferred(slidingwindow(var, 4))
            @test A == slidingwindow(var, 4, stride = 4)
            @test length(A) == A.count == 37
            @test A.size == 4
            @test A.stride == 4
            @test A.obsdim == LearnBase.default_obsdim(var)
            @test nobs(A) == length(A)
            @test @inferred(parent(A)) === var
        end
        @test slidingwindow((X,X), 5) == @inferred(slidingwindow((X,X), 5, (ObsDim.Last(),ObsDim.Last())))
        @test slidingwindow((X,X), 5) == @inferred(slidingwindow((X,X), 5, 5, (ObsDim.Last(),ObsDim.Last())))
        A = slidingwindow(X',5,obsdim=1)
        @test A == @inferred(slidingwindow(X',5,ObsDim.First()))
        @test A == @inferred(slidingwindow(X',5,5,ObsDim.First()))
        @test A == slidingwindow(X',5,obsdim=:first)
        @test A == slidingwindow(X',5,stride=5,obsdim=:first)
    end

    println("<HEARTBEAT>")

    @testset "typestability" begin
        @test_throws ErrorException @inferred(slidingwindow(X, 5, obsdim=2))
        for var in (vars..., tuples..., Xs, ys)
            @test typeof(@inferred(slidingwindow(var, 5))) <: MLDataPattern.UnlabeledSlidingWindow
            @test typeof(@inferred(slidingwindow(var, 5, 3))) <: MLDataPattern.UnlabeledSlidingWindow
            @test typeof(@inferred(slidingwindow(var, 5, ObsDim.Last()))) <: MLDataPattern.UnlabeledSlidingWindow
            @test typeof(@inferred(slidingwindow(var, 5, 3, ObsDim.Last()))) <: MLDataPattern.UnlabeledSlidingWindow
            @test typeof(@inferred(slidingwindow(var, 5))) <: MLDataPattern.UnlabeledSlidingWindow
            @test_throws ErrorException @inferred(slidingwindow(var, 5, obsdim=:last))
        end
        for tup in tuples
            @test typeof(@inferred(slidingwindow(tup,5,(fill(ObsDim.Last(),length(tup))...,)))) <: MLDataPattern.UnlabeledSlidingWindow
            @test typeof(@inferred(slidingwindow(tup,5,3,(fill(ObsDim.Last(),length(tup))...,)))) <: MLDataPattern.UnlabeledSlidingWindow
        end
        @test typeof(@inferred(slidingwindow(CustomType(),5))) <: MLDataPattern.UnlabeledSlidingWindow
        @test typeof(@inferred(slidingwindow(CustomType(), 5, ObsDim.Undefined()))) <: MLDataPattern.UnlabeledSlidingWindow
    end

    println("<HEARTBEAT>")

    @testset "AbstractArray interface" begin
        for var in (vars..., tuples..., Xs, ys)
            A = slidingwindow(var, 5)
            @test_throws BoundsError A[-1]
            @test_throws BoundsError A[length(A)+1]
            @test @inferred(nobs(A)) == 30
            @test @inferred(length(A)) == 30
            @test @inferred(size(A)) == (30,)
            @test @inferred(A[2:3]) == [datasubset(var,(1*5+1):(2*5)), datasubset(var,(2*5+1):(3*5))]
            @test eltype(A[2:3]) == typeof(datasubset(var,(1*5+1):(2*5)))
            @test @inferred(A[[1,3]]) == [datasubset(var,(0*5+1):(1*5)), datasubset(var,(2*5+1):(3*5))]
            @test @inferred(A[1]) == datasubset(var, 1:5)
            @test @inferred(A[30]) == datasubset(var, (29*5+1):(30*5))
            @test @inferred(A[30]) == datasubset(var, nobs(var)-4:nobs(var))
            @test A[end] == A[30]
            @test @inferred(getobs(A,1)) == getobs(var, 1:5)
            @test @inferred(getobs(A,2:3)) == [getobs(var,(1*5+1):(2*5)), getobs(var,(2*5+1):(3*5))]
            if !(var isa Tuple)
                @test eltype(getobs(A,1)) == eltype(var)
            end
            @test typeof(getobs(A,1)) == typeof(getobs(var, 1:5))
            @test eltype(getobs(A,2:3)) == typeof(getobs(var, 1:5))
            @test typeof(@inferred(collect(A))) <: Vector
            # custom stride
            A = slidingwindow(var, 5, 2)
            @test_throws BoundsError A[-1]
            @test_throws BoundsError A[length(A)+1]
            @test @inferred(nobs(A)) == 73
            @test @inferred(length(A)) == 73
            @test @inferred(size(A)) == (73,)
            @test A[end] == A[73]
            @test @inferred(getobs(A,1)) == getobs(var, 1:5)
            @test @inferred(A[2:3]) == [datasubset(var,3:3+4), datasubset(var,5:5+4)]
            if !(var isa Tuple)
                @test eltype(getobs(A,1)) == eltype(var)
            end
            @test typeof(@inferred(collect(A))) <: Vector
        end
    end

    println("<HEARTBEAT>")

    @testset "batchview" begin
        bv = batchview(slidingwindow(1:100, 5))
        @test_throws MethodError shuffleobs(bv)
        @test_throws MethodError splitobs(bv)
        @test_throws MethodError oversample(bv)
        @test_throws MethodError undersample(bv)
        @test_throws MethodError stratifiedobs(bv)
        for var in vars # (vars..., tuples..., Xs, ys)
            A = slidingwindow(var, 5)
            bv = @inferred batchview(A, 2)
            go = @inferred getobs(bv)
            @test go isa Array
            @test go == bv
            @test length(bv) == length(A) / 2
            @test eltype(bv) <: Array{<:SubArray}
            bs = @inferred bv[1]
            @test length(bs) == 5
            @test all(nobs.(bs) .== 2)
        end
        A = batchview(slidingwindow(1:20, 2), 5)
        @test length(A) == 2
        @test A[1] == [[1,3,5,7,9],[2,4,6,8,10]]
        @test A[2] == [[11,13,15,17,19],[12,14,16,18,20]]
        @test A[1] isa Array
    end

    println("<HEARTBEAT>")

    @testset "special values" begin
        A = @inferred slidingwindow(1:10, 2)
        @test length(A) == 5
        @test getobs(A) == [
            [1, 2],
            [3, 4],
            [5, 6],
            [7, 8],
            [9, 10],
        ]
        A = @inferred slidingwindow(1:10, 3)
        @test length(A) == 3
        @test getobs(A) == [
            [1, 2, 3],
            [4, 5, 6],
            [7, 8, 9],
        ]
        A = @inferred slidingwindow(1:10, 3, 2)
        @test length(A) == 4
        @test getobs(A) == [
            [1, 2, 3],
            [3, 4, 5],
            [5, 6, 7],
            [7, 8, 9],
        ]
        A = @inferred slidingwindow(1:10, 3, 4)
        @test length(A) == 2
        @test getobs(A) == [
            [1, 2, 3],
            [5, 6, 7],
        ]
    end
end

@testset "LabeledSlidingWindow" begin
    @test_throws UndefVarError LabeledSlidingWindow
    @test MLDataPattern.LabeledSlidingWindow <: AbstractVector
    @test MLDataPattern.LabeledSlidingWindow <: LearnBase.DataView
    @test !(MLDataPattern.LabeledSlidingWindow <: LearnBase.AbstractObsIterator)
    @test !(MLDataPattern.LabeledSlidingWindow <: LearnBase.AbstractBatchIterator)

    @testset "constructor" begin
        @test_throws DimensionMismatch slidingwindow(i->i, (rand(2,10),rand(9)), 1)
        @test_throws DimensionMismatch slidingwindow(i->i, (rand(2,10),rand(4,9,10),rand(9)), 1)
        @test_throws MethodError slidingwindow(i->i, EmptyType(), 1)
        @test_throws MethodError slidingwindow(i->i, EmptyType(), 1, ObsDim.Last())
        @test_throws MethodError slidingwindow(i->i, EmptyType(), 1, ObsDim.Undefined())
        @test_throws MethodError slidingwindow(i->i, EmptyType(), 1, ObsDim.Last())
        @test_throws MethodError slidingwindow(i->i, (EmptyType(),EmptyType()))
        @test_throws MethodError slidingwindow(i->i, CustomType(), 1, ObsDim.Last())
        @test_throws MethodError slidingwindow(i->i, EmptyType(), 1, obsdim=1)
        @test_throws MethodError slidingwindow(i->i, (EmptyType(),EmptyType()), 1)

        for var in (vars..., tuples..., Xs, ys)
            @test_throws MethodError slidingwindow(i->i, var...)
            @test_throws MethodError slidingwindow(i->i, var..., 1)
            @test_throws MethodError slidingwindow(i->i, var..., obsdim=:last)
            @test_throws MethodError slidingwindow(i->i, var..., 1, obsdim=:last)
            @test_throws ArgumentError slidingwindow(i->i, var, 151)
            @test_throws ArgumentError slidingwindow(i->i, var, 0)
            @test_throws ArgumentError slidingwindow(i->i, var, 1, 0)
            A = @inferred slidingwindow(identity, var, 5, 6)
            @test A == slidingwindow(i->i, var, 5, stride = 6)
            @test A == slidingwindow(i->i, var, 5, stride = 6, excludetarget=false)
            @test typeof(A) == typeof(slidingwindow(identity, var, 5, stride = 6, excludetarget=false))
            @test typeof(A) != typeof(slidingwindow(identity, var, 5, stride = 6, excludetarget=true))
            @test @inferred(length(A)) == A.count == 25
            @test A.size == 5
            @test A.stride == 6
            @test A.offset == 0
            @test A.obsdim == LearnBase.default_obsdim(var)
            @test @inferred(nobs(A)) == length(A)
            @test @inferred(parent(A)) === var
            A = @inferred(slidingwindow(i->i, var, 4))
            @test A == slidingwindow(i->i, var, 4, stride = 4)
            @test length(A) == A.count == 37
            @test A.size == 4
            @test A.stride == 4
            @test A.offset == 0
            @test A.obsdim == LearnBase.default_obsdim(var)
            @test nobs(A) == length(A)
            @test @inferred(parent(A)) === var
        end
        @test slidingwindow(i->i,(X,X), 5) == @inferred(slidingwindow(i->i,(X,X), 5, (ObsDim.Last(),ObsDim.Last())))
        @test slidingwindow(i->i,(X,X), 5) == @inferred(slidingwindow(i->i,(X,X), 5, Val(false)))
        @test slidingwindow(i->i,(X,X), 5) == @inferred(slidingwindow(i->i,(X,X), 5, 5, (ObsDim.Last(),ObsDim.Last())))
        @test slidingwindow(i->i,(X,X), 5) == @inferred(slidingwindow(i->i,(X,X), 5, 5, Val(false)))
        A = slidingwindow(i->i-1,X',5,obsdim=1)
        @test A.size == 5
        @test A.stride == 5
        @test A.offset == 1
        @test A.count == 29
        @test A.obsdim == LearnBase.ObsDim.First()
        @test A == @inferred(slidingwindow(i->i-1,X',5,ObsDim.First()))
        @test A == @inferred(slidingwindow(i->i-1,X',5,5,ObsDim.First()))
        @test A == slidingwindow(i->i-1,X',5,obsdim=:first)
        @test A == slidingwindow(i->i-1,X',5,stride=5,obsdim=:first)
    end

    println("<HEARTBEAT>")

    @testset "typestability" begin
        @test_throws ErrorException @inferred(slidingwindow(identity, X, 5, obsdim=2))
        for var in (vars..., tuples..., Xs, ys)
            @test typeof(@inferred(slidingwindow(identity, var, 5))) <: MLDataPattern.LabeledSlidingWindow
            @test typeof(@inferred(slidingwindow(identity, var, 5, 3))) <: MLDataPattern.LabeledSlidingWindow
            @test typeof(@inferred(slidingwindow(identity, var, 5, ObsDim.Last()))) <: MLDataPattern.LabeledSlidingWindow
            @test typeof(@inferred(slidingwindow(identity, var, 5, 3, ObsDim.Last()))) <: MLDataPattern.LabeledSlidingWindow
            @test typeof(@inferred(slidingwindow(identity, var, 5))) <: MLDataPattern.LabeledSlidingWindow
            @test_throws ErrorException @inferred(slidingwindow(identity, var, 5, obsdim=:last))
        end
        for tup in tuples
            @test typeof(@inferred(slidingwindow(identity,tup,5,(fill(ObsDim.Last(),length(tup))...,)))) <: MLDataPattern.LabeledSlidingWindow
            @test typeof(@inferred(slidingwindow(identity,tup,5,3,(fill(ObsDim.Last(),length(tup))...,)))) <: MLDataPattern.LabeledSlidingWindow
        end
        @test typeof(@inferred(slidingwindow(identity,CustomType(),5))) <: MLDataPattern.LabeledSlidingWindow
        @test typeof(@inferred(slidingwindow(identity,CustomType(), 5, ObsDim.Undefined()))) <: MLDataPattern.LabeledSlidingWindow
    end

    println("<HEARTBEAT>")

    @testset "AbstractArray interface" begin
        for var in (vars..., tuples..., Xs, ys)
            A = slidingwindow(i->i+5, var, 5)
            @test_throws BoundsError A[-1]
            @test_throws BoundsError A[length(A)+1]
            @test @inferred(nobs(A)) == 29
            @test @inferred(length(A)) == 29
            @test @inferred(size(A)) == (29,)
            ref = [
                (datasubset(var,(1*5+1):(2*5)), datasubset(var,2*5+1)),
                (datasubset(var,(2*5+1):(3*5)), datasubset(var,3*5+1))
            ]
            @test @inferred(A[2:3]) == ref
            @test eltype(A[2:3]) == typeof(ref[1])
            ref = [
                (datasubset(var,(0*5+1):(1*5)), datasubset(var,1*5+1)),
                (datasubset(var,(2*5+1):(3*5)), datasubset(var,3*5+1))
            ]
            @test @inferred(A[[1,3]]) == ref
            ref = (datasubset(var,(0*5+1):(1*5)), datasubset(var,1*5+1))
            @test @inferred(A[1]) == ref
            ref = (datasubset(var,(28*5+1):(29*5)), datasubset(var,29*5+1))
            @test @inferred(A[29]) == ref
            @test A[end] == A[29]
            ref = (getobs(var,1:5), getobs(var,6))
            @test @inferred(getobs(A,1)) == ref
            @test typeof(getobs(A,1)) == typeof(ref)
            ref = [
                (getobs(var,(1*5+1):(2*5)), getobs(var,2*5+1)),
                (getobs(var,(2*5+1):(3*5)), getobs(var,3*5+1))
            ]
            @test @inferred(getobs(A,2:3)) == ref
            @test typeof(getobs(A,2:3)) == typeof(ref)
            @test typeof(@inferred(collect(A))) <: Vector
            # custom stride
            A = slidingwindow(i->i-1, var, 5, 2)
            @test_throws BoundsError A[-1]
            @test_throws BoundsError A[length(A)+1]
            @test @inferred(nobs(A)) == 72
            @test @inferred(length(A)) == 72
            @test @inferred(size(A)) == (72,)
            @test A[end] == A[72]
            ref = (getobs(var,3:3+4), getobs(var,2))
            @test @inferred(getobs(A,1)) == ref
            ref = [
                (getobs(var,5:5+4), getobs(var,4)),
                (getobs(var,7:7+4), getobs(var,6)),
            ]
            @test @inferred(getobs(A,2:3)) == ref
            @test typeof(@inferred(collect(A))) <: Vector
            # TODO: exclude target
        end
    end

    println("<HEARTBEAT>")

    @testset "batchview" begin
        # number as target
        sw = slidingwindow(i->i+5, 1:31, 5)
        @test length(sw) == 6
        bv = @inferred batchview(sw, 2)
        @test_throws MethodError shuffleobs(bv)
        @test_throws MethodError splitobs(bv)
        @test_throws MethodError oversample(bv)
        @test_throws MethodError undersample(bv)
        @test_throws MethodError stratifiedobs(bv)
        @test length(bv) == 3
        @test nobs(bv) == 6
        tx, ty = @inferred bv[1]
        @test tx isa Array{<:SubArray}
        @test ty isa Array{<:SubArray}
        @test tx == [[1,6],[2,7],[3,8],[4,9],[5,10]]
        @test ty == [[6,11]]
        tx, ty = @inferred bv[2]
        @test tx == [[11,16],[12,17],[13,18],[14,19],[15,20]]
        @test ty == [[16,21]]
        tx, ty = @inferred bv[3]
        @test tx == [[21,26],[22,27],[23,28],[24,29],[25,30]]
        @test ty == [[26,31]]
        # also test with array as target
        sw = slidingwindow(i->i.+(5:6), 1:32, 5)
        @test length(sw) == 6
        bv = @inferred batchview(sw, 2)
        @test_throws MethodError shuffleobs(bv)
        @test_throws MethodError splitobs(bv)
        @test_throws MethodError oversample(bv)
        @test_throws MethodError undersample(bv)
        @test_throws MethodError stratifiedobs(bv)
        @test length(bv) == 3
        @test nobs(bv) == 6
        tx, ty = @inferred bv[1]
        @test tx isa Array{<:SubArray}
        @test ty isa Array{<:SubArray}
        @test tx == [[1,6],[2,7],[3,8],[4,9],[5,10]]
        @test ty == [[6,11],[7,12]]
        tx, ty = @inferred bv[2]
        @test tx == [[11,16],[12,17],[13,18],[14,19],[15,20]]
        @test ty == [[16,21],[17,22]]
        tx, ty = @inferred bv[3]
        @test tx == [[21,26],[22,27],[23,28],[24,29],[25,30]]
        @test ty == [[26,31],[27,32]]
    end

    println("<HEARTBEAT>")

    @testset "special values" begin
        A = @inferred slidingwindow(i->i+2, 1:10, 2)
        @test length(A) == 4
        @test getobs(A) == [
            ([1,2], 3),
            ([3,4], 5),
            ([5,6], 7),
            ([7,8], 9),
        ]
        A = @inferred slidingwindow(i->i+2, 1:10, 2, 3)
        @test length(A) == 3
        @test getobs(A) == [
            ([1,2], 3),
            ([4,5], 6),
            ([7,8], 9),
        ]
        A = @inferred slidingwindow(i->i+2:i+3, 1:10, 2)
        @test length(A) == 4
        @test getobs(A) == [
            ([1,2], [3,4]),
            ([3,4], [5,6]),
            ([5,6], [7,8]),
            ([7,8], [9,10]),
        ]
        A = @inferred slidingwindow(i->[i-2:i-1; i+1:i+2], 1:10, 1)
        @test length(A) == 6
        @test getobs(A) == [
            ([3], [1, 2, 4, 5]),
            ([4], [2, 3, 5, 6]),
            ([5], [3, 4, 6, 7]),
            ([6], [4, 5, 7, 8]),
            ([7], [5, 6, 8, 9]),
            ([8], [6, 7, 9, 10]),
        ]
        A = @inferred slidingwindow(i->i-1, 1:10, 2)
        @test length(A) == 4
        @test getobs(A) == [
            ([3,4], 2),
            ([5,6], 4),
            ([7,8], 6),
            ([9,10], 8),
        ]
        A = @inferred slidingwindow(i->i+1, 1:10, 3)
        @test length(A) == 3
        @test getobs(A) == [
            ([1,2,3], 2),
            ([4,5,6], 5),
            ([7,8,9], 8),
        ]
        A = @inferred slidingwindow(i->i+1, 1:10, 3, Val(true))
        @test length(A) == 3
        @test getobs(A) == [
            ([1,3], 2),
            ([4,6], 5),
            ([7,9], 8),
        ]
    end
end


text:
using GRAMPC
using Test
using SafeTestsets

@testset "GRAMPC.jl" begin
    @testset "C Interface" begin
        @safetestset "Error handling" begin include( "cinterface/test_error_handling.jl" ) end
    end
end

ext: jl


text:
#####
##### Callback placeholders
#####

"""
    Always{T}(val)

Functor that allways returns `val` when called.
"""
struct Always{T}
    val::T
end
(f::Always)(x...; kw...) = f.val

const donothing = Always(nothing)
const alwaysfalse = Always(false)

macro _ivdep_meta()
    return Expr(:loopinfo, Symbol("julia.ivdep"))
end

#####
##### printing
#####

# Print your message with out worrying about task switching.
# Suitable for running from withing finalizers.
"""
    safeprint(message; [force])

Print `message` in a way that is safe to do inside finalizers (can't trigger a task switch).
If keyword `force = false`, than message will only printed if `VERBOSE = true`.
"""
function safeprint(message; force = false)
    if (force || VERBOSE)
        ccall(:jl_safe_printf, Cvoid, (Cstring,), "$message\n")
    end
    return nothing
end

#####
##### internal macros
#####

cpu_pause() = ccall(:jl_cpu_pause, Cvoid, ())
safepoint() = ccall(:jl_gc_safepoint, Cvoid, ())

function _spinlock_impl(expr, canfail::Bool, locks)
    # If the inner expression can fail, then we need to wrap inner calls in a try-finally
    # block. If it is stated that `expr` cannot fail, then we can get some slightly better
    # code generation by removing these try-finally blocks.
    function inner(x, lock)
        if canfail
            return quote
                try
                    $x
                finally
                    unlock($lock)
                end
            end
        else
            return quote
                $x
                unlock($lock)
            end
        end
    end
    expr = esc(expr)
    locks = map(esc, locks)

    # Build expression from the inner most to the outermost to ensure locks are acquired
    # in the order given at the macro invocation.
    for lock in reverse(locks)
        v = gensym("lock")
        expr = quote
            $v = $lock
            while !trylock($v)
                cpu_pause()
                safepoint()
            end
            $(inner(expr, v))
        end
    end
    return expr
end

# a safe way to acquire locks from finalizers, where we can't wait (which switches tasks)
# Copied from CUDA.jl and Julia's "locks-mt.jl"
"""
    @spinlock lock1 lock2 ... lockN expr

Acquire locks in order than execute `expr`.
Acquiring the locks will not trigger a task switch.
"""
macro spinlock(l...)
    return _spinlock_impl(l[end], true, l[1:end-1])
end

"""
    @spinlock_nofail lock1 lock2 ... lockN expr

Acquire `locks` in order than execute `expr`.
Acquiring the locks will not trigger a task switch.
If `expr` throws an error, locks will not be released.
"""
macro spinlock_nofail(l...)
    return _spinlock_impl(l[end], false, l[1:end-1])
end

"""
    @requires lock1 lock2 ...

Assert that all locks are held. Only valid in debug mode.
"""
macro requires(locks...)
    lock_exprs = [:(Base.assert_havelock($(esc(lock)))) for lock in locks]

    return quote
        if DEBUG
            $(lock_exprs...)
        end
    end
end

macro checknothing(expr, action = :(return nothing))
    return quote
        x = $(esc(expr))
        x === nothing && $(esc(action))
        x
    end
end

"""
    @return_if_exists expr

If the result of `expr` is not `nothing`, then return that result.
Otherwise, continue execution.
"""
macro return_if_exists(expr)
    return quote
        x = $(esc(expr))
        x === nothing || return x
        x
    end
end

#####
##### Broadcast Destructuring
#####

# We need to search through the `bc` object for the first instance of T
# Once we find it, return it.
findT(::Type{T}, bc::Base.Broadcast.Broadcasted) where {T} = findT(T, bc.args)
findT(::Type{T}, x::Tuple) where {T} = findT(T, findT(T, first(x)), Base.tail(x))
findT(::Type{T}, x::U) where {T,U<:T} = x
findT(::Type{T}, x::SubArray{<:Any,<:Any,U}) where {T<:AbstractArray,U<:T} = findT(T, parent(x))
findT(::Type{T}, x) where {T} = nothing
findT(::Type{T}, ::Tuple{}) where {T} = nothing
findT(::Type{T}, x::U, rest) where {T,U<:T} = x
findT(::Type{T}, ::Any, rest) where {T} = findT(T, rest)

#####
##### ObjectCache
#####

mutable struct ObjectCache{T,Args}
    cache::Vector{T}
    head::Int
    args::Args
end

ObjectCache{T}() where {T} = ObjectCache{T,Tuple{}}(Vector{T}(), 0, ())

function Base.getindex(cache::ObjectCache{T}) where {T}
    head = cache.head
    if head > 0
        object = @inbounds(cache.cache[head])
        cache.head = head - 1
        return object
    end
    return T(cache.args...)
end

function return!(cache::ObjectCache{T}, x::T) where {T}
    head = cache.head
    if head < length(cache.cache)
        head += 1
        cache.cache[head] = x
    else
        push!(cache.cache, x)
    end
    return nothing
end

#####
##### Exceptions
#####

struct AllocationException <: Exception end

"""
$(TYPEDSIGNATURES)

Atomically store `v` to the location pointed to by `ptr` and return the old values
stored at `ptr`.
"""
function atomic_ptr_xchg!(ptr::Ptr{Ptr{T}}, v::Ptr{T}) where {T}
    return Base.llvmcall("""
        %ptr = inttoptr i64 %0 to i64*
        %rv = atomicrmw xchg i64* %ptr, i64 %1 acq_rel
        ret i64 %rv
        """,
        Ptr{T},
        Tuple{Ptr{Ptr{T}},Ptr{T}},
        ptr,
        v,
    )
end

#####
##### ToCached
#####

# Convenience Constructors
struct ToCached{M,S}
    manager::M
    status::S
    priority::CachedArrays.AllocationPriority
end

function tocached(m::M, s::S = CachedArrays.NotBusy(); priority = CachedArrays.ForceLocal) where {M,S}
    return ToCached{M,S}(m, s, priority)
end

function (f::ToCached)(::Type{T}, x...) where {T}
    return CachedArray{T}(undef, f.manager, x; f.status, f.priority)
end


text:
module TestJacobian
include("preamble.jl")

using DiffEqBase: ODEProblem
using StaticArrays: SVector
using SymEngine: diff, subs, symbols

using Bifurcations.Continuations: residual, residual_jacobian, get_prob_cache
using Bifurcations.Codim2: DiffEqCodim2Problem, NormalizingAS
using Bifurcations.Examples.PredatorPrey

∙(a, b) = sum(a .* b)
sjac(f, x) = [diff(f[i], x[j]) for i in 1:length(f), j in 1:length(x)]

nx = length(PredatorPrey.u0)
np = length(PredatorPrey.p)
f = generalized_f(PredatorPrey)
pe = [symbols("p_$i") for i in 1:np]
xe = [symbols("x_$i") for i in 1:nx]
ve = [symbols("v_$i") for i in 1:nx]
all_symbols = vcat(pe, xe, ve)
fe = f(xe, pe, 0)
fe = Vector(fe)  # just for safety
dfdx = sjac(fe, xe)

ue = vcat(xe, ve, pe[1], pe[3])
He = vcat(fe,
          dfdx * ve,
          ve ∙ ve - 1)
dHdu = sjac(He, ue)

make_prob(;
        p = (0.0, 3, 5, 3),
        u0 = SVector(0.0, 0.0),
        tspan = (0.0, 30.0),
        ode = ODEProblem{!(u0 isa SVector)}(PredatorPrey.f, u0, tspan, p),
        kwargs...) =
    DiffEqCodim2Problem(
        ode,
        (@lens _[1]),
        (@lens _[3]),
        ([-Inf, -Inf], [Inf, Inf]);
        augmented_system = NormalizingAS(),
        kwargs...)

function test_predator_prey(
        xr, vr, p1, p3;
        u0 = SVector(0.0, 0.0),
        prob = make_prob(u0 = u0),
        prob_cache = get_prob_cache(prob))

    pr = collect(Float64, prob.de_prob.p)
    pr[1] = p1
    pr[3] = p3

    ur = vcat(xr, vr, pr[1], pr[3])
    if prob.de_prob.u0 isa SVector
        ur = SVector(ur...)
    else
        utype = typeof(prob.de_prob.u0)
        if !(ur isa utype)
            ur = utype(ur)
        end
    end

    all_values = vcat(pr, xr, vr) :: Vector{Float64}
    evaluated = (e) -> Float64(subs(e, zip(all_symbols, all_values)...))

    H_desired = evaluated.(He)
    J_desired = evaluated.(dHdu)

    H1_actual = residual(ur, prob_cache)
    @test H1_actual ≈ H_desired

    H2_actual, J_actual = residual_jacobian(ur, prob_cache)
    @test H2_actual ≈ H_desired
    @test J_actual ≈ J_desired
end

rng = MersenneTwister(0)
@testset begin
    for args in [
            ([0, 0], [0, 1], 0.6, 5),
            ([0, 0], [-1, 0], 0.384806, 7.79613),
            ],
        u0 in [SVector(0.0, 0.0), [0.0, 0.0]]

        test_predator_prey(args...; u0=u0)
    end
    for _ in 1:100,
        u0 in [SVector(0.0, 0.0), [0.0, 0.0]]

        p1 = randn(rng)
        p3 = randn(rng)
        xr = randn(rng, nx)
        vr = normalize(randn(rng, nx))
        test_predator_prey(xr, vr, p1, p3; u0=u0)
    end
end

end  # module


text:
"""
    FeatureDescriptors.TestUtils

Provides a fake [`Descriptor`](@ref) for testing purposes and a `test_interface` function
to check that any [`Descriptor`](@ref) subtype implements the expected API.
"""
module TestUtils

using ..FeatureDescriptors
using Test

export FakeDescriptor, test_interface

"""
    FakeDescriptor <: Descriptor

A fake [`Descriptor`](@ref) for testing purposes only.
"""
abstract type FakeDescriptor <: Descriptor end

FeatureDescriptors.sources(::Type{<:FakeDescriptor}) = ["fake_table"]
FeatureDescriptors.quantity_key(::Type{<:FakeDescriptor}) = :quantity
FeatureDescriptors.categorical_keys(::Type{<:FakeDescriptor}) = [:category1, :category2]

"""
    TestUtils.test_interface(D::Type{<:Descriptor})

Test that a subtype of [`Descriptor`](@ref) implements the expected API.
"""
function test_interface(D)
    @testset "test_interface $D" begin
        @test FeatureDescriptors.sources(D) isa Vector{String}
        @test FeatureDescriptors.quantity_key(D) isa Symbol
        categories = FeatureDescriptors.categorical_keys(D)
        @test (isempty(categories) || categories isa Vector{Symbol})
        @test FeatureDescriptors.label(D) isa Symbol
    end
end

end

max_stars_repo_name: invenia/FeatureDescriptors

max_issues_repo_name: invenia/FeatureDescriptors

max_forks_repo_name: invenia/FeatureDescriptors


text:
@testset "broadcasted" begin
    @test gradcheck(x -> sum(tanh.(x)), rand(2))
end


text:
"""

    ccf(setting; starting_lambdas = nothing)

Perform signed gradient descent for clipped convex functions for a given regression setting.

# Arguments
- `setting::RegressionSetting`: RegressionSetting object with a formula and dataset.
- `starting_lambdas::Array{Float64,1}`: Starting values of weighting parameters used by signed gradient descent.
- `alpha::Float64`: Loss at which a point is labeled as an outlier (points with loss ≥ alpha will be called outliers).
- `max_iter::Int64`: Maximum number of iterations to run signed gradient descent.
- `beta::Float64`: Step size parameter.
- `tol::Float64`: Tolerance below which convergence is declared.

# Examples
```julia-repl
julia> reg0001 = createRegressionSetting(@formula(calls ~ year), phones);
julia> ccf(reg0001)
Dict{Any,Any} with 4 entries:
  "betas"     => [-63.4816, 1.30406]
  "outliers"  => [15, 16, 17, 18, 19, 20]
  "lambdas"   => [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  2.77556e-17, 2.77556e-17, 0…
  "residuals" => [-2.67878, -1.67473, -0.37067, -0.266613, 0.337444, 0.941501, 1.44556, 2.04962, 1…

```

# References
Barratt, S., Angeris, G. & Boyd, S. Minimizing a sum of clipped convex functions. Optim Lett 14, 2443–2459 (2020). https://doi.org/10.1007/s11590-020-01565-4

"""
function ccf(setting::RegressionSetting; starting_lambdas=nothing, alpha=nothing, p=3, max_iter=100, gamma=.1, tol=1e-4)
    X, y = @extractRegressionSetting setting
    return ccf(X, y, starting_lambdas=starting_lambdas, alpha=alpha, p=p, max_iter=max_iter, gamma=gamma, tol=tol)
end

"""

    ccf(X, y; starting_lambdas = nothing)

Perform signed gradient descent for clipped convex functions for a given regression setting.

# Arguments
- `X::Array{Float64, 2}`: Design matrix of the linear model.
- `y::Array{Float64, 1}`: Response vector of the linear model.
- `starting_lambdas::Array{Float64,1}`: Starting values of weighting parameters used by signed gradient descent.
- `alpha::Float64`: Loss at which a point is labeled as an outlier. If unspecified, will be chosen as p*mean(residuals.^2), where residuals are OLS residuals.
- `p::Float64`: Points that have squared OLS residual greater than p times the mean squared OLS residual are considered outliers.
- `max_iter::Int64`: Maximum number of iterations to run signed gradient descent.
- `beta::Float64`: Step size parameter.
- `tol::Float64`: Tolerance below which convergence is declared.

# References
Barratt, S., Angeris, G. & Boyd, S. Minimizing a sum of clipped convex functions. Optim Lett 14, 2443–2459 (2020). https://doi.org/10.1007/s11590-020-01565-4

"""
function ccf(X::Array{Float64,2}, y::Array{Float64,1}; starting_lambdas=nothing, alpha=nothing, p=3, max_iter=100, gamma=.1, tol=1e-4)
    n, p = size(X)

    if isnothing(starting_lambdas)
        starting_lambdas = ones(Float64, n) / 2
    end

    curr_lambdas = copy(starting_lambdas)
    old_lambdas = copy(starting_lambdas)
    curr_betas = nothing
    residuals = nothing

    for iter = 1:max_iter
        curr_betas = wls(X, y, curr_lambdas).betas
        residuals = X * curr_betas - y

        if isnothing(alpha)
            alpha = p * mean(residuals.^2)
        end

        @. curr_lambdas -= gamma * sign(residuals^2 - alpha)
        clamp!(curr_lambdas, 0., 1.)

        if norm(curr_lambdas - old_lambdas, Inf) <= tol
            break
        end

        old_lambdas .= curr_lambdas
    end

    result = Dict()
    result["betas"] = curr_betas
    result["lambdas"] = old_lambdas
    result["residuals"] = residuals
    result["outliers"] = findall(residuals.^2 .> alpha)

    return result
end

max_stars_repo_name: angeris/LinRegOutliers

max_issues_repo_name: angeris/LinRegOutliers

max_forks_repo_name: angeris/LinRegOutliers


text:
export Chain

"""
    Chain(layers...)

Chain multiple layers / functions together, so that they are called in sequence
on a given input.

```julia
m = Chain(x -> x^2, x -> x+1)
m(5) == 26
m = Chain(Dense(10, 5), Dense(5, 2))
x = rand(10)
m(x) == m[2](m[1](x))
```

`Chain` also supports indexing and slicing, e.g. `m[2]` or `m[1:end-1]`.
`m[1:3](x)` will calculate the output of the first three layers.

# Note
Chain is simply a Neural Network, but it is not a quantum state. To use it with
NeuralQuantum you should wrap it into a `PureStateAnsatz`.
"""
struct Chain{T<:Tuple}
  layers::T
  Chain(xs...) = new{typeof(xs)}(xs)
end

@forward Chain.layers Base.getindex, Base.length, Base.first, Base.last,
  Base.iterate, Base.lastindex

functor(c::Chain) = c.layers, ls -> Chain(ls...)

applychain(::Tuple{}, x) = x
applychain(fs::Tuple, x) = applychain(Base.tail(fs), first(fs)(x))

(c::Chain)(x) = applychain(c.layers, x)

Base.getindex(c::Chain, i::AbstractArray) = Chain(c.layers[i]...)

function Base.show(io::IO, c::Chain)
  print(io, "Chain(")
  join(io, c.layers, ", ")
  print(io, ")")
end

struct ChainCache{T<:Tuple} <: NNCache{Chain}
    caches::T
    valid::Bool
end

cache(l::Chain, arr_T, in_T, in_sz) = begin
    caches = []
    for layer = l.layers
        c = cache(layer, arr_T, in_T, in_sz)
        in_T, in_sz = layer_out_type_size(layer, in_T, in_sz)
        push!(caches, c)
    end
    ChainCache(Tuple(caches), false)
end

# Applychain with caches
applychain(::Tuple{}, ::Tuple{}, x) = x
applychain(fs::Tuple, fc::Tuple, x) = applychain(Base.tail(fs),
                                                 Base.tail(fc),
                                                 first(fs)(first(fc), x))
(c::Chain)(ch::ChainCache, x) = applychain(c.layers, ch.caches, x)

backpropchain(∇, ::Tuple{}, ::Tuple{}, δℒ) = δℒ
function backpropchain(∇, fs::Tuple, fc::Tuple, δℒ)
    a = reverse(Base.tail(reverse(∇)))
    b = reverse(Base.tail(reverse(fs)))
    c = reverse(Base.tail(reverse(fc)))

    backpropchain(a,b,c,
        backprop(last(∇), last(fs), last(fc), δℒ))
end

function logψ_and_∇logψ!(∇lnψ, net::Chain, c::ChainCache, σ)
    # forward pass
    lnψ = net(c, σ)

    # backward
    backpropchain(∇lnψ, net.layers, c.caches, 1.0)
    return lnψ
end

is_analytic(net::Chain) = true
out_type(net::Chain) = ComplexF64


text:
"""
    @gql_str(document, throw_on_error=true)

Create and optionally validate a GraphQL query string.

The string is parsed and semi-validated by GraphQLParser.jl.
Validation that does not need the schema from the server is performed.
For further information see the GraphQLParser documentation.

Parsing errors will always be thrown, but other validation errors can be turned off using the second argument.

An additional advantage of using this macro is that dollar signs do not need to be escaped (see example below).

# Examples

General usage

```julia-repl
julia> client = Client("https://countries.trevorblades.com");

julia> str = gql\"""
query(\$code: ID!){
    country(
        code:\$code
    ){
        name
    }
}
\""";

julia> variables = Dict("code" => "BR");

julia> GraphQLClient.execute(client, str; variables=variables)
GraphQLClient.GQLResponse{Any}
  data: Dict{String, Any}
          country: Dict{String, Any}
```

Parsing error

```julia-repl
julia> str = gql\"""
query(code: ID!){  # no \$ before variable name
    country(
        code:\$code
    ){
        name
    }
}
\""";

# ERROR: LoadError: ArgumentError: invalid GraphQL string at byte position 7 while parsing
#     Variable name must start with '\$'
#     query(code: ID!){  # no \$ before
#           ^
```

Validation error

```julia-repl
julia> str = gql\"""
{
    countries{
        name
    }
}

query{  # Can't have another operation when there is an anonymous operation
    countries{
        name
    }
}
\""";
# ERROR: LoadError: Validation Failed

# GraphQLParser.AnonymousOperationNotAlone
#       message: This anonymous operation must be the only defined operation.
#      location: Line 1 Column 1
```

Turning validation off

```julia-repl
julia> str = @gql_str \"""
{
    countries{
        name
    }
}

query{  # Can't have another operation when there is an anonymous operation
    countries{
        name
    }
}
\""" false
# No error
```
"""
macro gql_str(document, throw_on_error=true)
    is_valid_executable_document(document; throw_on_error)
    return document
end


text:
module LTVSourceReceiverModel

include("omniDirectionalLTISource/source.jl")
export LTISource

include("omniDirectionalLTISource/listener.jl")
export LTIListener

end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    BallSampleValidation(eestimator, sourceball, targetball;
                         tol=1e-4, maxiter=10)

Ball sample validation with error estimator `eestimator`, based on
samples collected from source and target data. `sourceball` is the
ball for the source and `targetball` is the ball for the target.

    BallSampleValidation(eestimator, sourceradius, targetradius;
                         tol=1e-4, maxiter=10)

Alternatively, specify the radii for Euclidean balls.
"""
struct BallSampleValidation{E<:AbstractErrorEstimator,
                            Bₛ<:BallNeighborhood,
                            Bₜ<:BallNeighborhood} <: AbstractErrorEstimator
  eestimator::E
  sourceball::Bₛ
  targetball::Bₜ
  tol::Float64
  maxiter::Int
end

BallSampleValidation(eestimator::AbstractErrorEstimator,
                     sball::BallNeighborhood,
                     tball::BallNeighborhood;
                     tol=1e-4, maxiter=10) =
  BallSampleValidation(eestimator, sball, tball, tol, maxiter)

BallSampleValidation(eestimator::AbstractErrorEstimator,
                     sradius::Real, tradius::Real;
                     tol=1e-4, maxiter=10) =
  BallSampleValidation(eestimator,
                       BallNeighborhood(sradius),
                       BallNeighborhood(tradius);
                       tol=tol, maxiter=maxiter)

function Base.error(solver::AbstractLearningSolver,
                    problem::LearningProblem,
                    eestimator::BallSampleValidation)
  # source and target data
  sdata = sourcedata(problem)
  tdata = targetdata(problem)
  ovars = outputvars(task(problem))

  # source and target balls
  sball = eestimator.sourceball
  tball = eestimator.targetball

  # source and target samplers
  ssampler = BallSampler(sball)
  tsampler = BallSampler(tball)

  # helper function
  function iid_estimate(ssamples, tsamples)
    subproblem = LearningProblem(ssamples, tsamples, task(problem))
    error(solver, subproblem, eestimator.eestimator)
  end

  # initialization with i.i.d. samples
  res = iid_estimate(sample(sdata, ssampler), sample(tdata, tsampler))

  # iteration parameters
  tol     = eestimator.tol
  maxiter = eestimator.maxiter
  err     = Dict(var => Inf for var in ovars)
  iter    = 0
  while maximum(values(err)) > tol && iter < maxiter
    # copy old result
    old = copy(res)
    iter += 1

    # estimate with more i.i.d. samples
    new = iid_estimate(sample(sdata, ssampler), sample(tdata, tsampler))

    # update result (online mean formula)
    for var in ovars
      res[var] = (iter*old[var] + new[var])/(iter + 1)
      err[var] = abs(res[var] - old[var]) / old[var]
    end
  end

  res
end


text:
import TriangleMesh, UnicodePlots, Plots
using Printf, LoopVectorization, LinearAlgebra, SparseArrays

include("CreateMeshFCFV.jl")
include("VisuFCFV.jl")
include("DiscretisationFCFV_Stokes.jl")
include("SolversFCFV_Stokes.jl")

#--------------------------------------------------------------------#

function SetUpProblem!(mesh, P, Vx, Vy, Sxx, Syy, Sxy, VxDir, VyDir, SxxNeu, SyyNeu, SxyNeu, SyxNeu, sx, sy)
    # Evaluate T analytic on cell faces
    @avx for in=1:mesh.nf
        eta      = 1
        x        = mesh.xf[in]
        y        = mesh.yf[in]
        VxDir[in] = x^2*(1 - x)^2*(4*y^3 - 6*y^2 + 2*y)
        VyDir[in] =-y^2*(1 - y)^2*(4*x^3 - 6*x^2 + 2*x)
        # Stress at faces
        p         =  x*(1-x)
        # Using pseudo-tractions
        dVxdx = x^2*(2*x - 2)*(4*y^3 - 6*y^2 + 2*y) + 2*x*(1 - x)^2*(4*y^3 - 6*y^2 + 2*y);
        dVxdy = x^2*(1 - x)^2*(12*y^2 - 12*y + 2);
        dVydx = -y^2*(1 - y)^2*(12*x^2 - 12*x + 2);
        dVydy = -y^2*(2*y - 2)*(4*x^3 - 6*x^2 + 2*x) - 2*y*(1 - y)^2*(4*x^3 - 6*x^2 + 2*x);
        SxxNeu[in] = - p + eta*dVxdx
        SyyNeu[in] = - p + eta*dVydy
        SxyNeu[in] =       eta*dVxdy
        SyxNeu[in] =       eta*dVydx
        # SxxNeu[in] = -8*p*y*(x - 1)*(2*y^2 - 3*y + 1) - p + 8*x^2*y*(x - 1)*(2*y^2 - 3*y + 1)
        # SyyNeu[in] = -p - 8*x*y^2*(y - 1)*(2*x^2 - 3*x + 1) - 8*x*y*(y - 1)^2*(2*x^2 - 3*x + 1)
        # SxyNeu[in] = p^2*(12*y^2 - 12*y + 2) + y^2*(y - 1)^2*(-12.0*x^2 + 12.0*x - 2.0)
        # SyxNeu[in] = SxyNeu[in]
    end
    # Evaluate T analytic on barycentres
    @avx for iel=1:mesh.nel
        x        = mesh.xc[iel]
        y        = mesh.yc[iel]
        p        =  x*(1-x)
        P[iel]   =  p
        Vx[iel]  =  x^2*(1 - x)^2*(4*y^3 - 6*y^2 + 2*y)
        Vy[iel]  = -y^2*(1 - y)^2*(4*x^3 - 6*x^2 + 2*x)
        Sxx[iel] = -8*p*y*(x - 1)*(2*y^2 - 3*y + 1) - p + 8*x^2*y*(x - 1)*(2*y^2 - 3*y + 1)
        Syy[iel] = -p - 8*x*y^2*(y - 1)*(2*x^2 - 3*x + 1) - 8*x*y*(y - 1)^2*(2*x^2 - 3*x + 1)
        Sxy[iel] = p^2*(12*y^2 - 12*y + 2) + y^2*(y - 1)^2*(-12.0*x^2 + 12.0*x - 2.0)
        sx[iel]  = -p^2*(24*y - 12) - 4*x^2*(4*y^3 - 6*y^2 + 2*y) - 8*x*(2*x - 2)*(4*y^3 - 6*y^2 + 2*y) - 2*x + 1.0*y^2*(2*y - 2)*(12*x^2 - 12*x + 2) + 2.0*y*(1 - y)^2*(12*x^2 - 12*x + 2) - 4*(1 - x)^2*(4*y^3 - 6*y^2 + 2*y) + 1
        sy[iel]  = -2*p*(1 - x)*(12*y^2 - 12*y + 2) - x^2*(2*x - 2)*(12*y^2 - 12*y + 2) + 1.0*y^2*(1 - y)^2*(24*x - 12) + 4*y^2*(4*x^3 - 6*x^2 + 2*x) + 8*y*(2*y - 2)*(4*x^3 - 6*x^2 + 2*x) + 4*(1 - y)^2*(4*x^3 - 6*x^2 + 2*x)
    end
    return
end

#--------------------------------------------------------------------#

function ComputeError( mesh, Vxe, Vye, Txxe, Tyye, Txye, Pe )
    eVx  = zeros(mesh.nel)
    eVy  = zeros(mesh.nel)
    eTxx = zeros(mesh.nel)
    eTyy = zeros(mesh.nel)
    eTxy = zeros(mesh.nel)
    eP   = zeros(mesh.nel)
    Vxa  = zeros(mesh.nel)
    Vya  = zeros(mesh.nel)
    Txxa = zeros(mesh.nel)
    Tyya = zeros(mesh.nel)
    Txya = zeros(mesh.nel)
    Pa   = zeros(mesh.nel)
    @avx for iel=1:mesh.nel
        x         = mesh.xc[iel]
        y         = mesh.yc[iel]
        Pa[iel]   =  x*(1-x)
        Vxa[iel]  =  x^2*(1 - x)^2*(4*y^3 - 6*y^2 + 2*y)
        Vya[iel]  = -y^2*(1 - y)^2*(4*x^3 - 6*x^2 + 2*x)
        Txxa[iel] = -8*Pa[iel]*y*(x - 1)*(2*y^2 - 3*y + 1) + 8*x^2*y*(x - 1)*(2*y^2 - 3*y + 1)
        Tyya[iel] = -8*x*y^2*(y - 1)*(2*x^2 - 3*x + 1) - 8*x*y*(y - 1)^2*(2*x^2 - 3*x + 1)
        Txya[iel] = Pa[iel]^2*(12*y^2 - 12*y + 2) + y^2*(y - 1)^2*(-12.0*x^2 + 12.0*x - 2.0)
        eVx[iel]  = Vxe[iel]  - Vxa[iel]
        eVy[iel]  = Vye[iel]  - Vya[iel]
        eTxx[iel] = Txxe[iel] - Txxa[iel]
        eTyy[iel] = Tyye[iel] - Tyya[iel]
        eTxy[iel] = Txye[iel] - Txya[iel]
        eP[iel]   = Pe[iel]   - Pa[iel]
    end
    errVx  = norm(eVx) /norm(Vxa)
    errVy  = norm(eVy) /norm(Vya)
    errTxx = norm(eTxx)/norm(Txxa)
    errTyy = norm(eTyy)/norm(Tyya)
    errTxy = norm(eTxy)/norm(Txya)
    errP   = norm(eP)  /norm(Pa)
    return errVx, errVy, errTxx, errTyy, errTxy, errP, Txxa, Tyya, Txya
end

#--------------------------------------------------------------------#

function StabParam(tau, dA, Vol, mesh_type, coeff)
    if mesh_type=="Quadrangles";        taui = tau;    end
    # if mesh_type=="UnstructTriangles";  taui = tau*dA; end
    if mesh_type=="UnstructTriangles";  taui = tau end
    return taui
end

#--------------------------------------------------------------------#

@views function main()

    println("\n******** FCFV STOKES ********")

    # Create sides of mesh
    xmin, xmax = 0, 1
    ymin, ymax = 0, 1
    n          = 4
    nx, ny     = 8*n, 8*n
    solver     = 1
    R          = 0.5
    inclusion  = 0
    # mesh_type  = "Quadrangles"
    mesh_type  = "UnstructTriangles"

    # Generate mesh
    if mesh_type=="Quadrangles"
        tau  = 20
        mesh = MakeQuadMesh( nx, ny, xmin, xmax, ymin, ymax, inclusion, R )
    elseif mesh_type=="UnstructTriangles"
        tau  = 20
        mesh = MakeTriangleMesh( nx, ny, xmin, xmax, ymin, ymax, inclusion, R )
    end
    println("Number of elements: ", mesh.nel)

    # Source term and BCs etc...
    Pa     = zeros(mesh.nel)
    Vxa    = zeros(mesh.nel)
    Vya    = zeros(mesh.nel)
    Sxxa   = zeros(mesh.nel)
    Syya   = zeros(mesh.nel)
    Sxya   = zeros(mesh.nel)
    sex    = zeros(mesh.nel)
    sey    = zeros(mesh.nel)
    VxDir  = zeros(mesh.nf)
    VyDir  = zeros(mesh.nf)
    SxxNeu = zeros(mesh.nf)
    SyyNeu = zeros(mesh.nf)
    SxyNeu = zeros(mesh.nf)
    SyxNeu = zeros(mesh.nf)
    gbar   = zeros(mesh.nel,mesh.nf_el,2)
    println("Model configuration :")
    @time SetUpProblem!(mesh, Pa, Vxa, Vya, Sxxa, Syya, Sxya, VxDir, VyDir, SxxNeu, SyyNeu, SxyNeu, SyxNeu, sex, sey)

    # Compute some mesh vectors
    println("Compute FCFV vectors:")
    @time ae, be, ze = ComputeFCFV(mesh, sex, sey, VxDir, VyDir, SxxNeu, SyyNeu, SxyNeu, SyxNeu, tau)

    # Assemble element matrices and RHS
    println("Compute element matrices:")
    @time Kuu_v, fu_v, Kup_v, fp = ElementAssemblyLoop(mesh, ae, be, ze, VxDir, VyDir, SxxNeu, SyyNeu, SxyNeu, SyxNeu, gbar, tau)

    # Assemble triplets and sparse
    println("Assemble triplets and sparse:")
    @time Kuu, fu, Kup = CreateTripletsSparse(mesh, Kuu_v, fu_v, Kup_v)

    # Solve for hybrid variable
    println("Linear solve:")
    @time Vxh, Vyh, Pe = StokesSolvers(mesh, Kuu, Kup, fu, fp, solver)

    # # Reconstruct element values
    println("Compute element values:")
    @time Vxe, Vye, Txxe, Tyye, Txye = ComputeElementValues(mesh, Vxh, Vyh, Pe, ae, be, ze, VxDir, VyDir, tau)

    # # Compute discretisation errors
    err_Vx, err_Vy, err_Txx, err_Tyy, err_Txy, err_P, Txxa, Tyya, Txya = ComputeError( mesh, Vxe, Vye, Txxe, Tyye, Txye, Pe )
    @printf("Error in Vx : %2.2e\n", err_Vx )
    @printf("Error in Vy : %2.2e\n", err_Vy )
    @printf("Error in Txx: %2.2e\n", err_Txx)
    @printf("Error in Tyy: %2.2e\n", err_Tyy)
    @printf("Error in Txy: %2.2e\n", err_Txy)
    @printf("Error in P  : %2.2e\n", err_P  )

    println(minimum(Pa))
    println(minimum(Pe))
    println(maximum(Pa))
    println(maximum(Pe))

    # Visualise
    println("Visualisation:")
    # @time PlotMakie( mesh, sex )
    @time PlotMakie( mesh, Pe, xmin, xmax, ymin, ymax, :viridis )
    # PlotElements( mesh )

end

main()


text:
mutable struct MyCallable
    f
end
(a::MyCallable)(x) = a.f(x)
@testset "callables" begin
    @testset "univariate" begin
    # FIXME
    end
    @testset "multivariate" begin
        function rosenbrock(x::Vector)
           return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
        end

        a = MyCallable(rosenbrock)

        optimize(a, rand(2), Optim.Options())
    end
end


text:
using SONames
using Test

@test !isempty(SONames.SONAMES[])
empty!(SONames.SONAMES[])
@test isempty(SONames.SONAMES[])
SONames.read()
@test !isempty(SONames.SONAMES[])

@test SONames.lookup("libc") == "/lib/libc.so.7"
@test SONames.lookup("libyoloswag") == ""


text:
# using CSV, Dates, HDF5, Statistics, Parameters

hydrooptions() = Dict(
    :gisregion => "Europe8",                    # "Europe8", "Eurasia38", "Scand3"

    :costclasses_min => [ 0,  50,  100],         # US $/MWh
    :costclasses_max => [50, 100, 9999],

    :storageclasses_min => [   0, 1e-6,  12],   # weeks (discharge time)
    :storageclasses_max => [1e-6,   12, 9e9]
)

mutable struct HydroOptions
    gisregion               ::String
    costclasses_min         ::Vector{Float64}
    costclasses_max         ::Vector{Float64}
    storageclasses_min      ::Vector{Float64}
    storageclasses_max      ::Vector{Float64}
end

HydroOptions() = HydroOptions("",[],[],[],[])

function HydroOptions(d::Dict{Symbol,Any})
    options = HydroOptions()
    for (key,val) in d
        setproperty!(options, key, val)
    end
    return options
end

function getregion(lon, lat, regions, lonrange=1:36000, latrange=1:18000)
    res = 0.01
    res2 = res/2
    lons = (-180+res2:res:180-res2)[lonrange]         # longitude values (pixel center)
    lats = (90-res2:-res:-90+res2)[latrange]          # latitude values (pixel center)
    flon = findfirst(round(lon-res2, digits=5) .< lons .<= round(lon+res2, digits=5))
    flat = findfirst(round(lat-res2, digits=5) .< lats .<= round(lat+res2, digits=5))
    (flon == nothing || flat == nothing) && return 0
    # println(lons[flon],", ",lats[flat])
    # println(flon, " ", flat)
    if regions[flon, flat] > 0
        return regions[flon, flat]
    else
        return maximum(regions[max(1,flon-2):min(length(lons),flon+2), max(1,flat-2):min(length(lats),flat+2)])
    end
end

function GIShydro(; optionlist...)
    options = HydroOptions(merge(hydrooptions(), optionlist))

    potential, existing, hydroplants, WECcapacity, WECpotential = readhydrodatabases()

    dayspermonth, nordicprofile, eleccost, capacity, reservoirenergy, dischargetime, monthlyinflow =
        potentialhydrovars(potential)

    countries, countrynames, regions, regionlist, lonrange, latrange = read_countries_and_regions(options)

    existingcapac, existinginflowcf =
        existinginflow(existing, countries, countrynames, regions, regionlist, lonrange, latrange, hydroplants, WECcapacity,
                        dayspermonth)

    potentialcapac, potentialinflowcf, potentialmeancost, potentialmeandischargetime =
        potentialinflow(options, potential, capacity, regions, regionlist, lonrange, latrange, monthlyinflow, eleccost,
                        dischargetime, dayspermonth)

    println("\nSaving...")

    @unpack gisregion = options
    mkpath(in_datafolder("output"))
    matopen(in_datafolder("output", "GISdata_hydro_$gisregion.mat"), "w") do file
        write(file, "existingcapac", existingcapac)
        write(file, "existinginflowcf", existinginflowcf)
        write(file, "potentialcapac", potentialcapac)
        write(file, "potentialinflowcf", potentialinflowcf)
        write(file, "potentialmeancost", potentialmeancost)
        write(file, "potentialmeandischargetime", potentialmeandischargetime)
    end
    nothing
end

function readhydrodatabases()
    println("\nReading hydro databases...")

    # lat,lon,COE,Production_GWh,Lake_surface_m2,Lake_volume_m3,
    #   Qm1,Qm2,Qm3,Qm4,Qm5,Qm6,Qm7,Qm8,Qm9,Qm10,Qm11,Qm12,Qm13,ContID,BasinID,SysID,CapCost
    # original headers:
    # lat,lon,COE ($/kWh),Production (GWh),Lake surface (m2),Lake volume (m3),
    #   Qm1,Qm2,Qm3,Qm4,Qm5,Qm6,Qm7,Qm8,Qm9,Qm10,Qm11,Qm12,Qm13,ContID,BasinID,
    #   SysID (1=DiversionalCanalPower/2=RiverPower), CapCost ($/kW)
    potential = DataFrame(CSV.File(in_datafolder("Hydro database (Gernaat) - potential.csv")))

    # GrandID,lat,lon,Production_kWh,m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11,m12,m13
    # original headers:
    # GrandID,lat,lon,Production_GWh,m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11,m12,m13
    existing = DataFrame(CSV.File(in_datafolder("Hydro database (Gernaat) - existing (GRanD).csv")))

    # country,country_long,name,gppd_idnr,capacity_mw,latitude,longitude,fuel1,fuel2,fuel3,fuel4,commissioning_year,owner,source,url,geolocation_source,year_of_capacity_data,generation_gwh_2013,generation_gwh_2014,generation_gwh_2015,generation_gwh_2016,estimated_generation_gwh
    elecplants = DataFrame(CSV.File(in_datafolder("WRI - Global Power Plant Database v1.10", "global_power_plant_database.csv")), copycols=true)

    # clean up wrong coordinates
    wrong = findall((elecplants.fuel1 .== "Hydro") .& (
            (elecplants.latitude.>90) .| (elecplants.latitude.<-90) .| (elecplants.longitude.>180) .| (elecplants.longitude.<-180)
        ))
    # elecplants.url[wrong]     # 8 results, check urls online
    # 1: DDMMSS in lat and lon
    # 2-8: lat & lon flipped
    w = wrong[1]
    elecplants.latitude[w], elecplants.longitude[w] = dms2deg(elecplants.latitude[w]), dms2deg(elecplants.longitude[w])
    w = wrong[2:8]
    elecplants.latitude[w], elecplants.longitude[w] = elecplants.longitude[w], elecplants.latitude[w]
    wrong = findall((elecplants.fuel1 .== "Hydro") .& (
            (elecplants.latitude.>90) .| (elecplants.latitude.<-90) .| (elecplants.longitude.>180) .| (elecplants.longitude.<-180)
        ))
    if !isempty(wrong)
        error("Cleanup of power plant database coordinates failed.")
    end

    # capacity_mw,latitude,longitude
    hydroplants = filter(r -> !ismissing(r.fuel1) && r.fuel1 == "Hydro", elecplants)[:,5:7]

    # Country, Capacity_MW, Pumped_MW, Other_MW, Generation_GWh, Generation_BP_GWh   % [0 = no data]
    # original headers:
    # Country, Total Hydropower Capacity (MW) in 2015, Pumped Storage Capacity (MW) in 2015, Excluding Pumped Storage (MW) in 2015, Estimated Net Hydropower Generation (GWh) in 2015, Consumption of Hydroelectricity (GWh) in 2015 (BP 2016)
    WECcapacity = CSV.File(in_datafolder("WEC hydro capacity 2015.csv")) |> DataFrame

    # Country, Undeveloped_GWh, Potential_GWh, Utilisation
    # original headers:
    # Country, Undeveloped (GWh/year), Total Potential (GWh/year), Current Utilisation (%)
    WECpotential = CSV.File(in_datafolder("WEC hydro potentials.csv")) |> DataFrame

    return potential, existing, hydroplants, WECcapacity, WECpotential
end

function potentialhydrovars(potential)
    dayspermonth = [Dates.daysinmonth(Date("2018-$m")) for m = 1:12]

    # No data on Nordic sites (hydro data limited to 60 degrees north), so take Nordic monthly inflow from a European model instead.
    nordicprofile = [2224, 1746, 2032, 4089, 15163, 23469, 16126, 9740, 9026, 9061, 5438, 3544]

    production = potential[:, :Production_GWh]              # GWh
    profile = Matrix(potential[:, 7:18])                    # m3/s  (columns Qm1-Qm12)
    eleccost = potential[:, :COE]                           # $/kWh
    capcost = potential[:, :CapCost]                        # $/kW
    type = potential[:, :SysID]
    reservoirsize = potential[:, :Lake_volume_m3] / 1e6     # Mm3

    cfriver = mean(profile, dims=2) ./ maximum(profile, dims=2)
    crf = CRF(0.1, 40)                                      # Gernaat assumes 10% discount rate and 40 year lifetime
    cf = capcost * crf/8760 ./ eleccost
    effic = 0.9*(type .== 1) .+ 0.7*(type .== 2)
    capacity = production ./ cf * 1000/8760                 # MW
    sortedflow = sort(profile, dims=2)                      # m3/s
    designflow = sortedflow[:, 9]                           # m3/s
    waterdensity = 997                                      # kg/m3
    grav = 9.81                                             # m/s2

    # My estimate of fallheight using eq 2-3 in the Gernaat paper.
    # Maybe ask for a new version of the database with this field included.
    fallheight = production*1e9/8760 ./ cf ./ effic ./ (designflow*waterdensity*grav)     # m

    # remove outliers
    f = findall(fallheight .> 2000)
    fallheight[f] .= 2000
    production[f] = fallheight[f] .* cf[f] .* effic[f] .* designflow[f] * waterdensity*grav*8760/1e9
    capacity[f] = production[f] ./ cf[f] * 1000/8760
    eleccost[f] = eleccost[f] .* potential[f, :Production_GWh] ./ production[f]
    capcost[f] = capcost[f] .* potential[f, :Production_GWh] ./ production[f]

    # J = kg m2/s2
    # m * m3/s * kg/m3 * m/s2 = J/s = W
    energyprofile = fallheight * dayspermonth' * 24/1e9 .* profile * waterdensity*grav    # GWh/month
    waterenergy = sum(energyprofile, dims=2)    # GWh/year

    # m * Mm3 * kg/m3 * m/s2 = MJ
    reservoirenergy = fallheight .* reservoirsize * waterdensity*grav / 3600 / 1000     # GWh
    dischargetime = reservoirenergy ./ capacity * 1000     # h
    # Matlab plots:  figure; hist(log10(reservoirenergy(reservoirenergy>0)),1000)
    # hist(log10(dischargetime(reservoirenergy>0)),1000)

    # Try calculating a monthly cf based on water inflow and turbine size.
    # Maybe not always accurate for large reservoirs since they can store between months,
    # but should work for run-of-river. Annual production and CF should still be fine.
    # (But why is my production estimate sometimes much higher than Gernaat's?)
    monthlyinflow = energyprofile .* effic
    monthlyprod = min.(monthlyinflow, capacity * dayspermonth' * 24/1000)   # GWh/month

    return dayspermonth, nordicprofile, eleccost, capacity, reservoirenergy, dischargetime, monthlyinflow
end

function read_countries_and_regions(options)
    println("Reading Global GADM country-level database...")
    countries, _, countrynames, _, _ = loadregions("Global_GADM0")
    numcountries = length(countrynames)

    @unpack gisregion = options
    println("Reading regions for $gisregion...")
    regions, offshoreregions, regionlist, lonrange, latrange = loadregions(gisregion)
    # bbox = [-90 -180; 90 180]

    return countries, countrynames, regions, regionlist, lonrange, latrange
end

function existinginflow(existing, countries, countrynames, regions, regionlist, lonrange, latrange, hydroplants, WECcapacity,
                        dayspermonth)
    println("\nCalculate monthly inflow for existing hydro...");

    existingprofile = Matrix(existing[:, 5:16])             # m3/s  (columns m1-m12)

    numcountries = length(countrynames)
    capacwri = zeros(numcountries,1)        # GW
    annualcf = zeros(numcountries,1)
    scalefactor_capacwri = zeros(numcountries,1)

    for i = 1:size(hydroplants,1)
        country = getregion(hydroplants[i,:longitude], hydroplants[i,:latitude], countries)
        i == 4867 && continue   # weird data point in the Indian Ocean
        if country > 0 && country != NOREGION
            capacwri[country] += hydroplants[i,:capacity_mw] / 1e3    # GW
        else
            println("Can't identify country: $(hydroplants[i,:])")
        end
    end

    # calculate scale factor for national hydro capacity (WRI capac/WEC capac)
    # calculate annual cf for each country from WEC data
    for i = 1:size(WECcapacity,1)-1     # last row is "World"
        cty = findfirst(countrynames .== Symbol(WECcapacity[i,:Country]))
        cty == nothing && error("Can't find country name: $(Symbol(WECcapacity[i,:Country]))")
        scalefactor_capacwri[cty] = capacwri[cty] / WECcapacity[i,:Capacity_MW] * 1000
        annualcf[cty] = WECcapacity[i,:Generation_GWh] / WECcapacity[i,:Capacity_MW] * 1000/8760
    end

    numreg = length(regionlist)
    existingcapac = zeros(numreg,1)         # GW
    existinginflow = zeros(numreg,12)       # GWh/month

    # calculate annual generation at each site using capacity and scale factor & CF from above
    # then get monthly inflow profile from nearest site in Gernaat/Grand
    # finally calculate monthly inflow for existing sites
    for i = 1:size(hydroplants,1)
        cty = getregion(hydroplants[i,:longitude], hydroplants[i,:latitude], countries)
        reg = getregion(hydroplants[i,:longitude], hydroplants[i,:latitude], regions, lonrange, latrange)
        (reg == 0 || reg == NOREGION || cty == 0 || cty == NOREGION) && continue
        scalefactor = scalefactor_capacwri[cty] > 0 ? scalefactor_capacwri[cty] : 1.0
        existingcapacity = hydroplants[i,:capacity_mw] / 1e3 / scalefactor      # GW
        annualgeneration = existingcapacity * annualcf[cty] * 8760              # GWh/year
        existingcapac[reg] += existingcapacity

        # approx distance in km to all existing plants in Gernaat/Grand
        dist = 111*sqrt.((existing[:,:lon] .- hydroplants[i,:longitude]).^2 + (existing[:,:lat] .- hydroplants[i,:latitude]).^2)
        d, ndx = findmin(dist)
        if regionlist[reg] == "NOR"
            inflowprofile = nordicprofile
        else
            inflowprofile = existingprofile[ndx,:]
        end
        monthlyprofile = inflowprofile .* dayspermonth
        monthlyprofile = monthlyprofile / sum(monthlyprofile)

        monthlygeneration = annualgeneration * monthlyprofile       # GWh/month
        existinginflow[reg,:] += monthlygeneration
    end

    existinginflowcf = existinginflow ./ existingcapac ./ dayspermonth' / 24

    return existingcapac, existinginflowcf
end

function potentialinflow(options, potential, capacity, regions, regionlist, lonrange, latrange, monthlyinflow, eleccost,
            dischargetime, dayspermonth)
    println("Calculate monthly inflow for potential hydro...")

    @unpack costclasses_min, costclasses_max, storageclasses_min, storageclasses_max, gisregion = options

    ncostclasses = length(costclasses_min)
    nstorageclasses = length(storageclasses_min)

    numreg = length(regionlist)
    potentialcapac = zeros(numreg,ncostclasses,nstorageclasses)        # GW
    potentialinflow = zeros(numreg,ncostclasses,nstorageclasses,12)    # GWh/month

    potentialmeancost = zeros(numreg,ncostclasses,nstorageclasses)             # $/kWh
    potentialmeandischargetime = zeros(numreg,ncostclasses,nstorageclasses)    # hours
    nobservations = zeros(numreg,ncostclasses,nstorageclasses)

    for i = 1:size(potential,1)
        reg = getregion(potential[i,:lon], potential[i,:lat], regions, lonrange, latrange)
        (reg == 0 || reg == NOREGION) && continue

        weeks = dischargetime[i] / 168
        storageclass = findfirst((weeks .>= storageclasses_min) .& (weeks .< storageclasses_max))
        cost = eleccost[i]*1000     # US $/MWh
        costclass = findfirst((cost .>= costclasses_min) .& (cost .< costclasses_max))

        if regionlist[reg] == "NOR"
            inflowprofile = nordicprofile/sum(nordicprofile) * sum(monthlyinflow[i,:])
        else
            inflowprofile = monthlyinflow[i,:]
        end

        potentialcapac[reg,costclass,storageclass] += capacity[i]/1000
        potentialinflow[reg,costclass,storageclass,:] += inflowprofile

        potentialmeancost[reg,costclass,storageclass] += eleccost[i]
        potentialmeandischargetime[reg,costclass,storageclass] += min(8760, dischargetime[i])
        nobservations[reg,costclass,storageclass] += 1
    end

    potentialinflowcf = potentialinflow ./ potentialcapac ./ reshape(dayspermonth, (1,1,1,12)) / 24
    potentialmeancost = potentialmeancost ./ nobservations
    potentialmeandischargetime = potentialmeandischargetime ./ nobservations

    return potentialcapac, potentialinflowcf, potentialmeancost, potentialmeandischargetime
end

# function mapurl(latlon1,latlon2)
#     url = "https://www.google.com/maps/dir/?api=1&origin=$(latlon1[1]),$(latlon1[2])" *
#             "&destination=$(latlon2[1]),$(latlon2[2])&travelmode=walking"
# end

max_stars_repo_name: henrikhodel/GlobalEnergyGIS

max_issues_repo_name: henrikhodel/GlobalEnergyGIS

max_forks_repo_name: henrikhodel/GlobalEnergyGIS


text:
# partial constraint development
using Test
using IncrementalInference
# going to introduce two new constraint mutable structs
import IncrementalInference: getSample
using Statistics

##

mutable struct DevelopPartial <: AbstractPrior
  x::Distribution
  partial::Tuple # should rather be static types or templates for performance Tuple{Int, Int} etc.
end
getSample(cf::CalcFactor{<:DevelopPartial}, N::Int=1) = (reshape(rand(cf.factor.x, N),1,N), )

mutable struct DevelopDim2 <: AbstractPrior
  x::Distribution
end
getSample(cf::CalcFactor{<:DevelopDim2}, N::Int=1) = (rand(cf.factor.x, N), )

##

N=100 # 50
fg = initfg()

v1 = addVariable!(fg,:x1,ContinuousEuclid{2}(),N=N)

pr = DevelopDim2(MvNormal([0.0;0.0], diagm([0.01;0.01]))) # *Matrix{Float64}(LinearAlgebra.I,2,2)))
f1  = addFactor!(fg,[:x1],pr)

dp = DevelopPartial(Normal(2.0, 1.0),(1,))
f2  = addFactor!(fg,[:x1],dp)

@testset "test evaluation of full constraint prior" begin
  pts = evalFactor(fg, f1, v1.label, N=N)
  @test size(pts,1) == 2
  @test size(pts,2) == N
  @test norm(Statistics.mean(pts,dims=2)[1] .- [0.0]) < 0.3
end

memcheck = getVal(v1)

@testset "test evaluation of partial constraint prior" begin

X1pts = getVal(v1)
pts = approxConv(fg, f2.label, :x1, N=N)

@test size(pts, 1) == 2
@test size(pts, 2) == N
@test norm(Statistics.mean(pts,dims=2)[1] .- [2.0]) < 0.75
# ensure the correct response from
@test norm(X1pts[1,:] - pts[1,:]) > 2.0
@test norm(X1pts[2,:] - pts[2,:]) < 1e-10
@test norm(X1pts - memcheck) < 1e-10

end

@testset "test solving of factor graph" begin
  getSolverParams(fg).N = N
  tree, smt, hist = solveTree!(fg)
  pts = getVal(fg, :x1)
  @test norm(Statistics.mean(pts,dims=2)[1] .- [0.0]) < 0.25
  @test norm(Statistics.mean(pts,dims=2)[2] .- [0.0]) < 0.25
end
# plotKDE(getBelief(fg, :x1),levels=3)

mutable struct DevelopPartialPairwise <: AbstractRelativeMinimize
  x::Distribution
  partial::Tuple
  DevelopPartialPairwise(x::Distribution) = new(x, (2,))
end
getSample(cf::CalcFactor{<:DevelopPartialPairwise}, N::Int=1) = (reshape(rand(cf.factor.x, N),1,N), )

function (dp::CalcFactor{<:DevelopPartialPairwise})(meas,
                                                    x1,
                                                    x2  )
  #
  # v0.21+
  return meas[1] - (x2[2]-x1[2])
end

v2 = addVariable!(fg,:x2,ContinuousEuclid{2},N=N)

dpp = DevelopPartialPairwise(Normal(10.0, 1.0))
f3  = addFactor!(fg,[:x1;:x2],dpp)

dp2 = DevelopPartial( Normal(-20.0, 1.0), (1,) )
f4  = addFactor!(fg,[:x2;], dp2, graphinit=false)
doautoinit!(fg, :x2)

# drawGraph(fg, show=true)

##

@testset "test evaluation of multiple simultaneous partial constraints" begin
global fg

ensureAllInitialized!(fg)
valx2 = getVal(fg, :x2)
pts = approxConv(fg, :x1x2f1, :x2, N=N) # evalFactor(fg, f3, v2.index, N=N)
@test size(pts,1) == 2
@test norm(Statistics.mean(pts,dims=2)[2] .- [10.0]) < 3.0
@test norm(valx2[1,:] - pts[1,:]) < 1e-5

pts = approxConv(fg, :x2f1, :x2, N=N) # evalFactor(fg, f4, v2.index, N=N)
@test size(pts,1) == 2
@test norm(Statistics.mean(pts,dims=2)[1] .- [-20.0]) < 0.75
@test (Statistics.std(pts,dims=2)[1] .- 1.0) < 0.4

end

# keep previous values to ensure funciton evaluation is modifying correct data fields

@warn "restore findRelatedFromPotential as testset!"
# @testset "test findRelatedFromPotential..." begin
# global v2, fg, f3, f4, N

thefac = getFactor(fg, :x1x2f1)

X2lpts = getVal(getVariable(fg, :x2))
keepaside, = findRelatedFromPotential(fg, thefac, :x2, N=N)
@test Ndim(keepaside) == 2
lpts = KernelDensityEstimate.getPoints(keepaside)
@test size(lpts,2) == N

@show X2lpts[2,95:100]
@show lpts[2,95:100]
@show getPoints(keepaside)

# DevelopPartialPairwise must only modify the second dimension of proposal distribution on X2
@test norm(X2lpts[1,:] - lpts[1,:]) < 1e-10
# @test norm(X2lpts[2,:] - lpts[2,:]) > 1e-10 # 10*N*0.5 # why so big?
memcheck = getVal(v2)
@test norm(X2lpts - memcheck) < 1e-10

X2lpts = getVal(v2)
p4, = findRelatedFromPotential(fg, f4, v2.label, N=N)
@test Ndim(p4) == 2
lpts = KernelDensityEstimate.getPoints(keepaside)
@test size(lpts,2) == N

# DevelopPartialPairwise must only modify the second dimension of proposal distribution on X2
@test norm(X2lpts[1,:] - lpts[1,:]) < 1e-10
@test norm(X2lpts[2,:] - lpts[2,:]) > 1e-10 # 10*N*0.5 # why so big?
memcheck = getVal(v2)
@test norm(X2lpts - memcheck) < 1e-10

# end

@testset "test belief prediction with partials..." begin

global v2, fg

# partial prior
X2pts = getVal(v2)
val, = predictbelief(fg, v2, [f4], N=N)
@test norm(X2pts[2,:] - val[2,:]) < 1e-10
@test 0.0 < norm(X2pts[1,:] - val[1,:])
@test norm(Statistics.mean(val[1,:]) .+ 20.0) < 0.75

# partial pairwise
X2pts = getVal(v2)
val, = predictbelief(fg, v2, [f3], N=N)
@test norm(X2pts[1,:] - val[1,:]) < 1e-10
@test 0.0 < norm(X2pts[2,:] - val[2,:])
@test abs(Statistics.mean(val[2,:] - getVal(v1)[2,:]) .- 10.0) < 0.75

# combination of partials
val, = predictbelief(fg, v2, [f3;f4], N=N)
# plotKDE(kde!(val),levels=3)
@test norm(Statistics.mean(val,dims=2)[1] .- [-20.0]) < 2.0
@test norm(Statistics.mean(val,dims=2)[2] .- [10.0]) < 2.0
@test (Statistics.std(val,dims=2)[1] .- 1.0) < 3.0
@test (Statistics.std(val,dims=2)[2] .- 1.0) < 3.0

getSolverParams(fg).N = N
tree, smt, hist = solveTree!(fg)

pts = getVal(fg, :x1)
@test norm(Statistics.mean(pts,dims=2)[1] .- [0.0]) < 0.5
@test norm(Statistics.mean(pts,dims=2)[2] .- [0.0]) < 0.5

pts = getVal(fg, :x2)
@test norm(Statistics.mean(pts,dims=2)[1] .- [-20.0]) < 2.0
@test norm(Statistics.mean(pts,dims=2)[2] .- [10.0]) < 2.0
@test (Statistics.std(pts,dims=2)[1]-1.0) < 3.0
@test (Statistics.std(pts,dims=2)[2]-1.0) < 3.0

end

# plotKDE(getBelief(fg, :x2),levels=3)

# spyCliqMat(tree, :x2)

#


text:
struct POMCPOWTree{B,A,O,RB}
    # action nodes
    n::Vector{Int}
    v::Vector{Float64}
    generated::Vector{Vector{Pair{O,Int}}}
    a_child_lookup::Dict{Tuple{Int,O}, Int} # may not be maintained based on solver params
    a_labels::Vector{A}
    n_a_children::Vector{Int}

    # observation nodes
    sr_beliefs::Vector{B} # first element is #undef
    total_n::Vector{Int}
    tried::Vector{Vector{Int}}
    o_child_lookup::Dict{Tuple{Int,A}, Int} # may not be maintained based on solver params
    o_labels::Vector{O}

    # root
    root_belief::RB

    function POMCPOWTree{B,A,O,RB}(root_belief, sz::Int=1000) where{B,A,O,RB}
        sz = min(sz, 100_000)
        return new(
            sizehint!(Int[], sz),
            sizehint!(Int[], sz),
            sizehint!(Vector{Pair{O,Int}}[], sz),
            Dict{Tuple{Int,O}, Int}(),
            sizehint!(A[], sz),
            sizehint!(Int[], sz),

            sizehint!(Array{B}(undef, 1), sz),
            sizehint!(Int[0], sz),
            sizehint!(Vector{Int}[Int[]], sz),
            Dict{Tuple{Int,A}, Int}(),
            sizehint!(Array{O}(undef, 1), sz),

            root_belief
        )
    end
end

@inline function push_anode!(tree::POMCPOWTree{B,A,O}, h::Int, a::A, n::Int=0, v::Float64=0.0, update_lookup=true) where {B,A,O}
    anode = length(tree.n) + 1
    push!(tree.n, n)
    push!(tree.v, v)
    push!(tree.generated, Pair{O,Int}[])
    push!(tree.a_labels, a)
    push!(tree.n_a_children, 0)
    if update_lookup
        tree.o_child_lookup[(h, a)] = anode
    end
    push!(tree.tried[h], anode)
    tree.total_n[h] += n
    return anode
end

struct POWTreeObsNode{B,A,O,RB} <: BeliefNode
    tree::POMCPOWTree{B,A,O,RB}
    node::Int
end

isroot(h::POWTreeObsNode) = h.node==1
@inline function belief(h::POWTreeObsNode)
    if isroot(h)
        return h.tree.root_belief
    else
        return StateBelief(h.tree.sr_beliefs[h.node])
    end
end
function sr_belief(h::POWTreeObsNode)
    if isroot(h)
        error("Tried to access the sr_belief for the root node in a POMCPOW tree")
    else
        return h.tree.sr_beliefs[h.node]
    end
end
n_children(h::POWTreeObsNode) = length(h.tree.tried[h.node])


text:
# CuArrays development often happens in lockstep with other packages, so try to match branches
if haskey(ENV, "GITLAB_CI")
  using Pkg
  function match_package(package, branch)
    try
      Pkg.add(PackageSpec(name=package, rev=String(branch)))
      @info "Installed $package from $branch branch"
    catch ex
      @warn "Could not install $package from $branch branch, trying master" exception=ex
      Pkg.add(PackageSpec(name=package, rev="master"))
      @info "Installed $package from master branch"
    end
  end

  branch = ENV["CI_COMMIT_REF_NAME"]
  for package in ("GPUArrays", "CUDAnative")
    match_package(package, branch)
  end
end

using Test

include("util.jl")

using Random
Random.seed!(1)

using CuArrays

using GPUArrays
import GPUArrays: allowscalar, @allowscalar

testf(f, xs...; kwargs...) = GPUArrays.TestSuite.compare(f, CuArray, xs...; kwargs...)

allowscalar(false)

@testset "CuArrays" begin

include("base.jl")
isdefined(CuArrays, :CUDNN)     && include("dnn.jl")
isdefined(CuArrays, :CUBLAS)    && include("blas.jl")
isdefined(CuArrays, :CUSOLVER)  && include("solver.jl")
isdefined(CuArrays, :CUFFT)     && include("fft.jl")
isdefined(CuArrays, :CURAND)    && include("rand.jl")

end


text:
#ДАНО: Робот - в произвольной клетке поля (без внутренних перегородок и маркеров)
#РЕЗУЛЬТАТ: Робот - в исходном положении, и все клетки по периметру внешней рамки промакированы
function granica(r::Robot)
    q = 0
    while isborder(r,HorizonSide(0)) == false
        move!(r,HorizonSide(0))
        q = q+1
    end
    for side in (HorizonSide(mod(i,4)) for i=1:4)
        putmarkers(r,side)
    end
    move!(r,HorizonSide(1))
    while ismarker(r) == false
        putmarker!(r)
        move!(r,HorizonSide(1))
    end
    move!(r,HorizonSide(3))
    for i in 1:q
        move!(r,HorizonSide(2))
    end
end
function putmarkers(r::Robot,side::HorizonSide)
    while isborder(r,side)==false
        move!(r,side)
        putmarker!(r)
    end
end

max_stars_repo_name: listova-nastya/mirea-progs

max_issues_repo_name: listova-nastya/mirea-progs

max_forks_repo_name: listova-nastya/mirea-progs


text:
using Spectra
using Test

# testing long correction function
x_for_long = [20.,21.,22.,23.,24.,25.]
y_for_long = [1.0,1.0,1.0,1.0,1.0,1.0]
h_plank = 6.626070040e-34 #Plank constant
k_bolt = 1.38064852e-23  #Boltzman constant
c_light = 299792458. #speed of light in vaccum
nu0 = 1.0./(514.532).*1e9 #laser wavenumber at 514.532
nu = 100.0.*x_for_long # cm-1 to m-1
T = 23.0+273.15 # the temperature in K

x_long,long_res,eselong = tlcorrection([x_for_long y_for_long],23.0,514.532) # using the function
t0 = nu0.^3.0 .*nu./((nu0.-nu).^4)
t1= 1 .- exp.(-h_plank.*c_light.*nu./(k_bolt.*T)) # c in m/s  : t1 dimensionless
long_calc= y_for_long.*t0.*t1 # pour les y
long_calc = long_calc./trapz(x_for_long,long_calc) # area normalisation

@test isapprox(long_res,long_calc)
@test isapprox(x_for_long,x_long)


text:
@testset "Scattering Transform Definition" begin
Transform = ScatteringTransform1d(2, [12, 1], [8, 8])
@test Scat.depth(Transform)==2
end


text:
module TestFF

using FF
using BitIntegers, Test

@testset "defenitions" begin
  @test @isdefined ff
  @test ff <: Integer
end

@testset "arithmetics" begin
  @test ff(1)/2 == ff(uint256"0x183227397098d014dc2822db40c0ac2e9419f4243cdcb848a1f0fac9f8000001")
  @test ff(-1) == ff(uint256"0x30644e72e131a029b85045b68181585d2833e84879b9709143e1f593f0000000")
end

end #module


text:
module GenomicTiles
using BioFeatures
using DataFrames
using IntervalTrees

#Tiles
export make_tiles
export calculate_tile_value
export set_tile_value
export set_tile_value1
export add_tiles_to_tiles_matrix
export add_tiles_to_predefined_tiles_matrix
export Tile
#Cluster
export cluster_tiles
export order_vector
#Plots
export plot_matrix
export plot_matrix2
export plot_cluster_bar_chart

include(Pkg.dir("GenomicTiles","src","Tiles.jl"))
include(Pkg.dir("GenomicTiles","src","Cluster.jl"))
include(Pkg.dir("GenomicTiles","src","Plots.jl"))
end


text:
"""
    sma(x::Array{Float64}; n::Int64=10)::Array{Float64}

Simple moving average (SMA)
"""
function sma(x::Array{Float64}; n::Int64=10)::Array{Float64}
    return runmean(x, n=n, cumulative=false)
end

"""
    trima(x::Array{Float64}; n::Int64=10, ma::Function=sma, args...)::Array{Float64}

Triangular moving average (TRIMA)
"""
function trima(x::Array{Float64}; n::Int64=10, ma::Function=sma)::Array{Float64}
    return ma(ma(x, n=n), n=n)
end

"""
    wma(x::Array{Float64}; n::Int64=10, wts::Array{Float64}=collect(1:n)/sum(1:n))::Array{Float64}

Weighted moving average (WMA)
"""
function wma(x::Array{Float64}; n::Int64=10, wts::Array{Float64}=collect(1:n)/sum(1:n))
    @assert n<size(x,1) && n>0 "Argument n out of bounds"
    out = fill(NaN, size(x,1))
    @inbounds for i = n:size(x,1)
        out[i] = (wts' * x[i-n+1:i])[1]
    end
    return out
end

function first_valid(x::Array{Float64})::Int
    if !isnan(x[1])
        return 1
    else
        @inbounds for i in 2:length(x)
            if !isnan(x[i])
                return i
            end
        end
    end
    return 0
end

"""
    ema(x::Array{Float64}; n::Int64=10, alpha::Float64=2.0/(n+1.0), wilder::Bool=false)::Array{Float64}

Exponential moving average (EMA)
"""
function ema(x::Array{Float64}; n::Int64=10, alpha::Float64=2.0/(n+1), wilder::Bool=false)
    @assert n<size(x,1) && n>0 "Argument n out of bounds."
    if wilder
        alpha = 1.0/n
    end
    out = zeros(size(x))
    i = first_valid(x)
    out[1:n+i-2] .= NaN
    out[n+i-1] = mean(x[i:n+i-1])
    @inbounds for i = n+i:size(x,1)
        out[i] = alpha * (x[i] - out[i-1]) + out[i-1]
    end
    return out
end

"""
    mma(x::Array{Float64}; n::Int64=10)::Array{Float64}

Modified moving average (MMA)
"""
function mma(x::Array{Float64}; n::Int64=10)
    return ema(x, n=n, alpha=1.0/n)
end

"""
    dema(x::Array{Float64}; n::Int64=10, alpha=2.0/(n+1), wilder::Bool=false)::Array{Float64}

Double exponential moving average (DEMA)
"""
function dema(x::Array{Float64}; n::Int64=10, alpha=2.0/(n+1), wilder::Bool=false)
    return 2.0 * ema(x, n=n, alpha=alpha, wilder=wilder) -
        ema(ema(x, n=n, alpha=alpha, wilder=wilder),
            n=n, alpha=alpha, wilder=wilder)
end

"""
    tema(x::Array{Float64}; n::Int64=10, alpha=2.0/(n+1), wilder::Bool=false)::Array{Float64}

Triple exponential moving average (TEMA)
"""
function tema(x::Array{Float64}; n::Int64=10, alpha=2.0/(n+1), wilder::Bool=false)
    return 3.0 * ema(x, n=n, alpha=alpha, wilder=wilder) -
        3.0 * ema(ema(x, n=n, alpha=alpha, wilder=wilder),
                  n=n, alpha=alpha, wilder=wilder) +
        ema(ema(ema(x, n=n, alpha=alpha, wilder=wilder),
                n=n, alpha=alpha, wilder=wilder),
            n=n, alpha=alpha, wilder=wilder)
end

"""
    mama(x::Array{Float64}; fastlimit::Float64=0.5, slowlimit::Float64=0.05)::Matrix{Float64}

MESA adaptive moving average (MAMA)
"""
function mama(x::Array{Float64}; fastlimit::Float64=0.5, slowlimit::Float64=0.05)::Matrix{Float64}
    n = size(x,1)
    out = zeros(Float64, n, 2)
    #smooth = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    smooth_1 = 0.0
    smooth_2 = 0.0
    smooth_3 = 0.0
    smooth_4 = 0.0
    smooth_5 = 0.0
    smooth_6 = 0.0
    smooth_7 = 0.0
    #detrend = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    detrend_1 = 0.0
    detrend_2 = 0.0
    detrend_3 = 0.0
    detrend_4 = 0.0
    detrend_5 = 0.0
    detrend_6 = 0.0
    detrend_7 = 0.0
    #Q1 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    Q1_1 = 0.0
    Q1_2 = 0.0
    Q1_3 = 0.0
    Q1_4 = 0.0
    Q1_5 = 0.0
    Q1_6 = 0.0
    Q1_7 = 0.0
    #I1 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    I1_1 = 0.0
    I1_2 = 0.0
    I1_3 = 0.0
    I1_4 = 0.0
    I1_5 = 0.0
    I1_6 = 0.0
    I1_7 = 0.0
    #I2 = [0.0, 0.0]
    I2_1 = 0.0
    I2_2 = 0.0
    #Q2 = [0.0, 0.0]
    Q2_1 = 0.0
    Q2_2 = 0.0
    #Re = [0.0, 0.0]
    Re_1 = 0.0
    Re_2 = 0.0
    #Im = [0.0, 0.0]
    Im_1 = 0.0
    Im_2 = 0.0
    #per = [0.0, 0.0]
    per_1 = 0.0
    per_2 = 0.0
    #sper = [0.0, 0.0]
    sper_1 = 0.0
    sper_2 = 0.0
    #phase = [0.0, 0.0]
    phase_1 = 0.0
    phase_2 = 0.0
    jI = 0.0
    jQ = 0.0
    dphase = 0.0
    alpha = 0.0
    a = 0.0962
    b = 0.5769
    c = 0.075
    d = 0.54
    @inbounds for i = 13:n
        # Smooth and detrend price movement ====================================
        smooth_7 = (4*x[i] + 3*x[i-1] + 2*x[i-2] + x[i-3]) * 0.1
        detrend_7 = (a*smooth_7+b*smooth_5-b*smooth_3-a*smooth_1) * (c*per_1+d)
        # Compute InPhase and Quandrature components ===========================
        Q1_7 = (a*detrend_7+b*detrend_5-b*detrend_3-a*detrend_1) * (c*per_1+d)
        I1_7 = detrend_4
        # Advance phase of I1 and Q1 by 90 degrees =============================
        jQ = (a*Q1_7+b*Q1_5-b*Q1_3-a*Q1_1) * (c*per_1+d)
        jI = (a*I1_7+b*I1_5-b*I1_3-a*I1_1) * (c*per_1+d)
        # Phasor addition for 3 bar averaging ==================================
        Q2_2 = Q1_7 + jI
        I2_2 = I1_7 - jQ
        # Smooth I & Q components before applying the discriminator ============
        Q2_2 = 0.2 * Q2_2 + 0.8 * Q2_1
        I2_2 = 0.2 * I2_2 + 0.8 * I2_1
        # Homodyne discriminator ===============================================
        Re_2 = I2_2 * I2_1 + Q2_2*Q2_1
        Im_2 = I2_2 * Q2_1 - Q2_2*I2_1
        Re_2 = 0.2 * Re_2 + 0.8*Re_1
        Im_2 = 0.2 * Im_2 + 0.8*Im_1
        if (Im_2 != 0.0) & (Re_2 != 0.0)
            per_2 = 360.0/atan(Im_2/Re_2)
        end
        if per_2 > 1.5 * per_1
            per_2 = 1.5*per_1
        elseif per_2 < 0.67 * per_1
            per_2 = 0.67 * per_1
        end
        if per_2 < 6.0
            per_2 = 6.0
        elseif per_2 > 50.0
            per_2 = 50.0
        end
        per_2 = 0.2*per_2 + 0.8*per_1
        sper_2 = 0.33*per_2 + 0.67*sper_1
        if I1_7 != 0.0
            phase_2 = atan(Q1_7/I1_7)
        end
        dphase = phase_1 - phase_2
        if dphase < 1.0
            dphase = 1.0
        end
        alpha = fastlimit / dphase
        if alpha < slowlimit
            alpha = slowlimit
        end
        out[i,1] = alpha*x[i] + (1.0-alpha)*out[i-1,1]
        out[i,2] = 0.5*alpha*out[i,1] + (1.0-0.5*alpha)*out[i-1,2]
        # Reset/increment array variables
        # smooth
        smooth_1 = smooth_2
        smooth_2 = smooth_3
        smooth_3 = smooth_4
        smooth_4 = smooth_5
        smooth_5 = smooth_6
        smooth_6 = smooth_7
        # detrend
        detrend_1 = detrend_2
        detrend_2 = detrend_3
        detrend_3 = detrend_4
        detrend_4 = detrend_5
        detrend_5 = detrend_6
        detrend_6 = detrend_7
        # Q1
        Q1_1 = Q1_2
        Q1_2 = Q1_3
        Q1_3 = Q1_4
        Q1_4 = Q1_5
        Q1_5 = Q1_6
        Q1_6 = Q1_7
        # I1
        I1_1 = I1_2
        I1_2 = I1_3
        I1_3 = I1_4
        I1_4 = I1_5
        I1_5 = I1_6
        I1_6 = I1_7
        # I2
        I2_1 = I2_2
        # Q2
        Q2_1 = Q2_2
        # Re
        Re_1 = Re_2
        # Im
        Im_1 = Im_2
        # per
        per_1 = per_2
        # sper
        sper_1 = sper_2
        # phase
        phase_1 = phase_2
    end
    out[1:32,:] .= NaN
    return out
end

"""
    hma(x::Array{Float64}; n::Int64=20)::Array{Float64}

Hull moving average (HMA)
"""
function hma(x::Array{Float64}; n::Int64=20)
    return wma(2 * wma(x, n=Int64(round(n/2.0))) - wma(x, n=n), n=Int64(trunc(sqrt(n))))
end

"""
Sine-weighted moving average
"""
function swma(x::Array{Float64}; n::Int64=10)
    @assert n<size(x,1) && n>0 "Argument n out of bounds."
    w = sin.(collect(1:n) * 180.0/6.0)  # numerator weights
    d = sum(w)  # denominator = sum(numerator weights)
    out = zeros(size(x))
    out[1:n-1] .= NaN
    @inbounds for i = n:size(x,1)
        out[i] = sum(w .* x[i-n+1:i]) / d
    end
    return out
end

"""
Kaufman adaptive moving average (KAMA)
"""
function kama(x::Array{Float64}; n::Int64=10, nfast::Float64=0.6667, nslow::Float64=0.0645)
    @assert n<size(x,1) && n>0 "Argument n out of bounds."
    @assert nfast>0.0 && nfast<1.0 "Argument nfast out of bounds."
    @assert nslow>0.0 && nslow<1.0 "Argument nslow out of bounds."
    dir = diffn(x, n=n)  # price direction := net change in price over past n periods
    vol = runsum(abs.(diffn(x,n=1)), n=n, cumulative=false)  # volatility/noise
    er = abs.(dir) ./ vol  # efficiency ratio
    ssc = er * (nfast-nslow) .+ nslow  # scaled smoothing constant
    sc = ssc .^ 2  # smoothing constant
    # initiliaze result variable
    out = zeros(size(x))
    i = ndims(x) > 1 ? findfirst(.!isnan.(x)).I[1] : findfirst(.!isnan.(x))
    out[1:n+i-2] .= NaN
    out[n+i-1] = mean(x[i:n+i-1])
    @inbounds for i = n+1:size(x,1)
        out[i] = out[i-1] + sc[i]*(x[i]-out[i-1])
    end
    return out
end

"""
    alma{Float64}(x::Array{Float64}; n::Int64=9, offset::Float64=0.85, sigma::Float64=6.0)::Array{Float64}

Arnaud-Legoux moving average (ALMA)
"""
function alma(x::Array{Float64}; n::Int64=9, offset::Float64=0.85, sigma::Float64=6.0)
    @assert n<size(x,1) && n>0 "Argument n out of bounds."
    @assert sigma>0.0 "Argument sigma must be greater than 0."
    @assert offset>=0.0 && offset<=1 "Argument offset must be in (0,1)."
    out = zeros(size(x))
    out[1:n-1] .= NaN
    m = floor(offset*(float(n)-1.0))
    s = float(n) / sigma
    w = exp.(-(((0.0:-1.0:-float(n)+1.0) .- m) .^ 2.0) / (2.0*s*s))
    wsum = sum(w)
    if wsum != 0.0
        w = w ./ wsum
    end
    @inbounds for i = n:length(x)
        out[i] = sum(x[i-n+1] .* w)
    end
    return out
end

function lagged(x::Array{Float64}, n::Int=1)::Array{Float64}
    if n > 0
        return [fill(NaN,n); x[1:end-n]]
    elseif n < 0
        return [x[(-n+1):end]; fill(NaN,-n)]
    else
        return x
    end
end

"""
    zlema(x::Array{Float64}; n::Int=10, ema_args...)::Array{Float64}

Zero-lag exponential moving average (ZLEMA)
"""
function zlema(x::Array{Float64}; n::Int=10, ema_args...)::Array{Float64}
    return ema(x+(x-lagged(x,round(Int, (n-1)/2.0))), n=n; ema_args...)
end


text:
using Revise, Parameters, KrylovKit
using BifurcationKit
using LinearAlgebra, SparseArrays, LinearMaps, DiffEqOperators, Setfield
const BK = BifurcationKit

using GLMakie
GLMakie.inline!(true)

contour3dMakie(x; k...) = GLMakie.contour(x;  k...)
contour3dMakie(x::AbstractVector; k...) = contour3dMakie(reshape(x,Nx,Ny,Nz); k...)
contour3dMakie(ax, x; k...) = (contour(ax, x;  k...))
contour3dMakie(ax, x::AbstractVector; k...) = contour3dMakie(ax, reshape(x,Nx,Ny,Nz); k...)
contour3dMakie!(ax, x; k...) = (AbstractPlotting.contour!(ax, x;  k...))
contour3dMakie!(ax, x::AbstractVector; k...) = contour3dMakie!(ax, reshape(x,Nx,Ny,Nz); k...)

function Laplacian3D(Nx, Ny, Nz, lx, ly, lz, bc = :Neumann)
	speye(n) = sparse(I, n, n)
	hx = 2lx/Nx; hy = 2ly/Ny; hz = 2lz/Nz
	D2x = CenteredDifference{1}(2, 2, hx, Nx)
	D2y = CenteredDifference{1}(2, 2, hy, Ny)
	D2z = CenteredDifference{1}(2, 2, hz, Nz)
	Qx = Neumann0BC(hx); Qy = Neumann0BC(hy); Qz = Neumann0BC(hz)

	_A = kron(speye(Ny), sparse(D2x * Qx)[1]) + kron(sparse(D2y * Qy)[1], speye(Nx))
	A = kron(speye(Nz), _A) + kron(kron(sparse(D2z * Qz)[1], speye(Ny)), speye(Nx))
	return sparse(A), D2x
end

function F_sh(u, p)
	@unpack l, ν, L1 = p
	return -(L1 * u) .+ (l .* u .+ ν .* u.^2 .- u.^3)
end

function dF_sh(u, p, du)
	@unpack l, ν, L1 = p
	return -(L1 * du) .+ (l .+ 2 .* ν .* u .- 3 .* u.^2) .* du
end

# function J_sh(u, p)
# 	@unpack l, ν, L1 = p
# 	return -L1 .+ spdiagm(0 => l .+ 2 .* ν .* u .- 3 .* u.^2)
# end

d2F_sh(u, p, dx1, dx2) = (2 .* p.ν .* dx2 .- 6 .* dx2 .* u) .* dx1
d3F_sh(u, p, dx1, dx2, dx3) = (-6 .* dx2 .* dx3) .* dx1
jet = (F_sh, (x, p) -> (dx -> dF_sh(x, p, dx)), d2F_sh, d3F_sh)

const TY = Float64
AF = Array{TY}
####################################################################################################
Nx = Ny = Nz = 22; N = Nx*Ny*Nz
lx = ly = lz = pi

X = -lx .+ 2lx/(Nx) * collect(0:Nx-1)
Y = -ly .+ 2ly/(Ny) * collect(0:Ny-1)
Z = -lz .+ 2lz/(Nz) * collect(0:Nz-1)

sol0 = [(cos(x) .* cos(y )) for x in X, y in Y, z in Z]
	sol0 .= sol0 .- minimum(vec(sol0))
	sol0 ./= maximum(vec(sol0))
	sol0 = sol0 .- 0.0
	sol0 .*= 1.7

Δ, D2x = @time Laplacian3D(Nx, Ny, Nz, lx, ly, lz, :Neumann)
L1 = (I + Δ)^2
par = (l = 0.1, ν = 1.2, L1 = L1)
Pr = @time cholesky(L1)
using SuiteSparse
LinearAlgebra.ldiv!(P::SuiteSparse.CHOLMOD.Factor{Float64}, v) = -(P \ v)
# LinearAlgebra.ldiv!(o, P::SuiteSparse.CHOLMOD.Factor{Float64}, v) = o .= P \ v

# rtol must be small enough to pass the folds and to get precise eigenvalues
ls = GMRESKrylovKit(verbose = 0, rtol = 1e-9, maxiter = 150, ishermitian = true, Pl = Pr)
####################################################################################################
struct SH3dEig{Ts, Tσ} <: BK.AbstractEigenSolver
	ls::Ts
	σ::Tσ
end

BifurcationKit.geteigenvector(eigsolve::SH3dEig, vecs, n::Union{Int, Array{Int64,1}}) = vecs[n]

function (sheig::SH3dEig)(J, nev::Int; verbosity = 0, kwargs...)
	σ = sheig.σ
	nv = 30
	Jshift = du -> J(du) .- σ .* du
	A = du -> sheig.ls(Jshift, du)[1]
	# we adapt the krylov dimension as function of the requested eigenvalue number
	vals, vec, info = KrylovKit.eigsolve(A, AF(rand(Nx*Ny*Nz)), nev, :LM, tol = 1e-12, maxiter = 20, verbosity = verbosity, ishermitian = true, krylovdim = max(nv, nev + nv))
	vals2 = 1 ./vals .+ σ
	Ind = sortperm(vals2, by = real, rev = true)
	return vals2[Ind], vec[Ind], true, info.numops
end

eigSH3d = SH3dEig((@set ls.rtol = 1e-9), 0.1)

optnew = NewtonPar(verbose = true, tol = 1e-8, maxIter = 20, linsolver= @set ls.verbose=0)
	@set! optnew.eigsolver = eigSH3d
	sol_hexa, hist, flag = @time BK.newton(F_sh,
		(x, p) -> (dx -> dF_sh(x, p, dx)),
		AF(vec(sol0)), par, optnew)
	println("--> norm(sol) = ", norm(sol_hexa, Inf64))

contour3dMakie(sol0)
contour3dMakie(sol_hexa)
###################################################################################################
optcont = ContinuationPar(dsmin = 0.0001, dsmax = 0.005, ds= -0.001, pMax = 0.15, pMin = -.1, newtonOptions = setproperties(optnew; tol = 1e-9, maxIter = 15), maxSteps = 146, detectBifurcation = 3, nev = 15, nInversion = 4, plotEveryStep  = 1)

	br, u1 = @time continuation(
		F_sh, #J_sh,
		(x, p) -> (dx -> dF_sh(x, p, dx)),
		AF(zeros(N)), par, (@lens _.l), optcont;
		plot = true, verbosity = 3,
		plotSolution = (ax, x, p) -> contour3dMakie(ax, x),
		recordFromSolution = (x, p) -> (n2 = norm(x), n8 = norm(x, 8)),
		normC = x -> norm(x, Inf),
		event = BK.FoldDetectEvent,
		tangentAlgo = BorderedPred(),

BK.plotBranch(br)
contour3dMakie(u1.u)
####################################################################################################
computeNormalForm(jet..., br, 2; issymmetric = true)

br1, = @time continuation(jet..., br, 3, setproperties(optcont; saveSolEveryStep = 10, detectBifurcation = 0, pMax = 0.1, plotEveryStep = 5, dsmax = 0.01);
	plot = true, verbosity = 3,
	δp = 0.005,
	verbosedeflation = false,
	tangentAlgo = BorderedPred(),
	issymmetric = true,
	plotSolution = (ax, x, p) -> contour3dMakie(ax, x),
	recordFromSolution = (x, p) -> (n2 = norm(x), n8 = norm(x, 8)),
	finaliseSolution = (z, tau, step, contResult; k...) -> begin
		if isnothing(br.eig) == true
			Base.display(contResult.eig[end].eigenvals)
		end
		true
	end,
	# callbackN = cb,
	normC = x -> norm(x, Inf))

BK.plotBranch(br,br1...)
BK.plotBranch(br1[15])

BK.plotBranch(br1...)

fig = Figure(resolution = (1200, 900))
	for i=1:min(25,length(br1))
		ix = div(i,5)+1; iy = i%5+1
		@show i, ix, iy
		ax = Axis3(fig[ix, iy], title = "$i", aspect = (1, 1, 1))
		hidedecorations!(ax, grid=false)
		contour3dMakie!(ax, br1[i].sol[2].x)
		ax.protrusions = (0, 0, 0, 10)
		# out = AbstractPlotting.contour!(ax, reshape(br1[i].sol[2].x, Nx, Ny, Nz))
		# @show out
		# Colorbar(fig[ix, iy], )
		# Colorbar(ax, hm)
	end
	display(fig)


text:
#= test/sketch.jl
=#

@testset "sketch" begin
  m = 20
  n = 10
  M = rand(ComplexF64, m, n)
  r = 5
  opts = LRAOptions(sketch_randn_niter=1)

  @time for (t, s) in ((:RandomGaussian,       :randn),
                       (:RandomSubset,         :sub  ),
                       (:SRFT,                 :srft ),
                       (:SparseRandomGaussian, :sprn ))
    opts.sketch = s
    for T in (Float32, Float64, ComplexF32, ComplexF64)
      let A
          println("  $t/$T")

          A = convert(Array{T}, T <: Real ? real(M) : M)

          S = sketch(:left,  :n, A, r, opts)
          @test size(S) == (r, n)
          S = sketch(:left,  :c, A, r, opts)
          @test size(S) == (r, m)
          S = sketch(:right, :n, A, r, opts)
          @test size(S) == (m, r)
          S = sketch(:right, :c, A, r, opts)
          @test size(S) == (n, r)
      end
    end
  end
end


text:
push!(LOAD_PATH, "../")

using Jive
runtests(@__DIR__, skip=["revise.jl"])

max_stars_repo_name: HenricoWitvliet/FWFTables