text:
using SIMD
using Test, InteractiveUtils

"""
    llvm_ir(f, args) :: String

Get LLVM IR of `f(args...)` as a string.
"""
llvm_ir(f, args) = sprint(code_llvm, f, Base.typesof(args...))

@testset "SIMD" begin

        # The vector we are testing. Ideally, we should be able to use any vector size
        # anywhere, but LLVM codegen bugs prevent us from doing so -- thus we make this
        # a parameter.
        global const nbytes = 32

        global const L8 = nbytes÷4
        global const L4 = nbytes÷8

        global const V8I32 = Vec{L8,Int32}
        global const V4F64 = Vec{L4,Float64}

        is_checking_bounds = Core.Compiler.inbounds_option() == :on

    @testset "Type properties" begin
        @test eltype(V8I32) === Int32
        @test eltype(V4F64) === Float64
        @test length(V8I32) == L8
        @test length(V4F64) == L4
        @test ndims(V8I32) == 1
        @test ndims(V4F64) == 1
        @test size(V8I32,1) == L8
        @test size(V4F64,1) == L4
        @test size(V8I32) == (L8,)
        @test size(V4F64) == (L4,)
    end

    @testset "Type conversion" begin

        global const v8i32 = ntuple(i->Int32(ifelse(isodd(i), i, -i)), L8)
        global const v4f64 = ntuple(i->Float64(ifelse(isodd(i), i, -i)), L4)

        @test string(V8I32(v8i32)) == "<8 x Int32>[" * string(v8i32)[2:end-1] * "]"
        @test string(V4F64(v4f64)) == "<4 x Float64>[" * string(v4f64)[2:end-1] * "]"

        @test convert(V8I32, V8I32(v8i32)) === V8I32(v8i32)
        @test convert(Vec{L8,Int64}, V8I32(v8i32)) ===
            Vec{L8, Int64}(convert(NTuple{L8,Int64}, v8i32))

        @test NTuple{L8,Int32}(V8I32(v8i32)) === Tuple(v8i32)
        @test NTuple{L4,Float64}(V4F64(v4f64)) === Tuple(v4f64)
        @test Tuple(V8I32(v8i32)) === Tuple(v8i32)
        @test Tuple(V4F64(v4f64)) === Tuple(v4f64)
    end

    @testset "Element-wise access" begin

        for i in 1:L8
            @test Tuple(setindex(V8I32(v8i32), 9.0, Val(i))) ===
                ntuple(j->Int32(ifelse(j==i, 9, v8i32[j])), L8)
            @test Tuple(setindex(V8I32(v8i32), 9.0, Val{i})) ===
                ntuple(j->Int32(ifelse(j==i, 9, v8i32[j])), L8)
            @test Tuple(setindex(V8I32(v8i32), 9.0, i)) ===
                ntuple(j->Int32(ifelse(j==i, 9, v8i32[j])), L8)

            @test V8I32(v8i32)[Val{i}] === v8i32[i]
            @test V8I32(v8i32)[i] === v8i32[i]
        end

        @test_throws BoundsError setindex(V8I32(v8i32), 0, Val(0))
        @test_throws BoundsError setindex(V8I32(v8i32), 0, Val{0})
        @test_throws BoundsError setindex(V8I32(v8i32), 0, Val(L8+1))
        @test_throws BoundsError setindex(V8I32(v8i32), 0, Val{L8+1})
        @test_throws BoundsError setindex(V8I32(v8i32), 0, 0)
        @test_throws BoundsError setindex(V8I32(v8i32), 0, L8+1)
        @test_throws BoundsError V8I32(v8i32)[Val(0)]
        @test_throws BoundsError V8I32(v8i32)[Val{0}]
        @test_throws BoundsError V8I32(v8i32)[Val(L8+1)]
        @test_throws BoundsError V8I32(v8i32)[Val{L8+1}]
        @test_throws BoundsError V8I32(v8i32)[0]
        @test_throws BoundsError V8I32(v8i32)[L8+1]

        for i in 1:L4
            @test Tuple(setindex(V4F64(v4f64), 9, Val(i))) ===
                ntuple(j->Float64(ifelse(j==i, 9.0, v4f64[j])), L4)
            @test Tuple(setindex(V4F64(v4f64), 9, Val{i})) ===
                ntuple(j->Float64(ifelse(j==i, 9.0, v4f64[j])), L4)
            @test Tuple(setindex(V4F64(v4f64), 9, i)) ===
                ntuple(j->Float64(ifelse(j==i, 9.0, v4f64[j])), L4)

            @test V4F64(v4f64)[Val(i)] === v4f64[i]
            @test V4F64(v4f64)[Val{i}] === v4f64[i]
            @test V4F64(v4f64)[i] === v4f64[i]
        end

        let
            v0 = zero(Vec{4, Float64})
            v1 = one(Vec{4, Float64})
            @test sum(v0*v0) == 0.0
            @test sum(v1*v1) == 4.0
        end
    end

    @testset "Integer arithmetic functions" begin

        global const v8i32b = map(x->Int32(x+1), v8i32)
        global const v8i32c = map(x->Int32(x*2), v8i32)

        notbool(x) = !(x>=typeof(x)(0))
        for op in (~, +, -, abs, notbool, sign, signbit)
            @test Tuple(op(V8I32(v8i32))) === map(op, v8i32)
        end

        for op in (
                +, -, *, ÷, %, ==, !=, <, <=, >, >=,
                copysign, div, flipsign, max, min, rem)
            @test Tuple(op(V8I32(v8i32), V8I32(v8i32b))) === map(op, v8i32, v8i32b)
        end

        global vifelsebool(x,y,z) = vifelse(x>=typeof(x)(0),y,z)
        for op in (vifelsebool, muladd)
            @test Tuple(op(V8I32(v8i32), V8I32(v8i32b), V8I32(v8i32c))) ===
                map(op, v8i32, v8i32b, v8i32c)
        end

        for op in (<<, >>, >>>)
            @test Tuple(op(V8I32(v8i32), Val(3))) === map(x->op(x,3), v8i32)
            @test Tuple(op(V8I32(v8i32), Val{3})) === map(x->op(x,3), v8i32)
            @test Tuple(op(V8I32(v8i32), Val(-3))) === map(x->op(x,-3), v8i32)
            @test Tuple(op(V8I32(v8i32), Val{-3})) === map(x->op(x,-3), v8i32)
            @test Tuple(op(V8I32(v8i32), 3)) === map(x->op(x,3), v8i32)
            @test Tuple(op(V8I32(v8i32), -3)) === map(x->op(x,-3), v8i32)
            @test Tuple(op(V8I32(v8i32), V8I32(v8i32))) === map(op, v8i32, v8i32)
        end

        @test Tuple(V8I32(v8i32)^0) === v8i32.^0
        @test Tuple(V8I32(v8i32)^1) === v8i32.^1
        @test Tuple(V8I32(v8i32)^2) === v8i32.^2
        @test Tuple(V8I32(v8i32)^3) === v8i32.^3
    end

    @testset "Floating point arithmetic functions" begin

        global const v4f64b = map(x->Float64(x+1), v4f64)
        global const v4f64c = map(x->Float64(x*2), v4f64)

        logabs(x) = log(abs(x))
        log10abs(x) = log10(abs(x))
        log2abs(x) = log2(abs(x))
        powi4(x) = x^4
        sqrtabs(x) = sqrt(abs(x))
        for op in (
                +, -,
                abs, ceil, inv, isfinite, isinf, isnan, issubnormal, floor, powi4,
                round, sign, signbit, sqrtabs, trunc)
            @test Tuple(op(V4F64(v4f64))) === map(op, v4f64)
        end
        function Base.isapprox(t1::Tuple,t2::Tuple)
            length(t1)==length(t2) &&
                all(Bool[isapprox(t1[i], t2[i]) for i in 1:length(t1)])
        end
        for op in (cos, exp, exp10, exp2, logabs, log10abs, log2abs, sin)
            rvec = Tuple(op(V4F64(v4f64)))
            rsca = map(op, v4f64)
            @test typeof(rvec) === typeof(rsca)
            @test isapprox(rvec, rsca)
        end

        @test isfinite(V4F64(0.0))[1]
        @test isfinite(V4F64(-0.0))[1]
        @test isfinite(V4F64(nextfloat(0.0)))[1]
        @test isfinite(V4F64(-nextfloat(0.0)))[1]
        @test isfinite(V4F64(1.0))[1]
        @test isfinite(V4F64(-1.0))[1]
        @test !isfinite(V4F64(Inf))[1]
        @test !isfinite(V4F64(-Inf))[1]
        @test !isfinite(V4F64(NaN))[1]

        @test !isinf(V4F64(0.0))[1]
        @test !isinf(V4F64(-0.0))[1]
        @test !isinf(V4F64(nextfloat(0.0)))[1]
        @test !isinf(V4F64(-nextfloat(0.0)))[1]
        @test !isinf(V4F64(1.0))[1]
        @test !isinf(V4F64(-1.0))[1]
        @test isinf(V4F64(Inf))[1]
        @test isinf(V4F64(-Inf))[1]
        @test !isinf(V4F64(NaN))[1]

        @test !isnan(V4F64(0.0))[1]
        @test !isnan(V4F64(-0.0))[1]
        @test !isnan(V4F64(nextfloat(0.0)))[1]
        @test !isnan(V4F64(-nextfloat(0.0)))[1]
        @test !isnan(V4F64(1.0))[1]
        @test !isnan(V4F64(-1.0))[1]
        @test !isnan(V4F64(Inf))[1]
        @test !isnan(V4F64(-Inf))[1]
        @test isnan(V4F64(NaN))[1]

        @test !issubnormal(V4F64(0.0))[1]
        @test !issubnormal(V4F64(-0.0))[1]
        @test issubnormal(V4F64(nextfloat(0.0)))[1]
        @test issubnormal(V4F64(-nextfloat(0.0)))[1]
        @test !issubnormal(V4F64(1.0))[1]
        @test !issubnormal(V4F64(-1.0))[1]
        @test !issubnormal(V4F64(Inf))[1]
        @test !issubnormal(V4F64(-Inf))[1]
        @test !issubnormal(V4F64(NaN))[1]

        @test !signbit(V4F64(0.0))[1]
        @test signbit(V4F64(-0.0))[1]
        @test !signbit(V4F64(nextfloat(0.0)))[1]
        @test signbit(V4F64(-nextfloat(0.0)))[1]
        @test !signbit(V4F64(1.0))[1]
        @test signbit(V4F64(-1.0))[1]
        @test !signbit(V4F64(Inf))[1]
        @test signbit(V4F64(-Inf))[1]
        @test !signbit(V4F64(NaN))[1]

        for op in (
                +, -, *, /, %, ^, ==, !=, <, <=, >, >=,
                copysign, flipsign, max, min, rem)
            @test Tuple(op(V4F64(v4f64), V4F64(v4f64b))) === map(op, v4f64, v4f64b)
        end

        @test Tuple(V4F64(v4f64)^0) === v4f64.^0
        @test Tuple(V4F64(v4f64)^1) === v4f64.^1
        @test Tuple(V4F64(v4f64)^2) === v4f64.^2
        @test Tuple(V4F64(v4f64)^3) === v4f64.^3

        for op in (fma, vifelsebool, muladd)
            @test Tuple(op(V4F64(v4f64), V4F64(v4f64b), V4F64(v4f64c))) ===
                map(op, v4f64, v4f64b, v4f64c)
        end

        v = V4F64(v4f64)
        @test v^5 === v * v * v * v * v

        # Make sure our dispatching rule does not select floating point `pow`.
        # See: https://github.com/eschnett/SIMD.jl/pull/43
        ir = llvm_ir(^, (V4F64(v4f64), 2))
        @test occursin("@llvm.powi.v4f64", ir)
    end

    @testset "Type promotion" begin

        for op in (
                ==, !=, <, <=, >, >=,
                &, |, ⊻, +, -, *, copysign, div, flipsign, max, min, rem)
            @test op(42, V8I32(v8i32)) === op(V8I32(42), V8I32(v8i32))
            @test op(V8I32(v8i32), 42) === op(V8I32(v8i32), V8I32(42))
        end
        @test vifelse(signbit(V8I32(v8i32)), 42, V8I32(v8i32)) ===
            vifelse(signbit(V8I32(v8i32)), V8I32(42), V8I32(v8i32))
        @test vifelse(signbit(V8I32(v8i32)), V8I32(v8i32), 42) ===
            vifelse(signbit(V8I32(v8i32)), V8I32(v8i32), V8I32(42))
        for op in (muladd,)
            @test op(42, 42, V8I32(v8i32)) ===
                op(V8I32(42), V8I32(42), V8I32(v8i32))
            @test op(42, V8I32(v8i32), V8I32(v8i32)) ===
                op(V8I32(42), V8I32(v8i32), V8I32(v8i32))
            @test op(V8I32(v8i32), 42, V8I32(v8i32)) ===
                op(V8I32(v8i32), V8I32(42), V8I32(v8i32))
            @test op(V8I32(v8i32), V8I32(v8i32), 42) ===
                op(V8I32(v8i32), V8I32(v8i32), V8I32(42))
            @test op(42, V8I32(v8i32), 42) ===
                op(V8I32(42), V8I32(v8i32), V8I32(42))
            @test op(V8I32(v8i32), 42, 42) ===
                op(V8I32(v8i32), V8I32(42), V8I32(42))
        end

        for op in (
                ==, !=, <, <=, >, >=,
                +, -, *, /, copysign, flipsign, max, min, rem)
            @test op(42, V4F64(v4f64)) === op(V4F64(42), V4F64(v4f64))
            @test op(V4F64(v4f64), 42) === op(V4F64(v4f64), V4F64(42))
        end
        @test vifelse(signbit(V4F64(v4f64)), 42, V4F64(v4f64)) ===
            vifelse(signbit(V4F64(v4f64)), V4F64(42), V4F64(v4f64))
        @test vifelse(signbit(V4F64(v4f64)), V4F64(v4f64), 42) ===
            vifelse(signbit(V4F64(v4f64)), V4F64(v4f64), V4F64(42))
        for op in (fma, muladd)
            @test op(42, 42, V4F64(v4f64)) ===
                op(V4F64(42), V4F64(42), V4F64(v4f64))
            @test op(42, V4F64(v4f64), V4F64(v4f64)) ===
                op(V4F64(42), V4F64(v4f64), V4F64(v4f64))
            @test op(V4F64(v4f64), 42, V4F64(v4f64)) ===
                op(V4F64(v4f64), V4F64(42), V4F64(v4f64))
            @test op(V4F64(v4f64), V4F64(v4f64), 42) ===
                op(V4F64(v4f64), V4F64(v4f64), V4F64(42))
            @test op(42, V4F64(v4f64), 42) ===
                op(V4F64(42), V4F64(v4f64), V4F64(42))
            @test op(V4F64(v4f64), 42, 42) ===
                op(V4F64(v4f64), V4F64(42), V4F64(42))
        end
    end

    @testset "Reduction operations" begin

        for op in (maximum, minimum, sum, prod)
            @test op(V8I32(v8i32)) === op(v8i32)
        end
        @test all(V8I32(v8i32)) == reduce(&, v8i32)
        @test any(V8I32(v8i32)) == reduce(|, v8i32)

        for op in (maximum, minimum, sum, prod)
            @test op(V4F64(v4f64)) === op(v4f64)
        end

        @test sum(Vec{3,Float64}(1)) === 3.0
        @test prod(Vec{5,Float64}(2)) === 32.0
    end

    @testset "Load and store functions" begin

        global const arri32 = valloc(Int32, L8, 2*L8) do i i end
        for i in 1:length(arri32)-(L8-1)
            @test vload(V8I32, arri32, i) === V8I32(ntuple(j->i+j-1, L8))
        end
        for i in 1:L8:length(arri32)-(L8-1)
            @test vloada(V8I32, arri32, i) === V8I32(ntuple(j->i+j-1, L8))
        end
        vstorea(V8I32(0), arri32, 1)
        vstore(V8I32(1), arri32, 2)
        for i in 1:length(arri32)
            @test arri32[i] == if i==1 0 elseif i<=(L8+1) 1 else i end
        end

        global const arrf64 = valloc(Float64, L4, 4*L4) do i i end
        for i in 1:length(arrf64)-(L4-1)
            @test vload(V4F64, arrf64, i) === V4F64(ntuple(j->i+j-1, L4))
        end
        for i in 1:4:length(arrf64)-(L4-1)
            @test vloada(V4F64, arrf64, i) === V4F64(ntuple(j->i+j-1, L4))
        end
        vstorea(V4F64(0), arrf64, 1)
        vstore(V4F64(1), arrf64, 2)
        for i in 1:length(arrf64)
            @test arrf64[i] == if i==1 0 elseif i<=(L4+1) 1 else i end
        end
    end

    @testset "Gather and scatter function" begin
        for (arr, VT) in [(arri32, V8I32), (arrf64, V4F64)]
            arr .= 1:length(arr)
            idxarr = floor.(Int, range(1, stop=length(arr), length=length(VT)))

            idx = Vec(Tuple(idxarr))
            @test vgather(arr, idx) === convert(VT, idx)
            @test vgathera(arr, idx) === convert(VT, idx)
            @test arr[idx] === convert(VT, idx)

            # Masked gather
            maskarr = zeros(Bool, length(VT))
            for i in 1:length(VT)
                maskarr[i] = true
                mask = Vec(Tuple(maskarr))
                @test vgather(arr, idx, mask) === VT(Tuple(idxarr .* maskarr))
                @test vgathera(arr, idx, mask) === VT(Tuple(idxarr .* maskarr))
                @test arr[idx, mask] === VT(Tuple(idxarr .* maskarr))
            end

            # Scatter
            varr = convert.(eltype(VT), idxarr .* 123)
            v = Vec(Tuple(varr))

            vscatter(v, fill!(arr, 0), idx)
            @test arr[idxarr] == varr
            vscattera(v, fill!(arr, 0), idx)
            @test arr[idxarr] == varr
            fill!(arr, 0)
            arr[idx] = v
            @test arr[idxarr] == varr

            # Masked scatter
            maskarr = zeros(Bool, length(VT))
            for i in 1:length(VT)
                maskarr[i] = true
                mask = Vec(Tuple(maskarr))
                vscatter(v, fill!(arr, 0), idx, mask)
                @test arr[idxarr] == varr .* maskarr
                vscattera(v, fill!(arr, 0), idx, mask)
                @test arr[idxarr] == varr .* maskarr
                fill!(arr, 0)
                arr[idx, mask] = v
                @test arr[idxarr] == varr .* maskarr
            end
        end
    end

    @testset "Index-based load/store" begin
        for (arr, VT) in [(arri32, V8I32), (arrf64, V4F64)]
            @testset "Vector ($VT)" begin
                fill!(arr.parent, 0)
                arr .= 1:length(arr)
                idx = VecRange{length(VT)}(1)
                @test arr[idx] === VT(Tuple(1:length(VT)))

                maskarr = zeros(Bool, length(VT))
                maskarr[1] = true
                mask = Vec(Tuple(maskarr))
                varr = zeros(length(VT))
                varr[1] = 1
                @test arr[idx, mask] === VT(Tuple(varr))

                @test_throws ArgumentError arr[idx, 1]
                @test_throws ArgumentError arr[idx, 1, mask]
                @test_throws ArgumentError arr[idx, mask, 1]
                @test_throws ArgumentError arr[idx, 1, mask, 1]

                lane = VecRange{length(VT)}(0)
                @test_throws BoundsError arr[lane]
                @test_throws BoundsError arr[lane + end]

                # Out-of-bound access (OK with valloc)
                slice = LinearIndices(arr)[end-length(VT)+2:end]
                varr = zeros(length(VT))
                varr[1:length(slice)] .= arr[slice]
                i = lane + first(slice)
                @test_throws BoundsError arr[i]
                f(x, i) = @inbounds x[i]
                if !is_checking_bounds
                    @test f(arr, i) === VT(Tuple(varr))
                end
            end

            arr .= 1:length(arr)
            @testset "$name" for (name, mat) in [
                        ("Matrix ($VT)", repeat(arr, outer=(1, 3))),
                        ("Matrix ($VT) with non-strided row",
                         view(repeat(arr, outer=(1, 5)), :, [2, 1, 5])),
                    ]
                idx = VecRange{length(VT)}(1)
                @test mat[idx, 1] === VT(Tuple(1:length(VT)))
                @test mat[idx, 2] === VT(Tuple(1:length(VT)))
                if mat isa SIMD.FastContiguousArray
                    @test mat[idx] === VT(Tuple(1:length(VT)))
                else
                    err = try
                        mat[idx]
                        nothing
                    catch err
                        err
                    end
                    @test err isa ArgumentError
                    @test occursin(
                        "Exactly 2 (non-mask) indices have to be specified.",
                        sprint(showerror, err))
                end

                maskarr = zeros(Bool, length(VT))
                maskarr[1] = true
                mask = Vec(Tuple(maskarr))
                varr = zeros(length(VT))
                varr[1] = 1
                @test mat[idx, 1, mask] === VT(Tuple(varr))
                if mat isa SIMD.FastContiguousArray
                    @test mat[idx, mask] === VT(Tuple(varr))
                else
                    err = try
                        mat[idx]
                        nothing
                    catch err
                        err
                    end
                    @test err isa ArgumentError
                    @test occursin(
                        "Exactly 2 (non-mask) indices have to be specified.",
                        sprint(showerror, err))
                end

                @test_throws ArgumentError mat[idx, 1, 1]
                @test_throws ArgumentError mat[idx, 1, 1, mask]
                @test_throws ArgumentError mat[idx, mask, 1]
                @test_throws ArgumentError mat[idx, 1, mask, 1]

                lane = VecRange{length(VT)}(0)
                @test_throws BoundsError mat[lane, 1]
                @test_throws BoundsError mat[lane + end, 1]
                if mat isa SIMD.FastContiguousArray
                    @test_throws BoundsError mat[lane + end]
                else
                    err = try
                        mat[lane + end]
                        nothing
                    catch err
                        err
                    end
                    @test err isa ArgumentError
                    @test occursin(
                        "Exactly 2 (non-mask) indices have to be specified.",
                        sprint(showerror, err))
                end

                # Out-of-bound access
                varr = collect(1:length(VT))
                i = lane + size(mat, 1) + 1
                @test_throws BoundsError mat[i, 1]
                f(x, i) = @inbounds x[i, 1]
                if !is_checking_bounds
                    @test f(mat, i) === VT(Tuple(varr))
                end
            end

            @testset "3D array ($VT)" begin
                arr .= 1:length(arr)
                arr3d = repeat(arr, outer=(1, 3, 5))
                idx = VecRange{length(VT)}(1)
                @test arr3d[idx, 1, 1] === VT(Tuple(1:length(VT)))
                @test arr3d[idx, 2, 1] === VT(Tuple(1:length(VT)))
                @test arr3d[idx, 1, 2] === VT(Tuple(1:length(VT)))
                @test arr3d[idx] === VT(Tuple(1:length(VT)))

                maskarr = zeros(Bool, length(VT))
                maskarr[1] = true
                mask = Vec(Tuple(maskarr))
                varr = zeros(length(VT))
                varr[1] = 1
                @test arr3d[idx, 1, 1, mask] === VT(Tuple(varr))
                @test arr3d[idx, mask] === VT(Tuple(varr))

                @test_throws ArgumentError arr3d[idx, 1, 1, 1]
                @test_throws ArgumentError arr3d[idx, 1, 1, 1, mask]
                @test_throws ArgumentError arr3d[idx, mask, 1]
                @test_throws ArgumentError arr3d[idx, 1, mask, 1]

                lane = VecRange{length(VT)}(0)
                @test_throws BoundsError arr3d[lane, 1, 1]
                @test_throws BoundsError arr3d[lane + end, 1, 1]
                @test_throws BoundsError arr3d[lane + end]

                # Out-of-bound access
                varr = collect(1:length(VT))
                i = lane + size(arr3d, 1) + 1
                @test_throws BoundsError arr3d[i, 1, 1]
                f(x, i) = @inbounds x[i, 1, 1]
                if !is_checking_bounds
                    @test f(arr3d, i) === VT(Tuple(varr))
                end
            end
        end
    end

    @testset "Real-world examples" begin

        function vadd!(xs::AbstractArray{T,1}, ys::AbstractArray{T,1},
                       ::Type{Vec{N,T}}) where {N,T}
            @assert length(ys) == length(xs)
            @assert length(xs) % N == 0
            lane = VecRange{N}(0)
            @inbounds for i in 1:N:length(xs)
                xs[lane + i] += ys[lane + i]
            end
        end

        let xs = valloc(Float64, L4, 4*L4) do i i end,
            ys = valloc(Float64, L4, 4*L4) do i 1 end
            vadd!(xs, ys, V4F64)
            @test xs == Float64[i+1 for i in 1:(4*L4)]
            # @code_native vadd!(xs, ys, V4F64)

            ir = llvm_ir(vadd!, (xs, ys, V4F64))
            @test occursin(r"( load <4 x double>.*){2}"s, ir)
            @test occursin(" store <4 x double>", ir)
            @test occursin(" fadd <4 x double>", ir)
        end

        function vsum(xs::AbstractArray{T,1}, ::Type{Vec{N,T}}) where {N,T}
            @assert length(xs) % N == 0
            sv = Vec{N,T}(0)
            lane = VecRange{N}(0)
            @inbounds for i in 1:N:length(xs)
                sv += xs[lane + i]
            end
            sum(sv)
        end

        let xs = valloc(Float64, L4, 4*L4) do i i end
            s = vsum(xs, V4F64)
            @test s === (x->(x^2+x)/2)(Float64(4*L4))
            # @code_native vsum(xs, V4F64)

            ir = llvm_ir(vsum, (xs, V4F64))
            @test occursin(" load <4 x double>", ir)
            @test occursin(" fadd <4 x double>", ir)
            @test occursin(r"( shufflevector <4 x double>.*){2}"s, ir)
        end

        function vadd_masked!(xs::AbstractArray{T,1}, ys::AbstractArray{T,1},
                              ::Type{Vec{N,T}}) where {N, T}
            @assert length(ys) == length(xs)
            limit = length(xs) - (N-1)
            vlimit = Vec(ntuple(i -> length(xs) - i + 1, Val(N)))
            lane = VecRange{N}(0)
            @inbounds for i in 1:N:length(xs)
                if i <= limit
                    xs[lane + i] += ys[lane + i]
                else
                    mask = Vec{N,Int}(i) <= vlimit
                    xs[lane + i, mask] = xs[lane + i, mask] + ys[lane + i, mask]
                end
            end
        end

        let xs = valloc(Float64, 4, 13) do i i end,
            ys = valloc(Float64, 4, 13) do i 1 end
            vadd_masked!(xs, ys, V4F64)
            @test xs == Float64[i+1 for i in 1:13]
            # @code_native vadd!(xs, ys, V4F64)

            ir = llvm_ir(vadd_masked!, (xs, ys, V4F64))
            @test occursin(r"(masked.load.v4f64.*){2}"s, ir)
            @test occursin("masked.store.v4f64", ir)
            @test occursin(" store <4 x double>", ir)
            @test occursin(" fadd <4 x double>", ir)
        end

        function vsum_masked(xs::AbstractArray{T,1}, ::Type{Vec{N,T}}) where {N,T}
            vlimit = Vec(ntuple(i -> length(xs) - i + 1, Val(N)))
            sv = Vec{N,T}(0)
            lane = VecRange{N}(0)
            @inbounds for i in 1:N:length(xs)
                mask = Vec{N,Int}(i) <= vlimit
                sv += xs[lane + i, mask]
            end
            sum(sv)
        end

        let xs = valloc(Float64, 4, 13) do i i end
            s = vsum_masked(xs, V4F64)
            @code_llvm vsum(xs, V4F64)
            @code_native vsum(xs, V4F64)
            @test s === sum(xs)

            ir = llvm_ir(vsum_masked, (xs, V4F64))
            @test occursin("masked.load.v4f64", ir)
            @test occursin(" fadd <4 x double>", ir)
            @test occursin(r"( shufflevector <4 x double>.*){2}"s, ir)
        end
    end

    @testset "Vector shuffles" begin

        for T in (Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Float32,Float64)
            a = Vec{4,T}((1,2,3,4))
            b = Vec{4,T}((5,6,7,8))
            @test shufflevector(a, b, Val((2,3,4,5))) === Vec{4,T}((3,4,5,6))
            @test shufflevector(a, b, Val((1,7,5,5))) === Vec{4,T}((2,8,6,6))
            @test shufflevector(a, b, Val(0:3)) === a
            @test shufflevector(a, b, Val(4:7)) === b
            @test shufflevector(a, Val((1,0,2,3))) === Vec{4,T}((2,1,3,4))
            @test shufflevector(a, b, Val((0,1,4,5,2,3,6,7))) === Vec{8,T}((1,2,5,6,3,4,7,8))
            @test shufflevector(shufflevector(a, b, Val((6,:undef,0,:undef))), Val((0,2))) === Vec{2,T}((7,1))
            @test isa(shufflevector(a, Val((:undef,:undef,:undef,:undef))), Vec{4,T})
            c = Vec{8,T}((1:8...,))
            d = Vec{8,T}((9:16...,))
            @test shufflevector(c, d, Val((0,1,8,15))) === Vec{4,T}((1,2,9,16))
            @test shufflevector(c, d, Val(1:2:15)) === Vec{8,T}((2:2:16...,))
        end

        let
            a = Vec{4,Bool}((true,false,true,false))
            b = Vec{4,Bool}((false,false,true,true))
            @test shufflevector(a, b, Val((2,3,4,5))) === Vec{4,Bool}((true,false,false,false))
        end

        for T in (Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Float32,Float64)
            a = Vec{4,T}((1,2,3,4))
            b = Vec{4,T}((5,6,7,8))
            @test shufflevector(a, b, Val{(2,3,4,5)}) === Vec{4,T}((3,4,5,6))
            @test shufflevector(a, b, Val{(1,7,5,5)}) === Vec{4,T}((2,8,6,6))
            @test shufflevector(a, b, Val{0:3}) === a
            @test shufflevector(a, b, Val{4:7}) === b
            @test shufflevector(a, Val{(1,0,2,3)}) === Vec{4,T}((2,1,3,4))
            @test shufflevector(a, b, Val{(0,1,4,5,2,3,6,7)}) === Vec{8,T}((1,2,5,6,3,4,7,8))
            @test shufflevector(shufflevector(a, b, Val{(6,:undef,0,:undef)}), Val{(0,2)}) === Vec{2,T}((7,1))
            @test isa(shufflevector(a, Val{(:undef,:undef,:undef,:undef)}), Vec{4,T})
            c = Vec{8,T}((1:8...,))
            d = Vec{8,T}((9:16...,))
            @test shufflevector(c, d, Val{(0,1,8,15)}) === Vec{4,T}((1,2,9,16))
            @test shufflevector(c, d, Val{1:2:15}) === Vec{8,T}((2:2:16...,))
        end

        let
            a = Vec{4,Bool}((true,false,true,false))
            b = Vec{4,Bool}((false,false,true,true))
            @test shufflevector(a, b, Val{(2,3,4,5)}) === Vec{4,Bool}((true,false,false,false))
        end
    end
end


text:
abstract Evaluator

type ProblemEvaluator <: Evaluator
  problem::OptimizationProblem
  fitness_scheme::FitnessScheme
  archive::Archive
  num_evals::Int64
  last_fitness

  ProblemEvaluator(problem::OptimizationProblem; archive = false,
    fitness_scheme = FloatFitness()) = begin
    archive = archive || TopListArchive(numdims(problem))
    new(problem, fitness_scheme, archive, 0, nothing)
  end
end

worst_fitness(e::Evaluator) = worst_fitness(e.fitness_scheme)
num_evals(e::Evaluator) = e.num_evals
numdims(e::Evaluator) = numdims(e.problem)
search_space(e::Evaluator) = search_space(e.problem)
describe(e::Evaluator) = "Problem: $(name(e.problem)) (dimensions = $(numdims(e)))"
problem_summary(e::Evaluator) = "$(name(e.problem))_$(numdims(e))d"
problem(e::Evaluator) = e.problem

function evaluate(e::Evaluator, candidate)
  e.last_fitness = eval1(candidate, e.problem)
  e.num_evals += 1
  add_candidate!(e.archive, e.last_fitness, candidate, e.num_evals)
  e.last_fitness
end

# A way to get the fitness of the last evaluated candidate. Leads to nicer
# code if we can first test if better than or worse than existing candidates
# and only want the fitness itself if the criteria fulfilled.
last_fitness(e::Evaluator) = e.last_fitness

function is_better(e::Evaluator, f1::Fitness, f2::Fitness)
  is_better(f1, f2, e.fitness_scheme)
end

function is_better(e::Evaluator, candidate, fitness::Fitness)
  is_better(evaluate(e, candidate), fitness, e.fitness_scheme)
end

function best_of(e::Evaluator, candidate1, candidate2)
  f1 = evaluate(e, candidate1)
  f2 = evaluate(e, candidate2)
  if is_better(f1, f2, e.fitness_scheme)
    return candidate1, f1
  else
    return candidate2, f2
  end
end

function rank_by_fitness(e::Evaluator, candidates)
  # Note that we re-evaluate all candidates here. This might be wasteful and
  # we should cache if evaluations are costly.
  fitness = [(c[1], c[2], evaluate(e, c[1])) for c=candidates]
  sort(fitness; by = (t) -> t[3])
end

function fitness_is_within_ftol(e::Evaluator, ftolerance; index = 1)
  try
    fitness_is_within_ftol(e.problem, ftolerance, best_fitness(e.archive); index = index)
  catch
    false # In case we have no information yet about the best fitness
  end
end

ext: jl


text:
using .Distributions
using .PyCall

"""Gaussian Bimodal"""
function gaussianBimodal(μ=0.0, σ=1.0, α=0.5)
    # μ is the mean of the satellite gaussian, σ is the width (SD) and α is the proportional probability (i.e. the mass of the satellite compared to the central gaussian)
    if rand() > α # Draw from the first distribution
        x = randn()
    else # Draw from the satellite
        x = σ*randn() + μ
    end
end

function gaussianBimodal(P::Process)
    # Sample from a distribution formed from a single unit gaussian at the origin, and a satellite gaussian with parameters (height, mean, width) you choose.
    seed(P.solver_rng)
    sol = [gaussianBimodal(parameter_function(P)(t)...) for t in P.transient_t0:P.savedt:P.tmax]
end

gaussianBimodalSim = Process(
    process = gaussianBimodal,
    parameter_profile = (constantParameter, constantParameter, constantParameter),
    parameter_profile_parameters = ((1.0,), (1.0,), (0.5)),
    transient_t0 = 0.0,
    t0 = 0.0,
    savedt = 1,
    tmax = 10000.0)
export gaussianBimodalSim

"""Bimodal Switching"""
function bimodalSwitching(t⃗, α::Function, δ::Function)
    v = zeros(length(t⃗))
    for t ∈ 2:length(t⃗)
        Y = rand(Binomial(1, α(t⃗[t])))
        v[t] = Y*v[t-1] + (1-Y)*(1-v[t-1])
    end
    x = randn(length(t⃗)) .+ v.*δ.(t⃗)
end

function bimodalSwitching(P::Process)
    # Draw from one Gaussian but at each time step have a probability 1-α of switching to another gaussian at distance δ away
    seed(P.solver_rng)
    p1, p2 = [P.parameter_profile[i](P.parameter_profile_parameters[i]...) for i ∈ 1:lastindex(P.parameter_profile)]
    sol = bimodalSwitching(P.transient_t0:P.savedt:P.tmax, p1, p2)
end

bimodalSwitchingSim = Process(
    process = bimodalSwitching,
    parameter_profile = (constantParameter, constantParameter),
    parameter_profile_parameters = ((0.5,), (3.0)),
    transient_t0 = 0.0,
    t0 = 0.0,
    savedt = 1,
    tmax = 10000.0)
export bimodalSwitchingSim

"""Shifty Noise"""
function shiftyNoise(P::Process)
    # Parameters (η, C)
    seed(P.solver_rng)
    (η, C) = parameter_functions(P)
    sol = [η(t)*randn() + C(t) for t in P.transient_t0:P.savedt:P.tmax]
end

shiftyNoiseSim = Process(
    process = shiftyNoise,
    parameter_profile = (constantParameter, stepNoise),
    parameter_profile_parameters = [(1.0,), ((0.0, 1000.0), 100.0, 2.0, 0.0)],
    transient_t0 = 0.0,
    t0 = 0.0,
    savedt = 1,
    tmax = 1000.0)
export shiftyNoiseSim

"""AR"""
function AR(P::Process)
    # Variable number of AR parameters (ϕ₁, ϕ₂, ϕ₃, ...)
    seed(P.solver_rng)
    ϕ⃗ = parameter_function(P)
    p = length(ϕ⃗(gett0(P)))
    ξ⃗ = randn(p+1)
    X = zeros(length(times(P)), p+1)
    X[1, 1:length(getX0(P))] = getX0(P)
    for t ∈ 2:length(times(P))
        X[t, :], ξ⃗ = AR(X[t-1, :], ξ⃗, forcevec(ϕ⃗(t)))
    end
    return X[Int.(times(P, transient=true) .- (gett0(P)-1)), 1:length(getX0(P))]
end

arSim = Process(
    process = AR,
    X0 = [0.0], # Number of initial conditions should really be 1 + num. of parameters, but if you do not specify they default to 0.0
    parameter_profile = (constant, constant, constant, constant),
    parameter_profile_parameters = ((0.1,), (0.1,), (0.1,), (0.1,)),
    t0 = 0,
    dt = 1,
    savedt = 1,
    tmax = 10000)
export arSim

"""
Wrap the main function of the python package colorednoise
"""
function colorednoise(β, N)
    cn.powerlaw_psd_gaussian(β, N)
end
export colorednoise

# """
# Colored noise with a power-law spectrum described by the exponent β
# """
# function coloredNoise(P::Process)
#     # Parameter β
#     seed(P.solver_rng) # Does this carry over to python? Not a huge deal if not
#     (β,) = parameter_functions(P)
#     sol = [η(t)*randn() + C(t) for t in P.transient_t0:P.savedt:P.tmax]
# end

shiftyNoiseSim = Process(
    process = shiftyNoise,
    parameter_profile = (constantParameter, stepNoise),
    parameter_profile_parameters = [(1.0,), ((0.0, 1000.0), 100.0, 2.0, 0.0)],
    transient_t0 = 0.0,
    t0 = 0.0,
    savedt = 1,
    tmax = 1000.0)
export shiftyNoiseSim


text:
test = "16,1,2,0,4,2,7,1,2,14"

function p(s::AbstractString)::Vector{Int}
    return parse.(Int, split(s, ','))
end

function minimize(f, v::AbstractVector)
    (min, max) = extrema(v)
    return minimum(c -> sum(x -> f(x, c), v), min:max)
end

function q1(a, b)
    return abs(a - b)
end

function q2(a, b)
    n = abs(a - b)
    return n * (n + 1) ÷ 2
end

let crabs = p(test)
    @assert minimize(q1, crabs) == 37
    @assert minimize(q2, crabs) == 168
end

crabs = p(readline("day07.in"))
println("Q1: ", minimize(q1, crabs))
println("Q2: ", minimize(q2, crabs))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
__precompile__(true)
module CameraUnits

import Unitful
using Unitful: @unit
using FixedPointNumbers

export γe, @cu_str

@unit γe       "γe"  PhotoElectron      1          false
@unit pco55_cu "cu"  PCO55              0.46γe     false
@unit pco42_cu "cu"  PCO42              0.46γe     false

# Base.convert(N0f16, x::PCO55) = reinterpret(N0f16, round(UInt16, ustrip(x)))
# Base.convert(N0f16, x::PCO42) = reinterpret(N0f16, round(UInt16, ustrip(x)))

include("cumacro.jl")

# Some gymnastics required here because if we precompile, we cannot add to
# Unitful.basefactors at compile time and expect the changes to persist to runtime.
const localunits = Unitful.basefactors
function __init__()
    merge!(Unitful.basefactors, localunits)
    Unitful.register(CameraUnits)
end

end # module


text:
using DelimitedFiles
using ForwardDiff

mutable struct SplinePath
    x_coefs ::Array{Float64, 2}
    y_coefs ::Array{Float64, 2}
    breaks  ::Array{Float64, 1}
end

"""
Load spline data consisting of coefficients for the polynomials
x_i(s) = a_i(s-s_i)³ + a_i(s-s_i)² + c_i(s-s_i) + d_i and
y_i(s) = e_i(s-s_i)³ + f_i(s-s_i)² + g_i(s-s_i) + h_i and breaks.
Polynomial coefficients returned have dimensions [N, 4] and breaks is a vector of length N+1.
"""
function load_splines()
    x_coefs = readdlm("../splines/nuscenes_scene_23_x_spline_coefs.csv", ',');
    y_coefs = readdlm("../splines/nuscenes_scene_23_y_spline_coefs.csv", ',');
    breaks = readdlm("../splines/nuscenes_scene_23_breaks.csv", ',')[:];

    x_coefs, y_coefs, breaks
end

"""
Performs a bisection search over interpolating spline intervals (input as the array 'breaks')
to find interval within which the input parameter 's' resides.
"""
function find_spline_interval(s::Float64, path::SplinePath)
    breaks = path.breaks

    if s < breaks[1] || s > breaks[end]
        if s < 0 && abs(s) < 1E-2
            println("The value s = $s is slightly below 0.")
            return 1
        end
        println("The value s = $s does not lie within the valid spline domain.")
        println("Path length has been overrun.")
        return
    end

    upper_idx_lim = length(breaks)
    lower_idx_lim = 1
    interval_idx = Int(div(upper_idx_lim - lower_idx_lim, 2))

    while true
        if breaks[interval_idx] <= s <= breaks[interval_idx+1]
            return interval_idx
        elseif s > breaks[interval_idx+1]
            lower_idx_lim = max(1, interval_idx + 1)
            interval_idx = Int(div(upper_idx_lim + lower_idx_lim, 2))
        elseif s < breaks[interval_idx]
            upper_idx_lim = min(length(breaks), interval_idx)
            interval_idx = Int(div(upper_idx_lim + lower_idx_lim, 2))
        else
            println("Unexpected case with s = ", s)
            throw(UndefVarError)
        end
    end
end

"""
Evaluates the x position of an interpolated spline path at arc length parameter s.
"""
function spline_x(s, path::SplinePath, spline_idx::Integer)
    x_coefs = path.x_coefs
    breaks = path.breaks

    x = x_coefs[spline_idx, 1] * (s - breaks[spline_idx]).^3 + x_coefs[spline_idx, 2] * (s - breaks[spline_idx]).^2 +
        x_coefs[spline_idx, 3] * (s - breaks[spline_idx]) + x_coefs[spline_idx, 4]
end

"""
Evaluates the x position of an interpolated spline path at arc length parameter s.
"""
function spline_y(s, path::SplinePath, spline_idx::Integer)
    y_coefs = path.y_coefs
    breaks = path.breaks

    y = y_coefs[spline_idx, 1] * (s - breaks[spline_idx]).^3 + y_coefs[spline_idx, 2] * (s - breaks[spline_idx]).^2 +
        y_coefs[spline_idx, 3] * (s - breaks[spline_idx]) + y_coefs[spline_idx, 4]
end

"""
Reprojects the vehicle position onto the path to find the appropriate arc length
parameter s. ds is resolution for arc length search.
TODO: Fix situation where multiple possible such s may be found.
TODO: Currently this only performs a local search. Add global search for when
distance between vehicle position and closes projection is greater than some
threshold.
"""
function find_best_s(q::Array{Float64,1}, path::SplinePath; ds=0.05, enable_global_search=false, sq_dist_tol=10^2)
    x = q[1]
    y = q[2]
    s = q[5]

    if s < 0
        return 0
    end

    spline_idx = find_spline_interval(s, path)
    x_proj = spline_x(s, path, spline_idx)
    y_proj = spline_y(s, path, spline_idx)
    current_dist = (x - x_proj)^2 + (y - y_proj)^2

    ss = [max(0,s-5):ds:s; s+ds:ds:s+5;] # TODO: this indexing looks a bit weird
    spline_idces = find_spline_interval.(ss, (path,))
    xs = spline_x.(ss, (path,), spline_idces)
    ys = spline_y.(ss, (path,), spline_idces)
    sq_distances = (xs .- x).^2 + (ys .- y).^2
    best_dist, best_idx = findmin(sq_distances)

    if best_dist > sq_dist_tol && enable_global_search
        println("Local search too narrow. Performing global search.")
        ss = [path_obj.breaks[1]:ds:path_obj.breaks[end];] # TODO: check this indexing
        spline_idces = find_spline_interval.(ss, (path,))
        xs = spline_x.(ss, (path,), spline_idces)
        ys = spline_y.(ss, (path,), spline_idces)
        sq_distances = (xs .- x).^2 + (ys .- y).^2
        best_dist, best_idx = findmin(sq_distances)
        if best_dist <=  sq_dist_tol
            println("Global search successfully found a better projection.")
        end
    end

    if best_dist >= current_dist
        s
    else
        ss[best_idx]
    end
end

"""
Evaluates the derivatives dx/ds and dy/ds for path at arc length parameter s0.
"""
function dpath_ds(s0, path::SplinePath, spline_idx::Integer)
    x_coefs = path.x_coefs
    y_coefs = path.y_coefs
    breaks = path.breaks

    ForwardDiff.derivative(s -> spline_x(s, path, spline_idx), s0),
        ForwardDiff.derivative(s -> spline_y(s, path, spline_idx), s0)
end

"""
Evaluates the angle of the tangent to the path at the reference point s0 with respect to the x-axis.
"""
function heading(s, path::SplinePath, spline_idx::Integer)
    x_coefs = path.x_coefs
    y_coefs = path.y_coefs
    breaks = path.breaks

    dx_ds, dy_ds = dpath_ds(s, path, spline_idx)
    atan(dy_ds, dx_ds)
end

max_stars_repo_name: StanfordASL/MATS

max_issues_repo_name: StanfordASL/MATS

max_forks_repo_name: StanfordASL/MATS


text:
__precompile__()

module GAMCSampler

using Distributions
using Klara

import Base: show

import Klara:
  RealLowerTriangular,
  format_iteration,
  format_percentage,
  generate_empty,
  initialize!,
  iterate!,
  recursive_mean!,
  reset!,
  sampler_state,
  set_diagnosticindices!,
  set_gmm,
  set_gmm!,
  tune!,
  tuner_state,
  variate_form

export
  GAMC,
  GAMCTune,
  GAMCTuner,
  GAMCState,
  MuvGAMCState,
  am_only_update!,
  cos_update!,
  exp_decay,
  iterate!,
  lin_decay,
  mod_update!,
  pow_decay,
  quad_decay,
  rand_exp_decay_update!,
  rand_lin_decay_update!,
  rand_pow_decay_update!,
  rand_quad_decay_update!,
  rand_update!,
  sampler_state,
  show,
  smmala_only_update!

include("stats/schedule.jl")
include("tuners/GAMCTuner.jl")
include("samplers/GAMC.jl")
include("samplers/iterate/GAMC.jl")

end # module


text:
"""
Baseclass for all aggregation modes.
"""
abstract type AggregateMode end

"""
    module AggMode

Types for aggregation of multiple observations.

- `AggMode.None()`
- `AggMode.Sum()`
- `AggMode.Mean()`
- `AggMode.WeightedSum(weights)`
- `AggMode.WeightedMean(weights)`
"""
module AggMode
    using ..LearnBase: AggregateMode

    """
        AggMode.None()

    Opt-out of aggregation. This is usually the default value.
    Using `None` will cause the element-wise results to be returned.
    """
    struct None <: AggregateMode end

    """
        AggMode.Sum()

    Causes the method to return the unweighted sum of the
    elements instead of the individual elements. Can be used in
    combination with `ObsDim`, in which case a vector will be
    returned containing the sum for each observation (useful
    mainly for multivariable regression).
    """
    struct Sum <: AggregateMode end

    """
        AggMode.Mean()

    Causes the method to return the unweighted mean of the
    elements instead of the individual elements. Can be used in
    combination with `ObsDim`, in which case a vector will be
    returned containing the mean for each observation (useful
    mainly for multivariable regression).
    """
    struct Mean <: AggregateMode end

    """
        AggMode.WeightedSum(weights; [normalize = false])

    Causes the method to return the weighted sum of all
    observations. The variable `weights` has to be a vector of
    the same length as the number of observations.
    If `normalize = true`, the values of the weight vector will
    be normalized in such as way that they sum to one.

    # Arguments

    - `weights::AbstractVector`: Vector of weight values that
      can be used to give certain observations a stronger
      influence on the sum.

    - `normalize::Bool`: Boolean that specifies if the weight
      vector should be transformed in such a way that it sums to
      one (i.e. normalized). This will not mutate the weight
      vector but instead happen on the fly during the
      accumulation.

      Defaults to `false`. Setting it to `true` only really
      makes sense in multivalue-regression, otherwise the result
      will be the same as for [`WeightedMean`](@ref).
    """
    struct WeightedSum{W<:AbstractVector} <: AggregateMode
        weights::W
        normalize::Bool
    end
    WeightedSum(weights::AbstractVector; normalize::Bool = false) = WeightedSum(weights, normalize)

    """
        AggMode.WeightedMean(weights; [normalize = true])

    Causes the method to return the weighted mean of all
    observations. The variable `weights` has to be a vector of
    the same length as the number of observations.
    If `normalize = true`, the values of the weight vector will
    be normalized in such as way that they sum to one.

    # Arguments

    - `weights::AbstractVector`: Vector of weight values that can
      be used to give certain observations a stronger influence
      on the mean.

    - `normalize::Bool`: Boolean that specifies if the weight
      vector should be transformed in such a way that it sums to
      one (i.e. normalized). This will not mutate the weight
      vector but instead happen on the fly during the
      accumulation.

      Defaults to `true`. Setting it to `false` only really makes
      sense in multivalue-regression, otherwise the result will
      be the same as for [`WeightedSum`](@ref).
    """
    struct WeightedMean{W<:AbstractVector} <: AggregateMode
        weights::W
        normalize::Bool
    end
    WeightedMean(weights::AbstractVector; normalize::Bool = true) = WeightedMean(weights, normalize)
end


text:
abstract type CropFrom end

struct FromOrigin <: CropFrom end
struct FromCenter <: CropFrom end
struct FromRandom <: CropFrom end

abstract type AbstractCrop <: Transform end

struct Crop{N, F<:CropFrom} <: AbstractCrop
    size::NTuple{N, Int}
    from::F
end

Crop(sz) = Crop(sz, FromOrigin())

function apply(crop::Crop, item::Item; randstate = getrandstate(crop))
    return apply(
        Project(CoordinateTransformations.IdentityTransformation()) |> crop,
        item;
        randstate = (nothing, randstate))
end

CenterCrop(sz) = Crop(sz, FromCenter())
RandomCrop(sz) = Crop(sz, FromRandom())

# The random state of a [`Crop`](#) consists of offsets from the origin.

getrandstate(crop::Crop{N, FromOrigin}) where N = Tuple(0. for _ in 1:N)
getrandstate(crop::Crop{N, FromCenter}) where N = Tuple(0.5 for _ in 1:N)
getrandstate(crop::Crop{N, FromRandom}) where N = Tuple(rand() for _ in 1:N)

struct PadDivisible <: AbstractCrop
    by::Int
end

struct CroppedProjectiveTransform{P<:ProjectiveTransform, C<:AbstractCrop} <: ProjectiveTransform
    tfm::P
    crop::C
end

function getrandstate(cropped::CroppedProjectiveTransform)
    return (getrandstate(cropped.tfm), getrandstate(cropped.crop))
end

function getprojection(
        cropped::CroppedProjectiveTransform,
        bounds;
        randstate = getrandstate(cropped))
    tfmstate, cropstate = randstate
    return getprojection(cropped.tfm, bounds; randstate = tfmstate)
end

function cropindices(
        cropped::CroppedProjectiveTransform{PT, C},
        P,
        bounds;
        randstate = getrandstate(cropped)) where {PT, C<:Crop}
    tfmstate, cropstate = randstate
    bounds_ = P.(bounds)
    return offsetcropindices(cropped.crop.size, boundsranges(bounds_), cropstate)
end

function cropindices(
        cropped::CroppedProjectiveTransform{PT, PadDivisible},
        P,
        bounds;
        randstate = getrandstate(cropped)) where {PT}
    tfmstate, cropstate = randstate
    bounds_ = P.(bounds)
    ranges = boundsranges(bounds_)

    sz = length.(ranges)
    pad = (cropped.crop.by .- (sz .% cropped.crop.by)) .% cropped.crop.by

    indices = UnitRange.(getindex.(ranges, 1), sz .+ pad)
    return indices
end

compose(tfm::ProjectiveTransform, crop::AbstractCrop) = CroppedProjectiveTransform(tfm, crop)

"""
    offsetcropindices(sz, indices, offsets)

Calculate indices for a crop of size `sz` out of `indices`.
Input and output indices are represented as tuples of ranges.

For every dimension i where `sz[i] < length(indices[i])`, offsets
the crop by `offsets[i]` times the difference between the two.
"""
function offsetcropindices(
        sz::NTuple{N, Int},
        indices::NTuple{N, <:AbstractRange},
        offsets::NTuple{N, <:Number}) where N
    mins = getindex.(indices, 1)
    diffs = length.(indices) .- sz

    startindices = floor.(Int, mins .+ (diffs .* offsets))
    endindices = startindices .+ sz .- 1

    return UnitRange.(startindices, endindices)
end


text:
"""
This type is used to store the solution to the problem presented
in the "Wald Friedman" notebook presented on the QuantEcon website.

Solution
----------
J : vector(Float64)
    Discretized value function that solves the Bellman equation
lb : scalar(Real)
    Lower cutoff for continuation decision
ub : scalar(Real)
    Upper cutoff for continuation decision
"""
type WFSolution
    J::Vector{Float64}
    lb::Real
    ub::Real
end

"""
This type is used to solve the problem presented in the "Wald Friedman"
notebook presented on the QuantEcon website.

Parameters
----------
c : scalar(Real)
    Cost of postponing decision
L0 : scalar(Real)
    Cost of choosing model 0 when the truth is model 1
L1 : scalar(Real)
    Cost of choosing model 1 when the truth is model 0
f0 : vector(Float64)
    A finite state probability distribution
f1 : vector(Float64)
    A finite state probability distribution
m : scalar(Int64)
    Number of points to use in function approximation
"""
immutable WaldFriedman
    c::Real
    L0::Real
    L1::Real
    f0::Vector{Float64}
    f1::Vector{Float64}
    m::Int64
    pgrid::LinSpace{Float64}
    sol::WFSolution
end

function WaldFriedman(c, L0, L1, f0, f1; m=25)
    pgrid = linspace(0.0, 1.0, m)

    # Renormalize distributions so nothing is "too" small
    f0 = clamp(f0, 1e-8, 1-1e-8)
    f1 = clamp(f1, 1e-8, 1-1e-8)
    f0 = f0 / sum(f0)
    f1 = f1 / sum(f1)
    J = zeros(m)
    lb = 0.
    ub = 0.

    WaldFriedman(c, L0, L1, f0, f1, m, pgrid, WFSolution(J, lb, ub))
end

"""
This function takes a value for the probability with which
the correct model is model 0 and returns the mixed
distribution that corresponds with that belief.
"""
function current_distribution(wf::WaldFriedman, p)
    return p*wf.f0 + (1-p)*wf.f1
end

"""
This function takes a value for p, and a realization of the
random variable and calculates the value for p tomorrow.
"""
function bayes_update_k(wf::WaldFriedman, p, k)
    f0_k = wf.f0[k]
    f1_k = wf.f1[k]

    p_tp1 = p*f0_k / (p*f0_k + (1-p)*f1_k)

    return clamp(p_tp1, 0, 1)
end

"""
This is similar to `bayes_update_k` except it returns a
new value for p for each realization of the random variable
"""
function bayes_update_all(wf::WaldFriedman, p)
    return clamp(p*wf.f0 ./ (p*wf.f0 + (1-p)*wf.f1), 0, 1)
end

"""
For a given probability specify the cost of accepting model 0
"""
function payoff_choose_f0(wf::WaldFriedman, p)
    return (1-p)*wf.L0
end

"""
For a given probability specify the cost of accepting model 1
"""
function payoff_choose_f1(wf::WaldFriedman, p)
    return p*wf.L1
end

"""
This function evaluates the expectation of the value function
at period t+1. It does so by taking the current probability
distribution over outcomes:

    p(z_{k+1}) = p_k f_0(z_{k+1}) + (1-p_k) f_1(z_{k+1})

and evaluating the value function at the possible states
tomorrow J(p_{t+1}) where

    p_{t+1} = p f0 / ( p f0 + (1-p) f1)

Parameters
----------
p : scalar
    The current believed probability that model 0 is the true
    model.
J : function (interpolant)
    The current value function for a decision to continue

Returns
-------
EJ : scalar
    The expected value of the value function tomorrow
"""
function EJ(wf::WaldFriedman, p, J)
    # Pull out information
    f0, f1 = wf.f0, wf.f1

    # Get the current believed distribution and tomorrows possible dists
    # Need to clip to make sure things don't blow up (go to infinity)
    curr_dist = current_distribution(wf, p)
    tp1_dist = bayes_update_all(wf, p)

    # Evaluate the expectation
    EJ = dot(curr_dist, J.(tp1_dist))

    return EJ
end

"""
For a given probability distribution and value function give
cost of continuing the search for correct model
"""
function payoff_continue(wf::WaldFriedman, p, J)
    return wf.c + EJ(wf, p, J)
end

"""
Evaluates the value function for a given continuation value
function; that is, evaluates

    J(p) = min( (1-p)L0, pL1, c + E[J(p')])

Uses linear interpolation between points
"""
function bellman_operator(wf::WaldFriedman, J)
    c, L0, L1, f0, f1 = wf.c, wf.L0, wf.L1, wf.f0, wf.f1
    m, pgrid = wf.m, wf.pgrid

    J_out = zeros(m)
    J_interp = LinInterp(pgrid, J)

    for (p_ind, p) in enumerate(pgrid)
        # Payoff of choosing model 0
        p_c_0 = payoff_choose_f0(wf, p)
        p_c_1 = payoff_choose_f1(wf, p)
        p_con = payoff_continue(wf, p, J_interp)

        J_out[p_ind] = min(p_c_0, p_c_1, p_con)
    end

    return J_out
end

"""
This function takes a value function and returns the corresponding
cutoffs of where you transition between continue and choosing a
specific model
"""
function find_cutoff_rule(wf::WaldFriedman, J)
    m, pgrid = wf.m, wf.pgrid

    # Evaluate cost at all points on grid for choosing a model
    p_c_0 = payoff_choose_f0(wf, pgrid)
    p_c_1 = payoff_choose_f1(wf, pgrid)

    # The cutoff points can be found by differencing these costs with
    # the Bellman equation (J is always less than or equal to p_c_i)
    lb = pgrid[searchsortedlast(p_c_1 - J, 1e-10)]
    ub = pgrid[searchsortedlast(J - p_c_0, -1e-10)]

    return lb, ub
end

function solve_model(wf; tol=1e-7)
    bell_op(x) = bellman_operator(wf, x)
    J =  compute_fixed_point(bell_op, zeros(wf.m), err_tol=tol, print_skip=5)

    wf.sol.J = J
    wf.sol.lb, wf.sol.ub = find_cutoff_rule(wf, J)
    return J
end

"""
This function takes an initial condition and simulates until it
stops (when a decision is made).
"""
function simulate(wf::WaldFriedman, f; p0=0.5)
    # Check whether vf is computed
    if sumabs(wf.sol.J) < 1e-8
        solve_model(wf)
    end

    # Unpack useful info
    lb, ub = wf.sol.lb, wf.sol.ub
    drv = DiscreteRV(f)

    # Initialize a couple useful variables
    decision_made = false
    decision = 0
    p = p0
    t = 0

    while !decision_made
        # Maybe should specify which distribution is correct one so that
        # the draws come from the "right" distribution
        k = rand(drv)
        t = t+1
        p = bayes_update_k(wf, p, k)
        if p < lb
            decision_made = true
            decision = 1
        elseif p > ub
            decision_made = true
            decision = 0
        end
    end

    return decision, p, t
end

"""
Uses the distribution f0 as the true data generating
process
"""
function simulate_tdgp_f0(wf::WaldFriedman; p0=0.5)
    decision, p, t = simulate(wf, wf.f0; p0=p0)

    if decision == 0
        correct = true
    else
        correct = false
    end

    return correct, p, t
end

"""
Uses the distribution f1 as the true data generating
process
"""
function simulate_tdgp_f1(wf::WaldFriedman; p0=0.5)
    decision, p, t = simulate(wf, wf.f1; p0=p0)

    if decision == 1
        correct = true
    else
        correct = false
    end

    return correct, p, t
end

"""
Simulates repeatedly to get distributions of time needed to make a
decision and how often they are correct.
"""
function stopping_dist(wf::WaldFriedman; ndraws=250, tdgp="f0")
    if tdgp=="f0"
        simfunc = simulate_tdgp_f0
    else
        simfunc = simulate_tdgp_f1
    end

    # Allocate space
    tdist = Array(Int64, ndraws)
    cdist = Array(Bool, ndraws)

    for i in 1:ndraws
        correct, p, t = simfunc(wf)
        tdist[i] = t
        cdist[i] = correct
    end

    return cdist, tdist
end


text: q=parse;r=reduce;v,b=0,join(reverse(digits(q(BigInt,"1"*readline(),base=16),base=2))[2:end]);function p(i);t(c)=q(Int,b[i[1]:((i[1]+=c)-1)],base=2);global v+=t(3);d,c=t(3),1;d==4&&return r((a,b)->16a+b,[t(4) for _=1:99 if c==1&&(c=t(1))<2]);a(v)=r([+,*,min,max,i,>,<,==][d+1],v);(t(1)==0&&(g=t(15)+i[1])!=0) ? (return[p(i) for _=1:99 if i[1]<g]|>a) : return[p(i) for _=1:t(11)]|>a;end;p([1])|>b->println("$v\n$b")

max_stars_repo_name: LiquidFun/adventofcode

max_issues_repo_name: LiquidFun/adventofcode

max_forks_repo_name: LiquidFun/adventofcode


text:
spg_get_major_version() = Base.convert(Int64,
    ccall((:spg_get_major_version, libsymspg), Cint, () ))
spg_get_minor_version() = Base.convert(Int64,
    ccall((:spg_get_minor_version, libsymspg), Cint, () ))
spg_get_micro_version() = Base.convert(Int64,
    ccall((:spg_get_micro_version, libsymspg), Cint, () ))


text:
module PropertyGraphMongo

using PropertyGraph
using Mongo
using LibBSON

using UUID

# package code goes here
export MongoStore
export MongoGraphLoader
export ChangeTrackerNotDefinedException

export persist

include("persistenceoperationresult.jl")
include("persistencebatchresult.jl")
include("mongostore.jl")
include("mongographloader.jl")
include("exceptions.jl")

end # module


text:
using MaskArrays

using Distributions

# Get a positive definite matrix for the covariance
z = randn(6,5)
Σ = z' * z

# Use that to build a distribution
d = MvNormal(Σ)

# Sample from it and "forget" some values
x = Vector{Union{Missing, Float64}}(undef, 5)
x .= rand(d)
x[[2,4]] .= missing

# Use the result to build a MaskArray
ma = maskarray(x)

using TupleVectors
using Sobol
using StatsFuns

# Now sample and log-weight the result
ω = SobolHypercube(2)
s = @with (;ω) 10000 begin
    # A very naive proposal
    x = 5 * norminvcdf(rand(ω))
    y = 5 * norminvcdf(rand(ω))
    imputed(ma) .= [x,y]

    # The log-weight
    ℓ = logpdf(d, ma)

    (;ℓ, x, y)
end

using StatsBase

# Use the log-weights to resample
rs = TupleVector(sample(s, Weights(exp.(s.ℓ)), 1000))

# And plot the result
using UnicodePlots
@with TupleVectors.unwrap(rs) begin
    scatterplot(x, y)
end

# julia> @with TupleVectors.unwrap(rs) begin
#            scatterplot(x, y)
#        end
#       ┌────────────────────────────────────────┐
#     4 │⠀⠀⠀⠁⠀⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
#       │⠀⠀⠀⠀⠀⠀⢀⠀⢀⡀⣀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
#       │⠀⠀⠀⠀⠀⠀⠀⠀⢀⠀⠀⢄⢀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
#       │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⢂⡱⣃⡀⡃⡄⡄⢼⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
#       │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠚⠡⢒⢽⣿⢯⣸⠄⠠⠀⠀⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
#       │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠄⠳⣽⣻⢿⢢⣷⠍⣂⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
#       │⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠥⢷⡿⣿⢷⡧⡿⣼⣮⣤⡤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤│
#       │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢼⡷⣳⠵⣿⢞⣚⢆⠧⠄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
#       │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠁⠈⣛⠾⣯⢝⢣⢽⡲⢇⡠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
#       │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠘⠘⠄⠟⢏⡧⢏⣷⢒⢦⠠⠀⠀⠀⠀⠀⠀⠀⠀⠀│
#       │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠊⠒⠫⣎⣿⢜⠪⠄⠂⢄⠀⠀⠀⠀⠀⠀│
#       │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠗⠌⠗⢆⠱⠂⠀⠀⠀⠀⠀⠀│
#       │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠣⠐⠡⡂⠀⠄⠁⠀⠀⠀⠀│
#       │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⠀⠀⠀│
#    -5 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠀⠀⠀⠀│
#       └────────────────────────────────────────┘
#       -4                                       5


text:
import CompositeTypes: components
import CompositeTypes.Display: displaysymbol

struct CompositeObject
    components
end
CompositeObject(components...) = CompositeObject(components)
components(obj::CompositeObject) = obj.components
displaysymbol(obj::CompositeObject) = 'C'
Display.object_parentheses(::CompositeObject) = false
Display.displaystencil(obj::CompositeObject) = composite_displaystencil(obj)
Base.show(io::IO, mime::MIME"text/plain", obj::CompositeObject) =
    composite_show(io, mime, obj)

struct StencilObject
    map
    object
end
Base.show(io::IO, mime::MIME"text/plain", obj::StencilObject) =
    composite_show(io, mime, obj)
Display.displaystencil(obj::StencilObject) =
    [Display.SymbolObject(obj.map), '(', obj.object, ')']

struct SetUnion
    components
end
components(obj::SetUnion) = obj.components
Display.combinationsymbol(object::SetUnion) = Display.Symbol('∪')
Display.displaystencil(obj::SetUnion) = composite_displaystencil(obj)

Base.show(io::IO, mime::MIME"text/plain", obj::SetUnion) =
    composite_show(io, mime, obj)


text:
# Basic Lead/Lag:

# syntax:
lead_lag_docstring(op, opped, prevnext, firstlast) = """\n
        $op(term, nsteps::Integer)

    This `@formula` term is used to introduce $opped variables.
    For example `$op(x,1)` effectively adds a new column containing
    the value of the `x` column from the $prevnext row.
    If there is no such row (e.g. because this is the $firstlast row),
    then the $opped column will contain `missing` for that entry.

    Note: this is only a basic row-wise $op operation.
    It is up to the user to ensure that data is sorted by the temporal variable,
    and that observations are spaced with regular time-steps.
    (Which may require adding extra-rows filled with `missing` values.)
"""

@doc lead_lag_docstring("lag", "lagged", "previous", "first") lag
@doc lead_lag_docstring("lead", "lead", "next", "last") lead

# struct for behavior
struct LeadLagTerm{T<:AbstractTerm, F<:Union{typeof(lead), typeof(lag)}} <: AbstractTerm
    term::T
    nsteps::Int
end

function apply_schema(t::FunctionTerm{F}, sch::Schema, ctx::Type) where F<:Union{typeof(lead), typeof(lag)}
    opname = string(nameof(F.instance))
    if length(t.args_parsed) == 1  # lag(term)
        term_parsed = first(t.args_parsed)
        nsteps = 1
    elseif length(t.args_parsed) == 2  # lag(term, nsteps)
        term_parsed, nsteps_parsed = t.args_parsed
        (nsteps_parsed isa ConstantTerm) ||
            throw(ArgumentError("$opname step must be a number (got $nsteps_parsed)"))
        nsteps = nsteps_parsed.n
    else
        throw(ArgumentError("`$opname` terms require 1 or 2 arguments."))
    end

    term = apply_schema(term_parsed, sch, ctx)
    return LeadLagTerm{typeof(term), F}(term, nsteps)
end

function modelcols(ll::LeadLagTerm{<:Any, F}, d::Tables.ColumnTable) where F
    original_cols = modelcols(ll.term, d)
    return F.instance(original_cols, ll.nsteps)
end

width(ll::LeadLagTerm) = width(ll.term)
function Base.show(io::IO, ll::LeadLagTerm{<:Any, F}) where F
    opname = string(nameof(F.instance))
    print(io, "$opname($(ll.term), $(ll.nsteps))")
end
function StatsBase.coefnames(ll::LeadLagTerm{<:Any, F}) where F
    opname = string(nameof(F.instance))
    coefnames(ll.term) .* "_$opname$(ll.nsteps)"
end


text:
"An AbstractGeometryTrait type for all geometries."
abstract type AbstractGeometryTrait end

"An AbstractGeometryCollectionTrait type for all geometrycollections."
abstract type AbstractGeometryCollectionTrait <: AbstractGeometryTrait end
"A GeometryCollection is a collection of `Geometry`s."
struct GeometryCollectionTrait <: AbstractGeometryCollectionTrait end

"An AbstractPointTrait for all points."
abstract type AbstractPointTrait <: AbstractGeometryTrait end
"A single point."
struct PointTrait <: AbstractPointTrait end

"An AbstractCurveTrait type for all curves."
abstract type AbstractCurveTrait <: AbstractGeometryTrait end
"An AbstractLineString type for all linestrings."
abstract type AbstractLineStringTrait <: AbstractCurveTrait end
"A LineStringTrait is a collection of straight lines between its `PointTrait`s."
struct LineStringTrait <: AbstractLineStringTrait end
"A LineTrait is [`LineStringTrait`](@ref) with just two points."
struct LineTrait <: AbstractLineStringTrait end
"A LinearRingTrait is a [`LineStringTrait`](@ref) with the same begin and endpoint."
struct LinearRingTrait <: AbstractLineStringTrait end

"A CircularStringTrait is a curve, with an odd number of points.
A single segment consists of three points, where the first and last are the beginning and end,
while the second is halfway the curve."
struct CircularStringTrait <: AbstractCurveTrait end
"A CompoundCurveTrait is a curve that combines straight [`LineStringTrait`](@ref)s and curved [`CircularStringTrait`](@ref)s."
struct CompoundCurveTrait <: AbstractCurveTrait end

"An AbstractSurfaceTrait type for all surfaces."
abstract type AbstractSurfaceTrait <: AbstractGeometryTrait end
"An AbstractCurvePolygonTrait type for all curved polygons."
abstract type AbstractCurvePolygonTrait <: AbstractSurfaceTrait end
"An [`AbstractCurvePolygonTrait`](@ref) that can contain either circular or straight curves as rings."
struct CurvePolygonTrait <: AbstractCurvePolygonTrait end
"An AbstractPolygonTrait type for all polygons."
abstract type AbstractPolygonTrait <: AbstractCurvePolygonTrait end
"An [`AbstractSurfaceTrait`](@ref) with straight rings either as exterior or interior(s)."
struct PolygonTrait <: AbstractPolygonTrait end
"A [`PolygonTrait`](@ref) that is triangular."
struct TriangleTrait <: AbstractPolygonTrait end
"A [`PolygonTrait`](@ref) that is rectangular and could be described by the minimum and maximum vertices."
struct RectangleTrait <: AbstractPolygonTrait end
"A [`PolygonTrait`](@ref) with four vertices."
struct QuadTrait <: AbstractPolygonTrait end
"A [`PolygonTrait`](@ref) with five vertices."
struct PentagonTrait <: AbstractPolygonTrait end
"A [`PolygonTrait`](@ref) with six vertices."
struct HexagonTrait <: AbstractPolygonTrait end

"An AbstractPolyhedralSurfaceTrait type for all polyhedralsurfaces."
abstract type AbstractPolyhedralSurfaceTrait <: AbstractSurfaceTrait end
"A PolyhedralSurfaceTrait is a connected surface consisting of [`PolygonTrait`](@ref)s."
struct PolyhedralSurfaceTrait <: AbstractPolyhedralSurfaceTrait end
"A TINTrait is a [`PolyhedralSurfaceTrait`](@ref) consisting of [`TriangleTrait`](@ref)s."
struct TINTrait <: AbstractPolyhedralSurfaceTrait end  # Surface consisting of Triangles

"An AbstractMultiPointTrait type for all multipoints."
abstract type AbstractMultiPointTrait <: AbstractGeometryCollectionTrait end
"A MultiPointTrait is a collection of [`PointTrait`](@ref)s."
struct MultiPointTrait <: AbstractMultiPointTrait end

"An AbstractMultiCurveTrait type for all multicurves."
abstract type AbstractMultiCurveTrait <: AbstractGeometryCollectionTrait end
"A MultiCurveTrait is a collection of [`CircularStringTrait`](@ref)s."
struct MultiCurveTrait <: AbstractMultiCurveTrait end
"An AbstractMultiLineStringTrait type for all multilinestrings."
abstract type AbstractMultiLineStringTrait <: AbstractMultiCurveTrait end
"A MultiLineStringTrait is a collection of [`LineStringTrait`](@ref)s."
struct MultiLineStringTrait <: AbstractMultiLineStringTrait end

"An AbstractMultiSurfaceTrait type for all multisurfaces."
abstract type AbstractMultiSurfaceTrait <: AbstractGeometryCollectionTrait end
"A MultiSurfaceTrait is a collection of [`AbstractSurfaceTrait`](@ref)s."
struct MultiSurfaceTrait <: AbstractMultiSurfaceTrait end
"An AbstractMultiPolygonTrait type for all multipolygons."
abstract type AbstractMultiPolygonTrait <: AbstractMultiSurfaceTrait end
"A MultiPolygonTrait is a collection of [`PolygonTrait`](@ref)s."
struct MultiPolygonTrait <: AbstractMultiPolygonTrait end

"An AbstractFeatureTrait for all features"
abstract type AbstractFeatureTrait end
"A FeatureTrait holds `geometries`, `properties` and an `extent`"
struct FeatureTrait <: AbstractFeatureTrait end

"An AbstractFeatureCollectionTrait for all feature collections"
abstract type AbstractFeatureCollectionTrait end
"A FeatureCollectionTrait holds objects of `FeatureTrait` and an `extent`"
struct FeatureCollectionTrait <: AbstractFeatureCollectionTrait end


text:
# This file was generated, do not modify it. # hide
mtm.model.model.reg = LassoRegressor()
mtm.model.range = range(model, :(reg.lambda), lower=500, upper=100_000, scale=:log10)
fit!(mtm, rows=train)

best_mdl = fitted_params(mtm).best_model
round(best_mdl.reg.lambda, sigdigits=4)

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
"""
A transient version of the `Gridap` `FEOperator` that depends on time
"""
abstract type TransientFEOperator{C<:OperatorType} <: GridapType end

"""
Returns the test space
"""
function get_test(op::TransientFEOperator)
  @abstractmethod
end

"""
Returns the (possibly) time-dependent trial space
"""
function get_trial(op::TransientFEOperator)
  @abstractmethod # time dependent
end

function allocate_residual(op::TransientFEOperator,uh,cache)
  @abstractmethod
end

function allocate_jacobian(op::TransientFEOperator,uh,cache)
  @notimplemented
end

"""
Idem as `residual!` of `ODEOperator`
"""
function residual!(
  b::AbstractVector,
  op::TransientFEOperator,
  t::Real,
  xh::Union{AbstractVector,Tuple{Vararg{AbstractVector}}},
  cache)
  @abstractmethod
end

"""
Idem as `jacobian!` of `ODEOperator`
"""
function jacobian!(
  A::AbstractMatrix,
  op::TransientFEOperator,
  t::Real,
  xh::Union{AbstractVector,Tuple{Vararg{AbstractVector}}},
  i::Int,
  γᵢ::Real,
  cache)
  @abstractmethod
end

"""
Idem as `jacobians!` of `ODEOperator`
"""
function jacobians!(
  A::AbstractMatrix,
  op::TransientFEOperator,
  t::Real,
  x::Tuple{Vararg{AbstractVector}},
  γ::Tuple{Vararg{Real}},
  cache)
  @abstractmethod
end

"""
Returns the assembler, which is constant for all time steps for a given FE
operator.

Note: adaptive FE spaces involve to generate new FE spaces and
corresponding operators, due to the ummutable approach in `Gridap`
"""
get_assembler(feop::TransientFEOperator) = @abstractmethod

# Default API

"""
Returns a `ODEOperator` wrapper of the `TransientFEOperator` that can be
straightforwardly used with the `ODETools` module.
"""
function get_algebraic_operator(feop::TransientFEOperator{C}) where C
  ODEOpFromFEOp{C}(feop)
end

OperatorType(::Type{<:TransientFEOperator{C}}) where C = C

# @fverdugo This function is just in case we need to override it in the future for some specialization.
# This default implementation is enough for the moment.
function allocate_cache(op::TransientFEOperator)
  nothing
end

function update_cache!(cache::Nothing,op::TransientFEOperator,t::Real)
  nothing
end

# Specializations

"""
Transient FE operator that is defined by a transient Weak form
"""
struct TransientFEOperatorFromWeakForm{C} <: TransientFEOperator{C}
  res::Function
  jacs::Tuple{Vararg{Function}}
  assem_t::Assembler
  trials::Tuple{Vararg{Any}}
  test::FESpace
  order::Integer
end

function TransientConstantFEOperator(m::Function,a::Function,b::Function,
  trial,test)
  res(t,(u,ut),v) = -1.0 * b(v)
  jac(t,(u,ut),du,v) = a(du,v)
  jac_t(t,(u,ut),dut,v) = m(dut,v)
  assem_t = SparseMatrixAssembler(trial,test)
  TransientFEOperatorFromWeakForm{Constant}(res,(jac,jac_t),assem_t,(trial,∂t(trial)),test,1)
end

function TransientConstantMatrixFEOperator(m::Function,a::Function,b::Function,
  trial,test)
  res(t,(u,ut),v) = m(ut,v) + a(u,v) - b(t,v)
  jac(t,(u,ut),du,v) = a(du,v)
  jac_t(t,(u,ut),dut,v) = m(dut,v)
  assem_t = SparseMatrixAssembler(trial,test)
  TransientFEOperatorFromWeakForm{ConstantMatrix}(res,(jac,jac_t),assem_t,(trial,∂t(trial)),test,1)
end

function TransientAffineFEOperator(m::Function,a::Function,b::Function,
  trial,test)
  res(t,(u,ut),v) = m(t,ut,v) + a(t,u,v) - b(t,v)
  jac(t,(u,ut),du,v) = a(t,du,v)
  jac_t(t,(u,ut),dut,v) = m(t,dut,v)
  assem_t = SparseMatrixAssembler(trial,test)
  TransientFEOperatorFromWeakForm{Affine}(res,(jac,jac_t),assem_t,(trial,∂t(trial)),test,1)
end

function TransientFEOperator(res::Function,jac::Function,jac_t::Function,
  trial,test)
  assem_t = SparseMatrixAssembler(trial,test)
  TransientFEOperatorFromWeakForm{Nonlinear}(res,(jac,jac_t),assem_t,(trial,∂t(trial)),test,1)
end

function TransientConstantFEOperator(m::Function,c::Function,a::Function,b::Function,
  trial,test)
  res(t,(u,ut,utt),v) = -1.0 * b(v)
  jac(t,(u,ut,utt),du,v) = a(du,v)
  jac_t(t,(u,ut,utt),dut,v) = c(dut,v)
  jac_tt(t,(u,ut,utt),dutt,v) = m(dutt,v)
  assem_t = SparseMatrixAssembler(trial,test)
  trial_t = ∂t(trial)
  trial_tt = ∂t(trial_t)
  TransientFEOperatorFromWeakForm{Constant}(
    res,(jac,jac_t,jac_tt),assem_t,(trial,trial_t,trial_tt),test,2)
end

function TransientConstantMatrixFEOperator(m::Function,c::Function,a::Function,b::Function,
  trial,test)
  res(t,(u,ut,utt),v) = m(utt,v) + c(ut,v) + a(u,v) - b(t,v)
  jac(t,(u,ut,utt),du,v) = a(du,v)
  jac_t(t,(u,ut,utt),dut,v) = c(dut,v)
  jac_tt(t,(u,ut,utt),dutt,v) = m(dutt,v)
  assem_t = SparseMatrixAssembler(trial,test)
  trial_t = ∂t(trial)
  trial_tt = ∂t(trial_t)
  TransientFEOperatorFromWeakForm{ConstantMatrix}(
    res,(jac,jac_t,jac_tt),assem_t,(trial,trial_t,trial_tt),test,2)
end

function TransientAffineFEOperator(m::Function,c::Function,a::Function,b::Function,
  trial,test)
  res(t,(u,ut,utt),v) = m(t,utt,v) + c(t,ut,v) + a(t,u,v) - b(t,v)
  jac(t,(u,ut,utt),du,v) = a(t,du,v)
  jac_t(t,(u,ut,utt),dut,v) = c(t,dut,v)
  jac_tt(t,(u,ut,utt),dutt,v) = m(t,dutt,v)
  assem_t = SparseMatrixAssembler(trial,test)
  trial_t = ∂t(trial)
  trial_tt = ∂t(trial_t)
  TransientFEOperatorFromWeakForm{Affine}(
    res,(jac,jac_t,jac_tt),assem_t,(trial,trial_t,trial_tt),test,2)
end

function TransientFEOperator(res::Function,jac::Function,jac_t::Function,
  jac_tt::Function,trial,test)
  assem_t = SparseMatrixAssembler(trial,test)
  trial_t = ∂t(trial)
  trial_tt = ∂t(trial_t)
  TransientFEOperatorFromWeakForm{Nonlinear}(
    res,(jac,jac_t,jac_tt),assem_t,(trial,trial_t,trial_tt),test,2)
end

function TransientFEOperator(res::Function,trial,test;order::Integer=1)
  jacs = ()
  for i in 1:order+1
    function jac_i(t,x,dxi,dv)
      function res_i(y)
        x[i] = y
        res(t,x,dv)
      end
      jacobian(res_i,x[i])
    end
    jacs = (jacs...,jac_i)
  end
  # jac(t,u,ut,du,dv) = jacobian(x->res(t,x,ut,dv),u)
  # jac_t(t,u,ut,dut,dv) = jacobian(xt->res(t,u,xt,dv),ut)
  TransientFEOperator(res,jacs...,trial,test)
end

function SparseMatrixAssembler(
  trial::Union{TransientTrialFESpace,TransientMultiFieldTrialFESpace},
  test::FESpace)
  SparseMatrixAssembler(evaluate(trial,nothing),test)
end

get_assembler(op::TransientFEOperatorFromWeakForm) = op.assem_t

get_test(op::TransientFEOperatorFromWeakForm) = op.test

get_trial(op::TransientFEOperatorFromWeakForm) = op.trials[1]

get_order(op::TransientFEOperatorFromWeakForm) = op.order

function allocate_residual(op::TransientFEOperatorFromWeakForm,uh::FEFunction,cache)
  V = get_test(op)
  v = get_fe_basis(V)
  xh = ()
  for i in 1:get_order(op)+1
    xh = (xh...,uh)
  end
  vecdata = collect_cell_vector(V,op.res(0.0,xh,v))
  allocate_vector(op.assem_t,vecdata)
end

function residual!(
  b::AbstractVector,
  op::TransientFEOperatorFromWeakForm,
  t::Real,
  xh::Tuple{Vararg{FEFunction}},
  cache)
  V = get_test(op)
  v = get_fe_basis(V)
  vecdata = collect_cell_vector(V,op.res(t,xh,v))
  assemble_vector!(b,op.assem_t,vecdata)
  b
end

function allocate_jacobian(op::TransientFEOperatorFromWeakForm,uh::FEFunction,cache)
  xh = ()
  for i in 1:get_order(op)+1
    xh = (xh...,uh)
  end
  _matdata = ()
  for i in 1:get_order(op)+1
    _matdata = (_matdata...,matdata_jacobian(op,0.0,xh,i,0.0))
  end
  matdata = vcat_matdata(_matdata)
  allocate_matrix(op.assem_t,matdata)
end

function jacobian!(
  A::AbstractMatrix,
  op::TransientFEOperatorFromWeakForm,
  t::Real,
  xh::Tuple{Vararg{FEFunction}},
  i::Integer,
  γᵢ::Real,
  cache)
  matdata = matdata_jacobian(op,t,xh,i,γᵢ)
  assemble_matrix_add!(A,op.assem_t, matdata)
  A
end

function jacobians!(
  A::AbstractMatrix,
  op::TransientFEOperatorFromWeakForm,
  t::Real,
  xh::Tuple{Vararg{FEFunction}},
  γ::Tuple{Vararg{Real}},
  cache)
  _matdata = ()
  for i in 1:get_order(op)+1
    _matdata = (_matdata...,matdata_jacobian(op,t,xh,i,γ[i]))
  end
  matdata = vcat_matdata(_matdata)
  assemble_matrix_add!(A,op.assem_t, matdata)
  A
end

function vcat_matdata(_matdata)
  term_to_cellmat_j = ()
  term_to_cellidsrows_j = ()
  term_to_cellidscols_j = ()
  for j in 1:length(_matdata)
    term_to_cellmat_j = (term_to_cellmat_j...,_matdata[j][1])
    term_to_cellidsrows_j = (term_to_cellidsrows_j...,_matdata[j][2])
    term_to_cellidscols_j = (term_to_cellidscols_j...,_matdata[j][3])
  end

  term_to_cellmat = vcat(term_to_cellmat_j...)
  term_to_cellidsrows = vcat(term_to_cellidsrows_j...)
  term_to_cellidscols = vcat(term_to_cellidscols_j...)

  matdata = (term_to_cellmat,term_to_cellidsrows, term_to_cellidscols)
end

function matdata_jacobian(
  op::TransientFEOperatorFromWeakForm,
  t::Real,
  xh::Tuple{Vararg{FEFunction}},
  i::Integer,
  γᵢ::Real)
  Uh = evaluate(get_trial(op),nothing)
  V = get_test(op)
  du = get_trial_fe_basis(Uh)
  v = get_fe_basis(V)
  matdata = collect_cell_matrix(Uh,V,γᵢ*op.jacs[i](t,xh,du,v))
end

# Tester

function test_transient_fe_operator(op::TransientFEOperator,uh)
  odeop = get_algebraic_operator(op)
  @test isa(odeop,ODEOperator)
  cache = allocate_cache(op)
  V = get_test(op)
  @test isa(V,FESpace)
  U = get_trial(op)
  U0 = U(0.0)
  @test isa(U0,FESpace)
  r = allocate_residual(op,uh,cache)
  @test isa(r,AbstractVector)
  residual!(r,op,0.0,(uh,uh),cache)
  @test isa(r,AbstractVector)
  J = allocate_jacobian(op,uh,cache)
  @test isa(J,AbstractMatrix)
  jacobian!(J,op,0.0,(uh,uh),1,1.0,cache)
  @test isa(J,AbstractMatrix)
  jacobian!(J,op,0.0,(uh,uh),2,1.0,cache)
  @test isa(J,AbstractMatrix)
  jacobians!(J,op,0.0,(uh,uh),(1.0,1.0),cache)
  @test isa(J,AbstractMatrix)
  cache = update_cache!(cache,op,0.0)
  true
end


text:
@testset "Constructors" begin
    λ, γ, δ = 1.5, 2.5, 3.5

    # direct args constructor
    glr     = GeneralizedLinearRegression()
    ols     = LinearRegression()
    ridge   = RidgeRegression(λ; fit_intercept=false)
    lasso   = LassoRegression(λ; penalize_intercept=true)
    logreg0 = LogisticRegression(penalty=:none)
    logreg1 = LogisticRegression(λ)
    logreg2 = LogisticRegression(λ, γ)
    mnreg2  = MultinomialRegression(λ, γ)
    hlreg   = RobustRegression(HuberRho(δ), λ)
    ladreg  = LADRegression(λ)

    @test isa(glr.loss, L2Loss)
    @test isa(glr.penalty, NoPenalty)
    @test glr.fit_intercept
    @test !glr.penalize_intercept

    @test isa(ols.loss, L2Loss)
    @test isa(ols.penalty, NoPenalty)
    @test ols.fit_intercept

    @test isa(ridge.loss, L2Loss)
    @test isa(ridge.penalty, ScaledPenalty{L2Penalty})
    @test !ridge.fit_intercept
    @test ridge.penalty.scale == λ

    @test isa(lasso.loss, L2Loss)
    @test isa(lasso.penalty, ScaledPenalty{L1Penalty})
    @test lasso.penalty.scale == λ
    @test lasso.penalize_intercept

    @test isa(logreg0.loss, LogisticLoss)
    @test isa(logreg0.penalty, NoPenalty)

    @test isa(logreg1.loss, LogisticLoss)
    @test isa(logreg1.penalty, ScaledPenalty{L2Penalty})
    @test logreg1.penalty.scale == λ

    @test isa(logreg2.loss, LogisticLoss)
    @test isa(logreg2.penalty, CompositePenalty)
    @test isa(logreg2.penalty.penalties[1], ScaledPenalty{L2Penalty})
    @test isa(logreg2.penalty.penalties[2], ScaledPenalty{L1Penalty})
    @test logreg2.penalty.penalties[1].scale == λ
    @test logreg2.penalty.penalties[2].scale == γ

    @test isa(mnreg2.loss, MultinomialLoss)
    @test isa(mnreg2.penalty.penalties[2], ScaledPenalty{L1Penalty})
    @test mnreg2.penalty.penalties[1].scale == λ
    @test mnreg2.penalty.penalties[2].scale == γ

    @test isa(hlreg.loss, RobustLoss{HuberRho{δ}})
    @test isa(hlreg.penalty, ScaledPenalty{L2Penalty})
    @test hlreg.penalty.scale == λ
    @test hlreg.fit_intercept

    @test isa(ladreg.loss, RobustLoss{QuantileRho{0.5}})
    @test isa(ladreg.penalty, ScaledPenalty{L2Penalty})
    @test ladreg.fit_intercept

    # ======

    ridge   = RidgeRegression(lambda=λ, fit_intercept=true)
    lasso   = LassoRegression(lambda=λ)
    logreg1 = LogisticRegression(lambda=λ)
    logreg2 = LogisticRegression(lambda=λ, gamma=γ)
    mnreg2  = MultinomialRegression(lambda=λ, gamma=γ)
    hlreg   = RobustRegression(rho=HuberRho(delta=δ), lambda=λ)
    qreg    = QuantileRegression(delta=δ, lambda=λ)

    @test ridge.penalty.scale == λ
    @test ridge.fit_intercept
    @test lasso.penalty.scale == λ
    @test logreg1.penalty.scale == λ
    @test logreg2.penalty.penalties[1].scale == λ
    @test logreg2.penalty.penalties[2].scale == γ
    @test mnreg2.penalty.penalties[1].scale == λ
    @test mnreg2.penalty.penalties[2].scale == γ
    @test hlreg.loss.rho isa HuberRho{δ}
    @test hlreg.penalty.scale == λ
    @test qreg.loss.rho isa QuantileRho{δ}
    @test qreg.penalty.scale == λ
end


text:
using SpikingNeuralNetworks
using ClearStacktrace
using Plots
unicodeplots()
SNN = SpikingNeuralNetworks
SNN.@load_units
#include("../current_search.jl")
using SpikeNetOpt
SNO = SpikeNetOpt
@testset "IZHI_spike_search" begin

    cell_type = "IZHI"
    ngt_spikes=rand(5:15)
    RS = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -65, d = 8))

    current_ = SNO.current_search(cell_type,RS,ngt_spikes)

    RS.I = [current_*nA]
    SNN.monitor(RS, [:v])
    SNN.monitor(RS, [:fire])
    #test_result(nspk, ngt_spikes, 1e-1)

    ALLEN_DURATION = 2000 * ms
    ALLEN_DELAY = 1000 * ms

    SNN.sim!([RS]; dt =1*ms, delay=ALLEN_DELAY,stimulus_duration=ALLEN_DURATION,simulation_duration = ALLEN_DURATION+ALLEN_DELAY+443ms)

    spikes = SNO.get_spikes(RS)
    spikes = [s*ms for s in spikes]
    nspk = size(spikes)[1]
    @test nspk==ngt_spikes

    #v = SNN.vecplot(RS, :v)
    #v |> display
end

#=
@testset "ADEXP_spike_search" begin

    cell_type = "ADEXP"
    ngt_spikes=rand(4:20)

    adparam = SNN.ADEXParameter(;a = 6.050246708405076, b = 7.308480222357973,
        cm = 803.1019662706587,
        v0= -63.22881649139353,
        τ_m=19.73777028610565,
        τ_w=351.0551915202058,
        θ=-39.232165554444265,
        delta_T=6.37124632135508,
        v_reset = -59.18792270568965,
        spike_delta = 16.33506432689027)

    E = SNN.AD(;N = 1, param=adparam)
    current_ = SNO.current_search(cell_type,E,ngt_spikes)
    E.I = [current_*nA]
    SNN.monitor(E, [:v])
    SNN.monitor(E, [:fire])

    ALLEN_DURATION = 2000 * ms
    ALLEN_DELAY = 1000 * ms
    SNN.sim!([E]; dt =1*ms, delay=ALLEN_DELAY,stimulus_duration=ALLEN_DURATION,simulation_duration = ALLEN_DURATION+ALLEN_DELAY+443ms)
    spikes = SNO.raster_synchp(E)
    spikes = [s*ms for s in spikes]
    nspk = size(spikes)[1]
    #test_result(nspk, ngt_spikes, 1e-1)
    @test nspk==ngt_spikes

    #v = SNN.vecplot(E, :v)
    #v |> display
end
=#
#=
@testset "IZHI" begin

    RS = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -65, d = 8))
    IB = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -55, d = 4))
    CH = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -50, d = 2))
    FS = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.1, b = 0.2, c = -65, d = 2))
    TC1 = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.25, c = -65, d = 0.05))
    TC2 = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.25, c = -65, d = 0.05))
    RZ = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.1, b = 0.26, c = -65, d = 2))
    LTS = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.1, b = 0.25, c = -65, d = 2))
    P = [RS, IB, CH, FS, TC1, TC2, RZ, LTS]

    SNN.monitor(P, [:v])
    SNN.monitor(P, [:fire])

    T = 2second
    for t = 0:T
        for p in [RS, IB, CH, FS, LTS]
            p.I = [10]
        end
        TC1.I = [(t < 0.2T) ? 0mV : 2mV]
        TC2.I = [(t < 0.2T) ? -30mV : 0mV]
        RZ.I =  [(0.5T < t < 0.6T) ? 10mV : 0mV]
        SNN.sim!(P, [], 0.1ms)

        #test_result(nspk, ngt_spikes, 1e-1)

        #v = vecplot(E, :v)
        #@show(v)
    end
    for p in P
        spikes = raster_synchp(p)
        spikes = [s*ms for s in spikes]
        nspk = size(spikes)[1]
        @test nspk>=1

        v = SNN.vecplot(p, :v)
        v |> display
    end
end
=#


text:
using Tofu
using Test

@testset "Tofu.jl" begin
    @test (◻.x)((x=1,)) == 1
    @test (◻.x + 2)((x=1,)) == 3
    @test_broken detect_ambiguities(Tofu, Base, imported=true, recursive=true) == []
end


text:
module OodleAnalytics

include("plotting/Plotting.jl")

end


text:
using PkgSRA, Test

@testset "sindy_turing_test.jl" begin
    include("sindy_turing_test.jl")
end
@testset "sindy_test.jl" begin
    include("sindy_test.jl")
end
@testset "regression_utils_test.jl" begin
    include("regression_utils_test.jl")
end