text:
using SphericalHarmonicArrays
import SphericalHarmonicArrays: ModeMismatchError
import SphericalHarmonicArrays: modes, shdims, shmodes

using OffsetArrays
using SphericalHarmonicModes
using Test
using LinearAlgebra
using Aqua

@testset "project quality" begin
    if VERSION >= v"1.6.0"
        Aqua.test_all(SphericalHarmonicArrays, unbound_args = false)
    else
        Aqua.test_all(SphericalHarmonicArrays, unbound_args = false, ambiguities = false)
    end
end

@testset "Constructors" begin
    mode_lm = LM(0:1, 0:0)
    mode_ml = ML(0:1, 0:0)
    mode_l₂l₁ = L2L1Triangle(0:2, 2,0:2)
    @testset "0dim" begin
        @test SHArray(zeros(), ()) == zeros()
        S = SHArray(zeros())
        @test S == zeros()
        @test deepcopy(S) == S
        @test deepcopy(S) !== S
    end
    function checkarrs(sa, arr::AbstractArray)
        @test sa == arr
        @test parent(sa) === arr
        @test SphericalHarmonicArrays.parenttype(sa) == typeof(arr)
        @test axes(sa) == axes(arr)
    end
    @testset "non SH" begin
        arr = zeros(ComplexF64, length(mode_lm))
        sa = SHArray(arr, axes(arr))
        checkarrs(sa, arr)
        sa = SHArray(arr, (:,))
        checkarrs(sa, arr)
        sad = deepcopy(sa)
        @test sad == sa
        @test modes(sad) == modes(sa)
        @test sad !== sa

        sa = SHArray(arr)
        checkarrs(sa, arr)
    end
    @testset "one SH axis" begin
        arr = zeros(ComplexF64, length(mode_lm))
        @testset "SHArray" begin
            sa = SHArray{ComplexF64}((mode_lm,))
            @test sa.modes == (mode_lm,)

            sa = SHArray(arr,(mode_lm,))
            @test sa.modes == (mode_lm,)
            sa = SHArray(arr, mode_lm)
            @test sa.modes == (mode_lm,)

            sa2 = SHArray(sa, (mode_lm,))
            @test parent(sa2) === parent(sa)
            @test sa2.modes == (mode_lm,)

            sa2 = SHArray(sa, mode_lm)
            @test parent(sa2) === parent(sa)
            @test sa2.modes == (mode_lm,)

            sa2d = deepcopy(sa2)
            @test sa2d == sa2
            @test modes(sa2d) == modes(sa2)
            @test sa2d !== sa2
        end
        @testset "SHVector" begin
            @test SHVector{ComplexF64}(mode_lm) == arr
            @test SHVector{Float64}(mode_lm) == real(arr)
            @test SHVector(arr, mode_lm) == arr

            sv = SHVector(arr, mode_lm)
            svundef = SHVector{eltype(arr)}(undef, mode_lm)
            @test svundef.modes == sv.modes
            @test axes(parent(svundef)) == axes(parent(sv))

            svd = deepcopy(sv)
            @test svd == sv
            @test modes(svd) == modes(sv)
            @test svd !== sv
        end
    end
    @testset "2D mixed axes" begin
        @testset "first SH" begin
            arr = zeros(ComplexF64, length(mode_lm),2)
            ax = (mode_lm, 1:2)

            function testSA(sa, arr, ax)
                @test parent(sa) == arr
                @test sa.modes == ax
            end

            sa = SHArray{ComplexF64}(ax)
            testSA(sa, arr, ax)

            sa = SHArray(arr, ax)
            testSA(sa, arr, ax)
            @test parent(sa) === arr

            sa = SHArray(arr, (mode_lm, :))
            testSA(sa, arr, ax)
            @test parent(sa) === arr
            sad = deepcopy(sa)
            @test sad == sa
            @test modes(sad) == modes(sa)
            @test sad !== sa

            sa = SHArray(arr, 1=>mode_lm)
            testSA(sa, arr, ax)
            @test parent(sa) === arr

            @test SHArray{ComplexF64}(ax) == sa
            @test SHArray{Float64}(ax) == SHArray(real(arr), ax)
            @test SHArray(arr, ax) == sa

            saundef = SHArray{ComplexF64}(undef, ax)
            @test saundef.modes == sa.modes
            @test axes(parent(saundef)) == axes(parent(sa))

            sa = SHArray{ComplexF64}((mode_lm, 2))
            @test axes(sa) == (1:length(mode_lm), 1:2)
        end
        @testset "second SH" begin
            arr = zeros(ComplexF64, 2,length(mode_lm),)
            ax = (1:2, mode_lm)
            @test SHArray{ComplexF64}(ax) == arr
            @test SHArray{Float64}(ax) == real(arr)
            @test SHArray(arr, ax) == arr

            @testset "Errors" begin
                @test_throws ArgumentError SHArray(arr,(1:2, LM(l_range(mode_lm))))
                @test_throws ArgumentError SHArray(arr,(1:3, LM(l_range(mode_lm))))
            end
        end
        @testset "both SH" begin
            arr = zeros(ComplexF64, length(mode_lm), length(mode_lm))
            ax = (mode_lm, mode_lm)
            @test SHArray{ComplexF64}(ax) == arr
            @test SHArray{Float64}(ax) == real(arr)
            @test SHArray(arr, ax) == arr

            ax = (mode_lm, mode_ml)
            sa = SHArray(arr, 1 => mode_lm, 2 => mode_ml)
            @test sa.modes === ax
            @test sa.parent === arr

            @testset "Errors" begin
                ax = (LM(l_range(mode_lm)),LM(l_range(mode_lm)))
                @test_throws ArgumentError SHArray(arr, ax)
            end
        end
        @testset "SHMatrix" begin
            function testSHMatrix(T, s, arr, mode1, mode2)
                @test s isa SHMatrix{T}
                @test s == arr
                @test parent(s) === arr
                @test shdims(s) == (1, 2)
                @test modes(s) == (mode1, mode2)
            end
            function testSHMatrix(T, s, mode1, mode2)
                @test s isa SHMatrix{T}
                @test shdims(s) == (1, 2)
                @test modes(s) == (mode1, mode2)
            end
            function testSHMatrix(mode1, mode2)
                arr = zeros(ComplexF64, length(mode1),length(mode2))

                sm1 = SHMatrix(arr, mode1, mode2)
                testSHMatrix(ComplexF64, sm1, arr, mode1, mode2)

                sm2 = SHMatrix(arr, (mode1, mode2))
                testSHMatrix(ComplexF64, sm2, arr, mode1, mode2)

                @test sm1 == sm2

                sm3 = SHMatrix{ComplexF64}(mode1, mode2)
                testSHMatrix(ComplexF64, sm3, mode1, mode2)
                @test sm3 == arr
                @test all(iszero, sm3)

                sm4 = SHMatrix{ComplexF64}((mode1, mode2))
                testSHMatrix(ComplexF64, sm4, mode1, mode2)
                @test sm4 == sm3
                @test all(iszero, sm4)

                sm5 = SHMatrix{Float64}(mode1, mode2)
                testSHMatrix(Float64, sm5, mode1, mode2)
                @test sm5 == real(arr)
                @test all(iszero, sm5)

                sm6 = SHMatrix{Float64}((mode1, mode2))
                testSHMatrix(Float64, sm6, mode1, mode2)
                @test sm6 == sm5
                @test all(iszero, sm6)
            end
            @testset "LM LM" begin
                testSHMatrix(mode_lm, mode_lm)
            end
            @testset "LM ML" begin
                testSHMatrix(mode_lm, mode_ml)
            end
            @testset "ML LM" begin
                testSHMatrix(mode_ml, mode_lm)
            end
            @testset "ML ML" begin
                testSHMatrix(mode_ml, mode_ml)
            end
            @testset "LM L2L1Triangle" begin
                testSHMatrix(mode_lm, mode_l₂l₁)
            end
            @testset "ML L2L1Triangle" begin
                testSHMatrix(mode_ml, mode_l₂l₁)
            end
            @testset "L2L1Triangle LM" begin
                testSHMatrix(mode_l₂l₁,mode_lm)
            end
            @testset "L2L1Triangle ML" begin
                testSHMatrix(mode_l₂l₁,mode_ml)
            end
            @testset "L2L1Triangle L2L1Triangle" begin
                testSHMatrix(mode_l₂l₁,mode_l₂l₁)
            end

            sm = SHMatrix{ComplexF64}((mode_lm, mode_lm))
            smundef = SHMatrix{ComplexF64}(undef, (mode_lm, mode_lm))
            @test smundef.modes == sm.modes
            @test axes(parent(smundef)) == axes(parent(sm))
        end
        @testset "Pairs" begin
            mode1 = LM(1:1, 0:1)
            mode2 = LM(2:2, 0:2)
            mode3 = LM(2:2)

            s = SHArray(zeros(map(length, (mode1, mode2))), 1=>mode1, 2=>mode2)
            @test s.modes == (mode1, mode2)

            s = SHArray(zeros(map(length, (mode1, mode2))), 2=>mode2, 1=>mode1)
            @test s.modes == (mode1, mode2)

            @test_throws ArgumentError SHArray(zeros(map(length, (mode1, mode2))), 2=>mode2, 3=>mode1)

            # type-unstable version
            s = SHArray(zeros(map(length, (mode1, mode3))), 1=>mode1, 2=>mode3)
            @test s.modes == (mode1, mode3)
        end
    end
    @testset "2D none SH" begin
        sa = SHArray{ComplexF64}((2, 2))
        @test axes(sa) == (1:2, 1:2)
        @test sa.modes == (1:2, 1:2)
    end
    @testset "ArrayInitializers" begin
        @testset "undef" begin
            @testset "SHVector" begin
                sha = SHVector{Vector}(undef,(mode_lm,))
                @test all(i -> !isassigned(sha, i), eachindex(sha))
                @test size(sha) == (length(mode_lm),)
                sha = SHVector{Vector}(undef, mode_lm)
                @test all(i -> !isassigned(sha, i), eachindex(sha))
            end
            @testset "SHMatrix" begin
               sha = SHMatrix{Vector}(undef,(mode_lm, mode_ml))
               @test all(i -> !isassigned(sha, i), eachindex(sha))
               @test size(sha) == (length(mode_lm),length(mode_ml))
               sha = SHMatrix{Vector}(undef, mode_lm, mode_ml)
               @test all(i -> !isassigned(sha, i), eachindex(sha))
            end
            @testset "SHArray" begin
               sha = SHArray{Vector,2}(undef,(1:2, mode_ml))
               @test all(i -> !isassigned(sha, i), eachindex(sha))
               @test size(sha) == (2, length(mode_ml))

               sha = SHArray{Vector}(undef,(1:2, mode_ml))
               @test all(i -> !isassigned(sha, i), eachindex(sha))
               @test size(sha) == (2, length(mode_ml))

               sha = SHArray{Vector,2}(undef, 2, mode_ml)
               @test all(i -> !isassigned(sha, i), eachindex(sha))
               @test size(sha) == (2, length(mode_ml))

               sha = SHArray{Vector}(undef, 2, mode_ml)
               @test all(i -> !isassigned(sha, i), eachindex(sha))
               @test size(sha) == (2, length(mode_ml))
            end
        end
        @testset "missing" begin
            sha = SHArray{Union{Vector, Missing}}(missing, 2, mode_ml)
            @test all(i -> isassigned(sha, i), eachindex(sha))
            @test size(sha) == (2, length(mode_ml))

            sha = SHArray{Union{Vector, Missing},2}(missing, 2, mode_ml)
            @test all(i -> isassigned(sha, i), eachindex(sha))
            @test size(sha) == (2, length(mode_ml))
        end
    end
    @testset "convert" begin
        A = SHArray(zeros(2), LM(1, ZeroTo(1)))
        T = SHArray{Float64, 1, Vector{Float64}, Tuple{LM{SingleValuedRange, ZeroTo{true}}}}
        @test convert(T, A) isa T
        @test typeof(A)(A) === A

        T = SHArray{Float64, 1, Vector{Float64}, typeof(axes(A))}
        @test convert(T, A) isa T

        A = ones(2)
        T = SHArray{Float64, 1, Vector{Float64}, typeof(axes(A))}
        @test convert(T, A) isa T
    end
end

@testset "Base functions" begin
    mode_lm = LM(0:1, 0:0)
    mode_ml = ML(0:1, 0:0)

    arr = zeros(length(mode_lm),length(mode_ml))
    sm = SHMatrix(arr,(mode_lm, mode_ml))

    @testset "parent" begin
        @test parent(sm) == arr
        @test parent(sm) === arr
    end
    @testset "axes" begin
        @test axes(sm, 1) == axes(arr, 1)
        @test axes(sm, 2) == axes(arr, 2)
    end
    @testset "size" begin
        @test size(sm, 1) == size(arr, 1)
        @test size(sm, 2) == size(arr, 2)
    end
    @testset "dataids" begin
        @test Base.dataids(sm) == Base.dataids(arr)
    end

    @testset "misc" begin
        @test strides(sm) == strides(arr)
        @test Base.elsize(typeof(sm)) == Base.elsize(typeof(arr))
        @test Base.unsafe_convert(Ptr{eltype(sm)}, sm) == Base.unsafe_convert(Ptr{eltype(arr)}, arr)
    end

    @testset "zeros and ones" begin
        for (DT, f) in ((:zeros, :zero), (:ones, :one))
            @eval begin
                mode_lm = LM(0:1, 0:0)

                sa = $DT(mode_lm, mode_lm)
                @test sa.modes == (mode_lm, mode_lm)
                @test eltype(sa) == Float64
                @test all(==($f(Float64)), sa)

                sa = $DT((mode_lm, mode_lm))
                @test sa.modes == (mode_lm, mode_lm)
                @test eltype(sa) == Float64
                @test all(==($f(Float64)), sa)

                sa = $DT(mode_lm, 2)
                @test sa.modes == (mode_lm, 1:2)
                @test eltype(sa) == Float64
                @test all(==($f(Float64)), sa)

                sa = $DT(ComplexF64, mode_lm, mode_lm)
                @test sa.modes == (mode_lm, mode_lm)
                @test eltype(sa) == ComplexF64
                @test all(==($f(ComplexF64)), sa)

                sa = $DT(ComplexF64, mode_lm, 1:2)
                @test sa.modes == (mode_lm, 1:2)
                @test eltype(sa) == ComplexF64
                @test all(==($f(ComplexF64)), sa)
            end
        end
    end
end

@testset "Indexing" begin
    mode_lm = LM(0:1, 0:0)
    @testset "0D" begin
        s = SHArray(zeros())
        @test s[] == 0
        s[] = 3
        @test s[] == 3
    end
    @testset "1D SH" begin
        v = SHVector{ComplexF64}(mode_lm)
        for i in eachindex(parent(v))
            parent(v)[i] = i
        end
        @testset "getindex" begin
            @test to_indices(v, (1,)) == (1,)
            @test to_indices(v, (big(1),)) == (1,)
            @testset "linear" begin
                for i in eachindex(v)
                    @test v[i] == i
                    @test v[Int32(i)] == i
                    @test v[Int64(i)] == i
                end
                lininds = LinearIndices(v)
                for i in eachindex(IndexCartesian(), v)
                    @test v[i] == lininds[i]
                end
                @test v[CartesianIndices(v)] == v
                @test v[CartesianIndices(v), 1] == v
            end
            @testset "mode_lm" begin
                for (ind,m) in enumerate(mode_lm)
                    @test v[m] == ind
                end
            end
            @testset "colon" begin
                @test v[:] == 1:size(v,1)
            end
        end
        @testset "setindex" begin
            @testset "linear" begin
                for i in eachindex(v)
                    v[i] = i^2
                    @test v[i] == i^2
                end
                for (ind, I) in enumerate(eachindex(IndexCartesian(),v))
                    v[I] = ind
                    @test v[I] == ind
                end
            end
            @testset "mode_lm" begin
                for (i, m) in enumerate(mode_lm)
                    v[m] = i
                    @test v[m] == i
                end
            end
            @testset "colon" begin
                v[:] .= 4
                @test all(@. v[:] == 4)
            end
        end
    end
    @testset "1D normal" begin
        v = SHArray(zeros(length(mode_lm)))
        @testset "getindex" begin
            @testset "linear" begin
                for i in eachindex(v)
                    @test v[i] == 0
                end
            end
            @testset "mode_lm" begin
                for m in mode_lm
                    @test_throws ArgumentError v[m]
                end
            end
            @testset "colon" begin
                @test all(@. v[:] == 0)
                @test v[:] == zeros(size(v))
            end
        end
        @testset "setindex" begin
            @testset "linear" begin
                for i in eachindex(v)
                    v[i] = i
                    @test v[i] == i
                end
            end
            @testset "mode_lm" begin
                for (i, m) in enumerate(mode_lm)
                    @test_throws ArgumentError v[m] = i
                end
            end
            @testset "colon" begin
                v[:] .= 4
                @test all(@. v[:] == 4)
            end
        end
    end
    @testset "2D mixed axes" begin
        @testset "first SH" begin
            a = SHArray{ComplexF64}((mode_lm, 1:2))
            @testset "getindex" begin
                @test to_indices(a, (1,)) == (1,)
                @test to_indices(a, (big(1),)) == (1,)
                @testset "linear" begin
                    for i in eachindex(IndexLinear(), a)
                        @test a[i] == 0
                    end
                    for i in eachindex(IndexCartesian(), a)
                        @test a[i] == 0
                    end
                end
                @testset "mode_lm" begin
                    for ind2 in axes(a, 2), m in mode_lm
                        @test a[m, ind2] == 0
                    end
                end
                @testset "colon" begin
                    for j in axes(a, 2)
                        @test a[:,j] == zeros(size(a, 1))
                    end
                    for i in axes(a, 1)
                        @test a[i,:] == zeros(size(a, 2))
                    end
                    @test a[:,:] == zeros(size(a))
                    @test a[:] == zeros(length(a))
                end
                @test a[CartesianIndices(a)] == a
                @test a[CartesianIndices(a), 1] == a
                @test_throws ArgumentError a[1, (0,0)]
            end
            @testset "setindex" begin
                @testset "linear" begin
                    for i in eachindex(a)
                        a[i] = i
                        @test a[i] == i
                    end
                end
                @testset "mode_lm" begin
                    for ind2 in axes(a, 2),(i, m) in enumerate(mode_lm)
                        a[m, ind2] = i
                        @test a[m, ind2] == i
                    end
                end
                @testset "colon" begin
                    for j in axes(a, 2)
                        @. a[:,j] = 1
                        @test a[:,j] == ones(size(a, 1))
                        a[:,j] += a[:,j]
                        @test a[:,j] == 2 .* ones(size(a, 1))
                    end
                    for i in axes(a, 1)
                        @. a[i,:] = 1
                        @test a[i,:] == ones(size(a, 2))
                        @. a[i,:] += a[i,:]
                        @test a[i,:] == 2 .* ones(size(a, 2))
                    end
                    @. a = 1
                    @test a == ones(size(a))
                    a[:,:] .= 1
                    @test a == ones(size(a))
                    a += a[:,:]
                    @test a == 2 .* ones(size(a))
                    a[:,:] += a[:,:]
                    @test a == 4 .* ones(size(a))
                    a[:] += a[:]
                    @test a == 8 .* ones(size(a))
                end
            end
        end
        @testset "second SH" begin
            a = SHArray{ComplexF64}((1:2, mode_lm))
            @testset "getindex" begin
                @testset "linear" begin
                    for i in eachindex(a)
                        @test a[i] == 0
                    end
                end
                @testset "mode_lm" begin
                    for (i, m) in enumerate(mode_lm),ind1 in axes(a, 1)
                        @test a[ind1, m] == 0
                    end
                end
                @testset "colon" begin
                    for j in axes(a, 2)
                        @test a[:,j] == zeros(size(a, 1))
                    end
                    for i in axes(a, 1)
                        @test a[i,:] == zeros(size(a, 2))
                    end
                    @test a[:,:] == zeros(size(a))
                    @test a[:] == zeros(length(a))
                end
            end
            @testset "setindex" begin
                @testset "linear" begin
                    for i in eachindex(a)
                        a[i] = i
                        @test a[i] == i
                    end
                end
                @testset "mode_lm" begin
                    for (i, m) in enumerate(mode_lm),ind1 in axes(a, 1)
                        a[ind1, m] = i
                        @test a[ind1, m] == i
                    end
                end
                @testset "colon" begin
                    for j in axes(a, 2)
                        @. a[:,j] = 1
                        @test a[:,j] == ones(size(a, 1))
                        a[:,j] += a[:,j]
                        @test a[:,j] == 2 .* ones(size(a, 1))
                    end
                    for i in axes(a, 1)
                        @. a[i,:] = 1
                        @test a[i,:] == ones(size(a, 2))
                        @. a[i,:] += a[i,:]
                        @test a[i,:] == 2 .* ones(size(a, 2))
                    end
                    @. a = 1
                    @test a == ones(size(a))
                    a[:,:] .= 1
                    @test a == ones(size(a))
                    a += a[:,:]
                    @test a == 2 .* ones(size(a))
                    a[:,:] += a[:,:]
                    @test a == 4 .* ones(size(a))
                    a[:] += a[:]
                    @test a == 8 .* ones(size(a))
                end
            end
        end
        @testset "both SH" begin
            a = SHArray{ComplexF64}((mode_lm, mode_lm))
            @testset "getindex" begin
                @testset "linear" begin
                    for i in eachindex(a)
                        @test a[i] == 0
                    end
                end
                @testset "mode_lm" begin
                    for m2 in mode_lm, m1 in mode_lm
                        @test a[m1, m2] == 0
                    end
                end
                @testset "colon" begin
                    for j in axes(a, 2)
                        @test a[:,j] == zeros(size(a, 1))
                    end
                    for i in axes(a, 1)
                        @test a[i,:] == zeros(size(a, 2))
                    end
                    @test a[:,:] == zeros(size(a))
                    @test a[:] == zeros(length(a))
                end
            end
            @testset "setindex" begin
                @testset "linear" begin
                    for i in eachindex(a)
                        a[i] = i
                        @test a[i] == i
                    end
                end
                @testset "mode_lm" begin
                    for (i, m2) in enumerate(mode_lm),m1 in mode_lm
                        a[m1, m2] = i
                        @test a[m1, m2] == i
                    end
                end
                @testset "colon" begin
                    for j in axes(a, 2)
                        @. a[:,j] = 1
                        @test a[:,j] == ones(size(a, 1))
                        a[:,j] += a[:,j]
                        @test a[:,j] == 2 .* ones(size(a, 1))
                    end
                    for i in axes(a, 1)
                        @. a[i,:] = 1
                        @test a[i,:] == ones(size(a, 2))
                        @. a[i,:] += a[i,:]
                        @test a[i,:] == 2 .* ones(size(a, 2))
                    end
                    @. a = 1
                    @test a == ones(size(a))
                    a[:,:] .= 1
                    @test a == ones(size(a))
                    a += a[:,:]
                    @test a == 2 .* ones(size(a))
                    a[:,:] += a[:,:]
                    @test a == 4 .* ones(size(a))
                    a[:] += a[:]
                    @test a == 8 .* ones(size(a))
                end
            end
        end
    end
    @testset "nested" begin
        s = SHVector{SHVector}(undef, LM(1:1));
        v = SHVector{ComplexF64}(LM(1:2))
        @testset "linear" begin
            for (ind, m) in enumerate(first(shmodes(s)))
                s[ind] = v
                @test s[m] == v
                @test s[ind] == v
            end
        end
        @testset "colon" begin
            for ind in eachindex(s)
                s[ind] = v
            end
            @test s[:] == [v for i in axes(s, 1)]
        end
    end
    @testset "colon and ranges" begin
        sa = SHArray{ComplexF64}((mode_lm, 2))
        for i in eachindex(sa)
            sa[i] = i
        end
        @test sa[:,1] == [1, 2]
        @test sa[:,2] == [3, 4]
        @test sa[1,:] == [1, 3]
        @test sa[(0, 0),:] == [1, 3]
        @test sa[2,:] == [2, 4]
        @test sa[(1, 0),:] == [2, 4]
    end
end

@testset "views" begin
    modes = LM(0:2)
    S = SHArray(1:length(modes), modes)
    Sv = @view S[:]
    for i in eachindex(S, Sv)
        @test S[i] == Sv[i]
    end

    S = SHArray(rand(2, length(modes)), :, modes)
    Sv = @view S[1:2, 1:2]
    for I in CartesianIndices(Sv)
        @test S[I] == Sv[I]
    end

    Sv = @view S[1, LM(1,-1)]
    @test Sv[] == S[1,2]
    Sv = @view S[1, LM(1:2,-1)]
    @test Sv == S[1,2:3]
end

@testset "similar" begin
    mode_lm = LM(0:1, 0:0)
    mode_ml = ML(0:1, 0:0)
    mode_l₂l₁ = L2L1Triangle(0:2, 2,0:2)

    function testsimilar(arr)
        SA = similar(arr)
        SAT = similar(arr, eltype(arr))
        SATM = similar(arr, eltype(arr), modes(arr))
        @test SA isa typeof(arr)
        @test SAT isa typeof(arr)
        @test size(parent(SA)) == size(parent(arr))
        @test modes(SA) == modes(arr)
        @test shdims(SA) == shdims(arr)
        @test modes(SAT) == modes(arr)
        @test shdims(SAT) == shdims(arr)
        @test modes(SATM) == modes(arr)
        @test shdims(SATM) == shdims(arr)

        @test axes(SA) == axes(arr)
        @test axes(SAT) == axes(arr)
        @test axes(SATM) == axes(arr)

        SA = similar(zeros(0), Int, modes(arr)...)
        @test eltype(SA) == Int
        @test axes(SA) == axes(arr)
        @test modes(SA) == modes(arr)
        @test shdims(SA) == shdims(arr)

        SA = similar(zeros(0), modes(arr)...)
        @test eltype(SA) == eltype(SA)
        @test axes(SA) == axes(arr)
        @test modes(SA) == modes(arr)
        @test shdims(SA) == shdims(arr)
    end

    SA = SHArray{ComplexF64}((mode_lm, 1:2, mode_ml))
    testsimilar(SA)
    v = SHVector{ComplexF64}(mode_lm)
    testsimilar(v)
    M = SHMatrix{ComplexF64}(mode_lm, mode_lm)
    testsimilar(M)

    SA = SHArray(Diagonal([1,2]), (LM(0:1, 0:0),1:2))
    @test parent(similar(SA)) isa Diagonal
    @test parent(similar(SA, Float64)) isa Diagonal{Float64}
end

@testset "nested" begin
    m = LM(0:2)
    sharr = SHVector{ComplexF64}(m)
    v = [sharr, sharr]
    n = length(v)
    @test typeof(v) == Array{typeof(sharr),1}
    modesv = LM(n:n, 1:n)
    shv = SHVector(v, modesv) # get the number of modes to match
    @test typeof(shv) == SHArray{typeof(sharr),1, typeof(v),Tuple{typeof(modesv)}}
    @test_throws ArgumentError SHVector(v, LM(0:0, 0:0))
end

@testset "broadcasting" begin
    @testset "unalias" begin
        s = SHArray(zeros(3, 3), (LM(1:1), LM(1:1)))
        @test Broadcast.broadcast_unalias(s, s) === s
    end
    @testset "same ndims" begin
        @testset "SHArray SHArray" begin
            @testset "first axis" begin
                s = SHArray{ComplexF64}((LM(0:1),-1:1)); @. s = 2
                t = similar(s); @. t = 4
                s2 = SHArray{ComplexF64}((LM(2:2,-2:1),-1:1)) # dimensions and mode types match but different modes
                s3 = SHArray{ComplexF64}((ML(2:2,-2:1),-1:1)) # dimensions and modes match but different mode types
                s4 = SHArray{ComplexF64}((LM(0:2),-1:1)) # dimensions don't match
                s5 = SHArray{ComplexF64}((-1:2, LM(1:1))); # sizes match but axes don't
                s6 = SHArray{ComplexF64}((LM(0:1),0:2)); # sizes match but axes don't

                u = s + t
                @test all(u .== s[1]+t[1])
                @test u == t + s

                u = @. s + t
                @test all(u .== s[1]+t[1])
                @test u == (@. t + s)

                u = @. s + 1
                @test all(u .== s[1]+1)
                @test u == (@. 1 + s)

                u = @. 1 + s + t
                @test all(u .== s[1]+t[1]+1)

                u = @. s*t
                @test all(u .== s[1]*t[1])
                @test u == (@. t*s)

                u = @. s*t + 1
                @test all(u .== s[1]*t[1]+1)

                @test_throws ModeMismatchError s + s2
                @test_throws ModeMismatchError s + s3
                @test_throws DimensionMismatch s + s4
                @test_throws DimensionMismatch s + s5
                @test_throws DimensionMismatch s + s6
            end
            @testset "second axis" begin
                s = SHArray{ComplexF64}((-1:1, LM(0:1))); @. s = 2
                t = similar(s); @. t = 4

                u = s + t
                @test all(u .== s[1]+t[1])
                @test (t + s) == u

                u = @. s + t
                @test all(u .== s[1]+t[1])
                @test u == (@. t + s)

                u = @. s + 1
                @test all(u .== s[1]+1)
                @test u == (@. 1 + s)

                u = @. 1 + s + t
                @test all(u .== s[1]+t[1]+1)

                u = @. s*t
                @test all(u .== s[1]*t[1])
                @test u == (@. t*s)

                u = @. s*t + 1
                @test all(u .== s[1]*t[1]+1)
            end
        end
        @testset "SHArray Array" begin
            s = SHArray{ComplexF64}((-1:1, LM(0:1))); @. s = 2
            t = zeros(axes(s)); @. t = 5

            u = s + t
            @test all(u .== s[1]+t[1])
            @test u == (t + s)

            u = @. s + t
            @test all(u .== s[1]+t[1])
            @test u == (@. t + s)

            u = @. s + 1
            @test all(u .== s[1]+1)
            @test u == (@. 1 + s)

            u = @. 1 + s + t
            @test all(u .== s[1]+t[1]+1)

            u = @. s*t
            @test all(u .== s[1]*t[1])
            @test u == (@. t*s)

            u = @. s*t + 1
            @test all(u .== s[1]*t[1]+1)
        end
    end
    @testset "different ndims" begin
        @testset "SHArray SHVector" begin
            s = SHArray{ComplexF64}((LM(0:1),-1:1)); @. s = 2
            t = SHVector{ComplexF64}(first(modes(s))); @. t = 4

            u = @. s + t
            @test all(u .== s[1]+t[1])
            @test modes(u) == modes(s)
            @test u == (@. t + s)

            u = @. s*t
            @test all(u .== s[1]*t[1])
            @test modes(u) == modes(s)
            @test u == (@. t*s)

            u = @. 1 + s + t
            @test all(u .== s[1]+t[1]+1)
            @test modes(u)==modes(s)

            u = @. s*t + 1
            @test all(u .== s[1]*t[1]+1)
            @test modes(u)==modes(s)
        end
        @testset "SHArray Vector" begin
            s = SHArray{ComplexF64}((LM(0:1),-1:1)); @. s = 2
            t = zeros(axes(s, 1)); @. t = 4

            u = @. s + t
            @test all(u .== s[1]+t[1])
            @test u == (@. t + s)

            u = @. s*t
            @test all(u .== s[1]*t[1])
            @test u == (@. t*s)

            u = @. 1 + s + t
            @test all(u .== s[1]+t[1]+1)

            u = @. s*t + 1
            @test all(u .== s[1]*t[1]+1)
        end
    end
end

@testset "accessor functions" begin
    mode_lm = LM(0:1, 0:0)

    @testset "SHArray" begin
        @testset "SHArrayOneAxis" begin
            arr = zeros(2, length(mode_lm),1)
            ax = (axes(arr, 1), mode_lm, axes(arr, 3))
            sa = SHArray(arr, ax)
            @test modes(sa) == ax
            @test shmodes(sa) == (mode_lm,)
            @test shdims(sa) == (2,)
        end
        @testset "general" begin
            arr = zeros(2, length(mode_lm), length(mode_lm))
            ax = (axes(arr, 1), mode_lm, mode_lm)
            sa = SHArray(arr, ax)
            @test modes(sa) == ax
            @test shmodes(sa) == (mode_lm, mode_lm)
            @test shdims(sa) == (2, 3)
        end
    end

    @testset "SHArrayOnlyFirstAxis" begin
        @testset "SHVector" begin
            sv = SHVector{ComplexF64}(mode_lm)
            @test modes(sv) == (mode_lm,)
            @test shmodes(sv) == (mode_lm,)
            @test shdims(sv) == (1,)
        end
        @testset "SHArray" begin
            arr = zeros(length(mode_lm),2, 1)
            ax = (mode_lm, axes(arr)[2:3]...)
            sa = SHArray(arr, ax)
            @test modes(sa) == ax
            @test shmodes(sa) == (mode_lm,)
            @test shdims(sa) == (1,)
        end
    end
end

@testset "permuteflipmodes" begin
    lmax = 2;
    modes = LM(0:lmax);
    v = SHArray(1:length(modes), modes)
    v2 = SphericalHarmonicArrays.permuteflipmodes(v)
    p = modeindex.((modes,), SphericalHarmonicModes.flip(modes))
    @test parent(v2) == v[p]

    modes = ML(0:lmax);
    v = SHArray(1:length(modes), modes)
    v2 = SphericalHarmonicArrays.permuteflipmodes(v)
    p = modeindex.((modes,), SphericalHarmonicModes.flip(modes))
    @test parent(v2) == v[p]
end

@testset "show" begin
    io = IOBuffer()
    showerror(io, SphericalHarmonicArrays.ModeMismatchError(LM(1:2), LM(1:4)))
    showerror(io, SphericalHarmonicArrays.ModeMismatchError(LM(1:2), ML(1:4)))

    sa = SHArray(zeros(2))
    Base.showarg(io, sa, false)

    take!(io)
    d = Diagonal([1, 2])
    s = SHArray(d, (1:2, LM(1:1, 0:1)))
    Base.print_matrix(io, d)
    showd = String(take!(io))
    Base.print_matrix(io, s)
    shows = String(take!(io))
    @test shows == showd

    @testset for j in 1:2, i in 1:2
        srep = Base.replace_in_print_matrix(s, i, j, string(s[i, j]))
        drep = Base.replace_in_print_matrix(s, i, j, string(d[i, j]))
        srep == drep
    end
end


text:
"""
"""
function fit!(
        transformer::ImmutablePackageSingleLabelPredictionTransformer,
        varargs...;
        kwargs...
        )
    if length(varargs) == 1
        return varargs[1]
    else
        return varargs
    end
end

"""
"""
function predict(
        transformer::ImmutablePackageSingleLabelPredictionTransformer,
        single_labelpredictions::AbstractVector,
        varargs...
        )
    result = DataFrames.DataFrame()
    label_name = transformer.single_label_name
    result[label_name] = single_labelpredictions
    return result
end

"""
"""
function predict_proba(
        transformer::ImmutablePackageSingleLabelPredictionTransformer,
        varargs...;
        kwargs...
        )
    if length(varargs) == 1
        return varargs[1]
    else
        return varargs
    end
end

ext: jl


text:
# jobassign.jl
"""
job assignment related functions
"""

using StatsBase

# assigntable = nothing

# function setassigntable(atbl,itbl,gtbl)
#     global assigntable = atbl
#     global intervaltable = itbl
#     global grouptable = gtbl # basically not used
#     return
# end

function isassignable(x, lst)
    if x in lst
        return true
    else
        return false
    end
end

function findassignnum(x,atbl)
    k = []
    for i in 1:length(atbl)
        if isassignable(x,atbl[i])
            push!(k,i)
        end
    end
    k
end

# make correct assignment
function assignjob!(jtbl::DVector{Int},atbl)
    for i in 1:length(jtbl)
        for j in 1:length(jtbl[i])
            x = jtbl[i][j]
            if !isassignable(x,atbl[i])
                rw = findassignnum(x,atbl)
                k2 = 0
                for r in rw
                    k2 = findfirst(x->x==0, jtbl[r])
                    if k2 > 0
                        jtbl[i][j], jtbl[r][k2] = jtbl[r][k2], jtbl[i][j]
                        break
                    end
                end
            end
        end # j
    end # i
    jtbl
end

# penalty list for bad assignment
function listbadassign(jlst::Vector{Int}, alst::Vector{Int})
    st = vcat(alst,0)
    [(x in st) ? 0 : 1 for x in jlst]
end
function listbadassign(tbl::DVector{Int}, atbl::DVector{Int})
    [listbadassign(t,a) for (t,a) in zip(tbl,atbl)]
end

max_stars_repo_name: Yoshinobu-Ishizaki/ProcGA

max_issues_repo_name: Yoshinobu-Ishizaki/ProcGA

max_forks_repo_name: Yoshinobu-Ishizaki/ProcGA


text:
"""
**QuadratureVI**

Variational Inference solver by approximating gradients via numerical integration via Quadrature

```julia
QuadratureVI(ϵ::T=1e-5,nGaussHermite::Integer=20,optimizer::Optimizer=Adam(α=0.1))
```

**Keyword arguments**

    - `ϵ::T` : convergence criteria
    - `nGaussHermite::Int` : Number of points for the integral estimation
    - `optimizer::Optimizer` : Optimizer used for the variational updates. Should be an Optimizer object from the [GradDescent.jl](https://github.com/jacobcvt12/GradDescent.jl) package. Default is `Adam()`
"""
mutable struct QuadratureVI{T<:Real} <: NumericalVI{T}
    ϵ::T #Convergence criteria
    nIter::Integer #Number of steps performed
    optimizer_η₁::LatentArray{Optimizer} #Learning rate for stochastic updates
    optimizer_η₂::LatentArray{Optimizer} #Learning rate for stochastic updates
    nPoints::Int64 #Number of points for the quadrature
    Stochastic::Bool #Use of mini-batches
    nSamples::Int64 #Number of samples of the data
    nSamplesUsed::Int64 #Size of mini-batches
    MBIndices::Vector #Indices of the minibatch
    ρ::T #Stochastic Coefficient
    HyperParametersUpdated::Bool #To know if the inverse kernel matrix must updated
    ∇η₁::LatentArray{AbstractVector{T}}
    ∇η₂::LatentArray{AbstractArray{T}}
    ∇μE::LatentArray{AbstractVector{T}}
    ∇ΣE::LatentArray{AbstractVector{T}}
    function QuadratureVI{T}(ϵ::T,nPoints::Integer,nIter::Integer,optimizer::Optimizer,Stochastic::Bool,nSamplesUsed::Integer=1) where T
        return new{T}(ϵ,nIter,[optimizer],[optimizer],nPoints,Stochastic,1,nSamplesUsed)
    end
end

function QuadratureVI(;ϵ::T=1e-5,nGaussHermite::Integer=1000,optimizer::Optimizer=Adam(α=0.1)) where {T<:Real}
    QuadratureVI{T}(ϵ,nGaussHermite,0,optimizer,false)
end

"""
**QuadratureSVI**

Stochastic Variational Inference solver by approximating gradients via numerical integration via Quadrature

```julia
QuadratureSVI(nMinibatch::Integer;ϵ::T=1e-5,nGaussHermite::Integer=20,optimizer::Optimizer=Adam(α=0.1))
```
    -`nMinibatch::Integer` : Number of samples per mini-batches

**Keyword arguments**

    - `ϵ::T` : convergence criteria, which can be user defined
    - `nGaussHermite::Int` : Number of points for the integral estimation (for the QuadratureVI)
    - `optimizer::Optimizer` : Optimizer used for the variational updates. Should be an Optimizer object from the [GradDescent.jl](https://github.com/jacobcvt12/GradDescent.jl) package. Default is `Adam()`
"""
function QuadratureSVI(nMinibatch::Integer;ϵ::T=1e-5,nGaussHermite::Integer=20,optimizer::Optimizer=Adam(α=0.1)) where {T<:Real}
    QuadratureVI{T}(ϵ,nGaussHermite,0,optimizer,false,nMinibatch)
end

function compute_grad_expectations!(model::VGP{<:Likelihood,<:QuadratureVI})
    for k in 1:model.nLatent
        for i in 1:model.nSample
            grad_quad(model.likelihood,model.y[k][i],model.μ[k][i],model.Σ[k][i,i])
        end
    end
end

function compute_grad_expectations!(model::SVGP{<:Likelihood,<:QuadratureVI})
    μ = model.κ.*model.μ; Σ = opt_diag(model.κ.*model.Σ,model.κ)
    for k in 1:model.nLatent
        for i in 1:model.nSample
            model.inference.∇μE[k][i], model.inference.∇ΣE[k][i] = grad_quad(model.likelihood,model.y[k][i],μ[k][i],Σ[k][i,i])
        end
    end
end

function grad_quad(likelihood::Likelihood,y::Real,μ::Real,σ²::Real,nPoints::Int) where {T<:Real}
    e = expectation(Normal(μ,sqrt(σ²)),n=nPoints)
    p = e(x->pdf(likelihood,y,x),d)
    dE = e(x->gradpdf(likelihood,y,x))/p
    d²E = e(x->hessiandiagpdf(likelihood,y,x))/p
    return dE, d²E - dE^2
end

function compute_log_expectations(model::VGP{<:Likelihood,<:QuadratureVI})
end

function compute_log_expectations(model::SVGP{<:Likelihood,<:QuadratureVI})
end


text:
export Bundle
export Clifford

abstract type Bundle end

"""
    Represents a Clifford bundle, AKA the Hopf fibration of S³.

fields: total and base.
"""
struct Clifford <: Bundle
    total::S³
    base::S²
end

"""
    Clifford(w, z, A, B)

Construct a Clifford bundle with the given Complex numbers: `w`, `z`, `A` and `B`. The ratio
A:B distinguishes the fibers from one another.
"""
Clifford(w::Complex, z::Complex, A::Complex, B::Complex) = begin
    @assert(isapprox(abs(w)^2 + abs(z)^2, 1),
            "The equation |w|² + |z|² = 1 must hold for the given w and z.")
    @assert(!isapprox(abs(A), 0) || !isapprox(abs(B), 0),
            "Aw + Bz = 0, either of A or B can be zero, but not both.")
    Clifford(ComplexPlane(w, z), ComplexLine(A / B))
end

"""
    Clifford(w, z, B)

Construct a Clifford bundle with the given Complex numbers: `w`, `z` and `B`. Here, `B`
on its own identifies a point in the base space and we compute A on the fly.
"""
Clifford(w::Complex, z::Complex, B::Complex) = begin
    # Aw + Bz = 0
    # Aw = -Bz
    # A = -Bz/w
    Clifford(w, z, -B * z / w, B)
end


text:
using MPI, MPIArrays

MPI.Init()
const comm = MPI.COMM_WORLD
const rank = MPI.Comm_rank(comm)
const nb_procs = MPI.Comm_size(comm)

N = 20 # size of the matrix

# Initial tiled distribution:
A = MPIArray{Float64}(comm, (nb_procs, 1), N, N)
forlocalpart!(lp -> fill!(lp, rank), A)
sync(A)

if rank == 0
  using Plots
  heatmap(getblock(A[:,:]))
  savefig("plot.png")
end

# Clean up
free(A)
MPI.Finalize()


text:
function calc_kld(mean, stddev, sample)
    posterior = DistributionsAD.TuringMvNormal(mean, abs.(stddev))
    posterior_loglikelihood = DistributionsAD.loglikelihood(posterior, sample)

    # Calculate prior loglikelihood
    μ2 = zeros(size(mean)[1])
    σ2 = ones(size(mean)[1])
    prior = DistributionsAD.TuringMvNormal(μ2, σ2)
    prior_loglikelihood = DistributionsAD.loglikelihood(prior, sample)

    return (posterior_loglikelihood - prior_loglikelihood)
end

@testset "Trainable Distributions" begin
    @testset "TrainableMvNormal" begin
        mean =
            Float32.([
                1.3070704,
                0.21997994,
                0.045152083,
                -0.90033674,
                0.27203867,
                0.017634762,
            ])
        stddev =
            Float32.([
                0.18203342,
                0.5463153,
                1.3376112,
                -0.19215275,
                0.5753563,
                -0.19471616,
            ])
        sample =
            Float32.([
                1.0493519,
                -0.14344203,
                0.27405795,
                -1.0525246,
                0.2674973,
                -0.3156111,
            ])

        base_grads = gradient(params(mean, stddev)) do
            kl = calc_kld(mean, stddev, sample)
            return kl
        end

        l = TrainableMvNormal((2, 3))
        Flux.Zygote.@ignore copyto!(l.mean, mean)
        Flux.Zygote.@ignore copyto!(l.stddev, stddev)
        Flux.Zygote.@ignore copyto!(l.sample, sample)
        struct_grads = gradient(params(l)) do
            kl = normal_kl_divergence(l)
            return kl
        end
        for param in zip(params(l), params(mean, stddev))
            @test isapprox(struct_grads[param[1]], base_grads[param[2]])
            @test size(param[1]) == size(struct_grads[param[1]])
        end

        grads = gradient(params(l)) do
            weight = l()
            kl = scale_mixture_kl_divergence(l)
            return kl
        end
        for param in params(l)
            @test size(param) == size(grads[param])
        end
    end
end


text:
using LinearAlgebra
using CompScienceMeshes
using SauterSchwab3D

const pI   = point(6,5,3)
const pII  = point(5,2,3)
const pIII = point(7,1,0)
const pIV  = point(3,4,0)

const qII  = point(3,8,4)
const qIII = point(7,5,6)

const P = simplex(pI,pII,pIII,pIV)
const Q = simplex(pI,qII,qIII)

const sing = SauterSchwab3D.singularity_detection(P,Q)

Accuracy2 = 16
cv_ref = CommonVertex5D(sing,SauterSchwab3D._legendre(Accuracy2,0.0,1.0))

function integrand(x,y)
   ((x-pI)'*(y-pIV))*return(exp(-im*1*norm(x-y))/(4pi*norm(x-y)))
end

function INTEGRAND(u,v)
   n1 = neighborhood(P,u)
   n2 = neighborhood(Q,v)
   x = cartesian(n1)
   y = cartesian(n2)
   output = integrand(x,y)*jacobian(n1)*jacobian(n2)
   return(output)
end

print("Ref: ")
ref = sauterschwab_parameterized(INTEGRAND, cv_ref)
println(ref)
println()

res_tp =[]
res_sp = []
res_gm = []
n1 = []
n2 = []
n3 = []

#Grauss tensor product
for i in 2:1:14
   Accuracy = i
   cv = CommonVertex5D(sing,SauterSchwab3D._legendre(Accuracy,0.0,1.0))

   int_tp = sauterschwab_parameterized(INTEGRAND, cv)

   num_pts = 2*length(cv.qps)^5
   push!(n1,num_pts)

   push!(res_tp,int_tp)
end

#ShunnHam simplex product
for i in 2:1:7
   Accuracy = i
   cv_s = CommonVertex5D_S(sing,(SauterSchwab3D._shunnham3D(Accuracy),
                            SauterSchwab3D._shunnham2D(Accuracy)))

   int_sp = sauterschwab_parameterized(INTEGRAND, cv_s)

   num_pts = 2*length(cv_s.qps[1])*length(cv_s.qps[2])
   push!(n2,num_pts)
   push!(res_sp,int_sp)
end

#GrundmannMoeller simplex product
for i in 2:1:12
   Accuracy = i
   cv_gm = CommonVertex5D_S(sing,(SauterSchwab3D._grundmannMoeller3D(2*Accuracy-1),
                             SauterSchwab3D._grundmannMoeller2D(2*Accuracy-1)))

   int_gm = sauterschwab_parameterized(INTEGRAND, cv_gm)

   num_pts = 2*length(cv_gm.qps[1])*length(cv_gm.qps[2])
   push!(n3,num_pts)
   push!(res_gm,int_gm)
end

err_tp = norm.(res_tp.-ref)/norm(ref)
err_sp = norm.(res_sp.-ref)/norm(ref)
err_gm = norm.(res_gm.-ref)/norm(ref)

println(err_tp)
println(err_sp)
println(err_gm)

using Plots
plot( yaxis=:log, xaxis=:log, fontfamily="Times" , width=600, height=400)
plot!(n1,err_tp, label="Gauss Tensor-Product",markershape=:circle)
plot!(n2,err_sp, label="Simplex Tensor-Product",markershape=:rect)
#plot!(n3,err_gm, label="Simplex-Product GM",markershape=:x)
plot!(xlims=(1e1,1e5),ylims=(1e-7,1))
p = plot!(xlabel="#Quad. points/Func. evals.", ylabel="Rel. Error.", title="Common Vertex 5D",legend=:bottomleft)
savefig(p,"CommonVetrtex5DPlot.png")

using BenchmarkTools

ref = -0.3263644363080338 + 0.46042867580530505im

cv = CommonVertex5D(sing,SauterSchwab3D._legendre(6,0.0,1.0))

int_tp = sauterschwab_parameterized(INTEGRAND, cv)

num_pts = 2*length(cv.qps)^5
println(num_pts)
err_tp = norm.(int_tp.-ref)/norm(ref)
println(err_tp)
@time for i in 1:100 sauterschwab_parameterized(INTEGRAND, cv) end

cv_s = CommonVertex5D_S(sing,(SauterSchwab3D._shunnham3D(6),
SauterSchwab3D._shunnham2D(6)))

int_sp = sauterschwab_parameterized(INTEGRAND, cv_s)

num_pts = 2*length(cv_s.qps[1])*length(cv_s.qps[2])
println(num_pts)
err_sp = norm.(int_sp.-ref)/norm(ref)
println(err_sp)
@time for i in 1:100 sauterschwab_parameterized(INTEGRAND, cv_s) end

#=
cv_gm = CommonVertex5D_S((SauterSchwab3D._grundmannMoeller3D(2*6-1),
SauterSchwab3D._grundmannMoeller2D(2*6-1)))

int_gm = sauterschwab_parameterized(INTEGRAND, cv_gm)

num_pts = 2*length(cv_gm.qps[1])*length(cv_gm.qps[2])
println(num_pts)
err_gm = norm.(int_gm.-ref)/norm(ref)
println(err_gm)
=#


text:
using Parquet
using Test

function test_codec()
    println("testing reading bitpacked run (old scheme)...")
    let data = UInt8[0x05, 0x39, 0x77]
        io = PipeBuffer(data)
        decoded = Parquet.read_bitpacked_run_old(io, Int32(8), Val{UInt8(3)}())
        @test decoded == Int32[0:7;]
    end
    println("passed.")

    println("testing reading bitpacked run...")
    let data = UInt8[0x88, 0xc6, 0xfa]
        io = PipeBuffer(data)
        decoded = Parquet.read_bitpacked_run(io, 8, Val{UInt8(3)}())
        @test decoded == Int32[0:7;]
    end
    println("passed.")
end

test_codec()


text:
@testset "Solvers" begin
  grid1D = RegularGrid{Float64}(100)
  grid2D = RegularGrid{Float64}(100,100)

  @testset "Kriging" begin
    problem1D = EstimationProblem(data1D, grid1D, :value)
    problem2D = EstimationProblem(data2D, grid2D, :value)
    solver = Kriging(:value => (variogram=GaussianVariogram(range=35.),))

    solution1D = solve(problem1D, solver)
    solution2D = solve(problem2D, solver)

    # basic checks
    result = digest(solution2D)
    @test Set(keys(result)) == Set([:value])
    @test isapprox(result[:value][:mean][26,26], 1., atol=1e-8)
    @test isapprox(result[:value][:mean][51,76], 0., atol=1e-8)
    @test isapprox(result[:value][:mean][76,51], 1., atol=1e-8)

    if ismaintainer || istravis
      @testset "Plot recipe" begin
        function plot_solution(fname, solution)
          plot(solution, size=(800,400))
          png(fname)
        end
        plot_sol1D(fname) = plot_solution(fname, solution1D)
        plot_sol2D(fname) = plot_solution(fname, solution2D)
        refimg1D = joinpath(datadir,"KrigingSolution1D.png")
        refimg2D = joinpath(datadir,"KrigingSolution2D.png")
        @test test_images(VisualTest(plot_sol1D, refimg1D), popup=!istravis, tol=0.1) |> success
        @test test_images(VisualTest(plot_sol2D, refimg2D), popup=!istravis, tol=0.1) |> success
      end
    end
  end

  @testset "SeqGaussSim" begin
    @testset "Conditional" begin
      problem1D = SimulationProblem(data1D, grid1D, :value, 1)
      problem2D = SimulationProblem(data2D, grid2D, :value, 1)
      solver = SeqGaussSim(:value => (variogram=GaussianVariogram(range=35.),))

      # Random.seed!(2017)
      # solution1D = solve(problem1D, solver)
      # solution2D = solve(problem2D, solver)

      # # basic checks
      # result = digest(solution2D)
      # @test Set(keys(result)) == Set([:value])
      # @test result[:value][1][26,26] == 1.
      # @test result[:value][1][51,76] == 0.
      # @test result[:value][1][76,51] == 1.

      if ismaintainer || istravis
        @testset "Plot recipe" begin
          function plot_solution(fname, solution)
            plot(solution, size=(800,400))
            png(fname)
          end
          plot_sol1D(fname) = plot_solution(fname, solution1D)
          plot_sol2D(fname) = plot_solution(fname, solution2D)
          refimg1D = joinpath(datadir,"SeqGaussSimSolution1D.png")
          refimg2D = joinpath(datadir,"SeqGaussSimSolution2D.png")
          #@test test_images(VisualTest(plot_sol1D, refimg1D), popup=!istravis) |> success
          #@test test_images(VisualTest(plot_sol2D, refimg2D), popup=!istravis) |> success
        end
      end
    end

    @testset "Unconditional" begin
      problem1D = SimulationProblem(grid1D, :value => Float64, 1)
      problem2D = SimulationProblem(grid2D, :value => Float64, 1)
      solver = SeqGaussSim(:value => (variogram=GaussianVariogram(range=35.),))

      #solution2D = solve(problem2D, SeqGaussSim())
      # TODO: test solution correctness
    end
  end

  @testset "CookieCutter" begin
    problem2D = SimulationProblem(grid2D, Dict(:facies => Int, :property => Float64), 3)

    # TODO: test solution correctness
  end
end


text:
# The simplest mechanistic/process-based species distribution model
#
# (For herb or grass of some kind...)
#
# Load packages
using GeoData, Plots, Dates
# We will use scientific units!
using Unitful, UnitfulRecipes

# First, load a series of worldclim stacks for each month
#
# Choose the average temperature and precipitation datasets
jan_stack = GeoStack(WorldClim{Climate}, (:tavg, :prec); month=Jan) |> replace_missing
july_stack = GeoStack(WorldClim{Climate}, (:tavg, :prec); month=July) |> replace_missing

plot(jan_stack)
jan_temp, jan_prec = jan_stack
july_temp, july_prec = july_stack

# Model 1: our first model! ################################
#
# Define the simplest (e.g. worst) possible model.
# Growth is the product of temperature and precipitation.
model1(temp, prec) = temp * prec

# Run and plot both months
model1.(jan_temp, jan_prec) |> plot
model1.(july_temp, july_prec) |> plot

# Model 2: I cant stand the rain ###########################
#
# Limit the positive effect of precipitation above 50mm a month.
# Use precipitaiton as a scalar between zero and one
# (still not very mechanistic...)
model2(temp, prec) = temp * (min(prec, 50) / 50)

# Run and plot both months
plot(model2.(jan_temp, jan_prec); title="Jan")
plot(model2.(july_temp, july_prec); title="July")

# Model 3: All grown up ####################################
#
# Use a parametrised growth-rate function from the literature
#
# Define a Parent-Tardeua plant-growth response model
function parent_tardieua(ΔH_A, α, t0, t)
    ex = exp(-ΔH_A / (u"R" * t))
    return t * ex / (1 + ex^(α * (1 - (t / t0))))
end

# Define our parameters as a struct (object)
struct ParentTardieu{Δ,Al,T,S}
plot(july_temp)
    ΔH_A::Δ
    α::Al
    t0::T
    scale::S
end

# Define a response function for the parameters at a given temperature
temp_response(p::ParentTardieu, t) = p.scale * parent_tardieua(p.ΔH_A, p.α, p.t0, t)

# Define some default parameters
#
# We use units specified with the string macro u"_"
ΔH_A = 63.5u"kJ*mol^-1"
α = 3.5
t0 = 300.0u"K"
scale = 1 / parent_tardieua(ΔH_A, α, t0, t0)

params = ParentTardieu(ΔH_A, α, t0, scale)

# Plot the temperature respponse
temperatures = (0.0u"°C":1.0u"K":40.0u"°C") .|> u"K"
plot(t -> temp_response(params, t), temperatures)

# Define the model
model3(p, temp, prec) = (1 + temp_response(p, u"K"(temp * u"°C")))u"1/m" * (min(prec, 50) / 50)

# Run and plot both months
plot(model3.(Ref(params), jan_temp, jan_prec); title="Growth rate Jan")
plot(model3.(Ref(params), july_temp, july_prec); title="Growth rate July")

# What else could we do to improve this model and code?

max_stars_repo_name: EcoJulia/Julia-PhD-course-Copenhagen

max_issues_repo_name: EcoJulia/Julia-PhD-course-Copenhagen

max_forks_repo_name: EcoJulia/Julia-PhD-course-Copenhagen


text:
# Generate a certificate and key if they do not exist

rel(p::String) = joinpath(dirname(@__FILE__), p)
if !isfile(rel("keys/server.crt"))
  @static if is_unix()
    run(`mkdir -p $(rel("keys"))`)
    run(`openssl req -x509 -nodes -days 365 -newkey rsa:2048 -keyout
      $(rel("keys/server.key")) -out $(rel("keys/server.crt"))`)
  end
end

# Simple HTTPS Server

using MbedTLS, HttpServer

http = HttpHandler() do req, res
    Response("Hello Secure World!")
end

server = Server(http)
cert = MbedTLS.crt_parse_file(rel("keys/server.crt"))
key = MbedTLS.parse_keyfile(rel("keys/server.key"))

run(server, port=8002, ssl=(cert, key))


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -5.483106975000851 <= q <= 6.193029789058127)

                     add_NL_constraint(m, :(gelu(0.7547513043737739 + -0.48113681936364117*gelu(-0.10262045764995253 + -0.12743000480166566*gelu(0.5778005478095731 + 0.8852363638482865*gelu(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + -0.8793054881197468*gelu(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))) + -0.1298782680512347*gelu(0.07154777431849713 + -0.5862593464659436*gelu(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + 0.9939436563912345*gelu(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))) + -0.9820460002458025*gelu(0.1263157800600383 + 0.1285425867505281*gelu(0.5778005478095731 + 0.8852363638482865*gelu(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + -0.8793054881197468*gelu(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))) + -0.6128023058992138*gelu(0.07154777431849713 + -0.5862593464659436*gelu(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + 0.9939436563912345*gelu(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))) + gelu(-0.6861294460251406 + 0.36167652828680286*gelu(-0.10262045764995253 + -0.12743000480166566*gelu(0.5778005478095731 + 0.8852363638482865*gelu(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + -0.8793054881197468*gelu(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))) + -0.1298782680512347*gelu(0.07154777431849713 + -0.5862593464659436*gelu(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + 0.9939436563912345*gelu(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))) + -0.5440699588803759*gelu(0.1263157800600383 + 0.1285425867505281*gelu(0.5778005478095731 + 0.8852363638482865*gelu(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + -0.8793054881197468*gelu(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))) + -0.6128023058992138*gelu(0.07154777431849713 + -0.5862593464659436*gelu(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + 0.9939436563912345*gelu(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
import Base: unique, unique!, ==
using ModelingToolkit: value, operation, arguments, istree, get_observed

"""
$(TYPEDEF)

A basis over the states with parameters , independent variable  and possible exogenous controls.
It extends an `AbstractSystem` as defined in `ModelingToolkit.jl`. `f` can either be a Julia function which is able to use ModelingToolkit variables or
a vector of `eqs`.
It can be called with the typical SciML signature, meaning out of place with `f(u,p,t)`
or in place with `f(du, u, p, t)`. If control inputs are present, it is assumed that no control corresponds to
zero for all inputs. The corresponding function calls are `f(u,p,t,inputs)` and `f(du,u,p,t,inputs)` and need to
be specified fully.

If `linear_independent` is set to `true`, a linear independent basis is created from all atom function in `f`.

If `simplify_eqs` is set to `true`, `simplify` is called on `f`.

Additional keyworded arguments include `name`, which can be used to name the basis, and
`observed` for defining observeables.

# Fields
$(FIELDS)

# Example

```julia
using ModelingToolkit
using DataDrivenDiffEq

@parameters w[1:2] t
@variables u[1:2](t)

Ψ = Basis([u; sin.(w.*u)], u, parameters = p, iv = t)
```

## Note

The keyword argument `eval_expression` controls the function creation
behavior. `eval_expression=true` means that `eval` is used, so normal
world-age behavior applies (i.e. the functions cannot be called from
the function that generates them). If `eval_expression=false`,
then construction via GeneralizedGenerated.jl is utilized to allow for
same world-age evaluation. However, this can cause Julia to segfault
on sufficiently large basis functions. By default eval_expression=false.

"""
mutable struct Basis <: AbstractBasis
    """The equations of the basis"""
    eqs::Vector{Equation}
    """Dependent (state) variables"""
    states::Vector
    """Control variables"""
    controls::Vector
    """Parameters"""
    ps::Vector
    """Observed"""
    observed::Vector
    """Independent variable"""
    iv::Num
    """Internal function representation of the basis"""
    f::Function
    """Name of the basis"""
    name::Symbol
    """Internal systems"""
    systems::Vector{Basis}
end

# Helper to build function
# TODO eval -> Runtime generated
function _build_ddd_function(rhs, states, parameters, iv, eval_expression::Bool = false)

    if eval_expression
        f_oop, f_iip = eval.(build_function(rhs, value.(states), value.(parameters), [value(iv)], expression = Val{true}))
    else
        f_oop, f_iip = build_function(rhs, value.(states), value.(parameters), [value(iv)], expression = Val{false})
    end

    function f(
        u::AbstractVector{T} where T,
        p::AbstractVector{T} where T,
        t::T where T
    )
        return f_oop(u, p, t)
    end

    function f(
        du::AbstractVector{T} where T,
        u::AbstractVector{T} where T,
        p::AbstractVector{T} where T,
        t::T where T
    )
        return f_iip(du, u, p, t)
    end

    function f(
        x::AbstractMatrix{T} where T,
        p::AbstractVector{T} where T,
        t::AbstractVector{T} where T
    )
        @assert size(x, 2) == length(t) "Measurements and time points must be of equal length!"

        return hcat([f(x[:,i], p, t[i]) for i in 1:size(x, 2)]...)

    end

    function f(
        y::AbstractMatrix{T} where T,
        x::AbstractMatrix{T} where T,
        p::AbstractVector{T} where T,
        t::AbstractVector{T} where T
    )
        @assert size(x, 2) == length(t) "Measurements and time points must be of equal length!"
        @assert size(x, 2) == size(y, 2) "Measurements and preallocated output must be of equal length!"

        @inbounds for i = 1:size(x, 2)
            @views f(y[:, i], x[:, i], p, t[i])
        end

        return
    end

    # Dispatch on DiffEqBase.NullParameters
    f(u,p::DiffEqBase.NullParameters, t) = f(u,[], t)
    f(du, u,p::DiffEqBase.NullParameters, t) = f(du, u,[], t)

    # And on the controls
    f(u,p,t,input) = f(u,p,t)
    f(du,u,p,t,input) = f(du,u,p,t)

    return f
end

function _build_ddd_function(
    rhs,
    states,
    parameters,
    iv,
    controls,
    eval_expression::Bool = false,
)

    isempty(controls) &&
        return _build_ddd_function(rhs, states, parameters, iv, eval_expression)

    # Assumes zero control is zero!

    if eval_expression
        c_oop, c_iip =
            eval.(
                build_function(
                    rhs,
                    value.(states),
                    value.(parameters),
                    [value(iv)],
                    value.(controls),
                    expression = Val{true},
                ),
            )
    else
        c_oop, c_iip = build_function(
            rhs,
            value.(states),
            value.(parameters),
            [value(iv)],
            value.(controls),
            expression = Val{false},
        )
    end

    function f(
        u::AbstractVector{T} where T,
        p::AbstractVector{T} where T,
        t::T where T,
        c::AbstractVector{T} where T = zeros(eltype(u), size(controls)...),
    )
        return c_oop(u, p, t, c)
    end

    function f(
        du::AbstractVector{T} where T,
        u::AbstractVector{T} where T,
        p::AbstractVector{T} where T,
        t::T where T,
        c::AbstractVector{T} where T= zeros(eltype(u), size(controls)...),
    )
        return c_iip(du, u, p, t, c)
    end

    function f(
        x::AbstractMatrix{T} where T,
        p::AbstractVector{T} where T,
        t::AbstractVector{T} where T
    )
        @assert size(x, 2) == length(t) "Measurements and time points must be of equal length!"

        return hcat([f(x[:,i], p, t[i]) for i in 1:size(x, 2)]...)

    end

    function f(
        y::AbstractMatrix{T} where T,
        x::AbstractMatrix{T} where T,
        p::AbstractVector{T} where T,
        t::AbstractVector{T} where T
    )
        @assert size(x, 2) == length(t) "Measurements and time points must be of equal length!"
        @assert size(x, 2) == size(y, 2) "Measurements and preallocated output must be of equal length!"

        @inbounds for i = 1:size(x, 2)
            @views f(y[:, i], x[:, i], p, t[i])
        end

        return
    end

    function f(
        x::AbstractMatrix{T} where T,
        p::AbstractVector{T} where T,
        t::AbstractVector{T} where T,
        u::AbstractMatrix{T} where T
    )
        @assert size(x, 2) == length(t) "Measurements and time points must be of equal length!"
        @assert size(x, 2) == size(u, 2) "Measurements and inputs must be of equal length!"

        return hcat([f(x[:,i], p, t[i], u[:, i]) for i in 1:size(x, 2)]...)

    end

    function f(
        y::AbstractMatrix{T} where T,
        x::AbstractMatrix{T} where T,
        p::AbstractVector{T} where T,
        t::AbstractVector{T} where T,
        u::AbstractMatrix{T} where T
    )
        @assert size(x, 2) == length(t) "Measurements and time points must be of equal length!"
        @assert size(x, 2) == size(y, 2) "Measurements and preallocated output must be of equal length!"
        @assert size(x, 2) == size(u, 2) "Measurements and inputs must be of equal length!"

        @inbounds for i = 1:size(x, 2)
            @views f(y[:, i], x[:, i], p, t[i], u[:, i])
        end

        return
    end

    # Dispatch on DiffEqBase.NullParameters
    f(u,p::DiffEqBase.NullParameters, t) = f(u,[], t)
    f(du, u,p::DiffEqBase.NullParameters, t) = f(du, u,[], t)

    return f
end

## Constructors

function Basis(eqs::AbstractVector, states::AbstractVector;
    parameters::AbstractVector = [], iv = nothing,
    controls::AbstractVector = [], observed::AbstractVector = [],
    name = gensym(:Basis),
    simplify = false, linear_independent = false,
    eval_expression = false,
    kwargs...)
    iv === nothing && (iv = Variable(:t))
    iv = value(iv)
    eqs = scalarize(eqs)
    states, controls, parameters, observed = value.(scalarize(states)), value.(scalarize(controls)), value.(scalarize(parameters)), value.(scalarize(observed))

    eqs = [eq for eq in eqs if ~isequal(Num(eq),zero(Num))]

    if linear_independent
        eqs_ = create_linear_independent_eqs(eqs, simplify)
    else
        eqs_ = simplify ? ModelingToolkit.simplify.(eqs) : eqs
    end

    unique!(eqs, !simplify)

    f = _build_ddd_function(eqs, states, parameters, iv, controls, eval_expression)

    eqs = [Variable(:φ,i) ~ eq for (i,eq) ∈ enumerate(eqs_)]

    return Basis(eqs, states, controls, parameters, observed, iv, f, name, Basis[])
end

function Basis(eqs::AbstractVector{Equation}, states::AbstractVector;
    parameters::AbstractVector = [], iv = nothing,
    controls::AbstractVector = [], observed::AbstractVector = [],
    name = gensym(:Basis),
    simplify = false, linear_independent = false,
    eval_expression = false,
    kwargs...)

    lhs = [x.lhs for x in eqs]
    # We filter out 0s
    eqs_ = [Num(x.rhs) for x in eqs if ~isequal(Num(x),zero(Num))]
    if linear_independent
        eqs_ = create_linear_independent_eqs(eqs_, simplify)
    else
        eqs_ = simplify ? ModelingToolkit.simplify.(eqs_) : eqs_
    end

    isnothing(iv) && (iv = Num(Variable(:t)))
    unique!(eqs_, !simplify)

    f = _build_ddd_function(eqs_, states, parameters, iv, controls, eval_expression)

    eqs = [lhs[i] ~ eq for (i,eq) ∈ enumerate(eqs_)]

    return Basis(eqs, value.(states), value.(controls), value.(parameters), value.(observed), value(iv), f, name, Basis[])

end

function Basis(f::Function, states::AbstractVector; parameters::AbstractVector = [], controls::AbstractVector = [],
     iv = nothing, kwargs...)

    isnothing(iv) && (iv = Num(Variable(:t)))

    try
        eqs = isempty(controls) ? f(states, parameters, iv) : f(states, parameters, iv, controls)
        return Basis(eqs, states, parameters = parameters, iv = iv, controls = controls; kwargs...)
    catch e
        rethrow(e)
    end
end

function Basis(n::Int)
    @parameters t
    @variables u[1:n](t)
    return Basis(u, u, iv = t)
end

## Printing

@inline function Base.print(io::IO, x::AbstractBasis)
    state = states(x)
    ps = parameters(x)
    Base.printstyled(io, "Model $(nameof(x)) with $(length(x)) equations\n"; bold=true)
    print(io, "States :")
    if length(state) < 5
        for xi in state
            print(io, " ", xi)
        end
    else
        print(io, " ", length(state))
    end

    if !isempty(ps)
         print(io, "\nParameters :")
         if length(ps) < 5
             for p_ in ps
                 print(io, " ", p_)
             end
         else
             print(io, " ", length(ps))
         end
     end

    println(io, "\nIndependent variable: $(independent_variable(x))")
    println(io, "Equations")
    for (i,eq) ∈ enumerate(equations(x))
        if i < 5 || i == length(x)
        println(io, "$(eq.lhs) = $(eq.rhs)")
        elseif i == 5
            println(io, "...")
        else
            continue
        end
    end
end

@inline function Base.println(io::IO, x::AbstractBasis, fullview::DataType = Val{false})
    fullview == Val{false} && return print(io, x)

    state = states(x)
    ps = parameters(x)
    Base.printstyled(io, "Model $(nameof(x)) with $(length(x)) equations\n"; bold=true)
    print(io, "States :")
    if length(state) < 5
        for xi in state
            print(io, " ", xi)
        end
    else
        print(io, " ", length(state))
    end

    if !isempty(ps)
         print(io, "\nParameters :")
         if length(ps) < 5
             for p_ in ps
                 print(io, " ", p_)
             end
         else
             print(io, " ", length(ps))
         end
     end

    println(io, "\nIndependent variable: $(independent_variable(x))")
    println(io, "Equations")
    for (i,eq) ∈ enumerate(equations(x))
        println(io, "$(eq.lhs) = $(eq.rhs)")
    end
end

## Getters

"""
    $(SIGNATURES)

    Returns the internal function representing the dynamics of the `Basis`. This can be called either inplace or out-of-place
    with the typical SciML signature `f(u,p,t)` or `f(du,u,p,t)`. If control variables are defined, the function can also be called
    by `f(u,p,t,control)` or `f(du,u,p,t,control)` and assumes `control .= 0` if no control is given.
"""
function dynamics(b::AbstractBasis)
    return get_f(b)
end

"""
Returns the control variables of the `Basis`.

$(SIGNATURES)
"""
ModelingToolkit.controls(b::AbstractBasis) = b.controls

## Callable

# Fallback
(b::AbstractBasis)(args...) = b.f(args...)

# OOP
function (b::AbstractBasis)(x::AbstractVector{T} where T, p::AbstractVector{T} where T = parameters(b),
    t::T where T <: Number = independent_variable(b))
    return b.f(x,p,t)
end

function (b::AbstractBasis)(x::AbstractVector{T} where T, p::AbstractVector{T} where T,
    t::T where T <: Number , u::AbstractVector{T} where T)
    return b.f(x,p,t, u)
end

function (b::AbstractBasis)(x::AbstractMatrix{T} where T)
    t = independent_variable(b)
    return b.f(x,parameters(b),[t for i in 1:size(x,2)])
end

function (b::AbstractBasis)(x::AbstractMatrix{T} where T, p::AbstractVector{T} where T,)
    t = independent_variable(b)
    return b.f(x,p,[t for i in 1:size(x,2)])
end

function (b::AbstractBasis)(x::AbstractMatrix{T} where T, p::AbstractVector{T} where T,
    t::AbstractVector{T} where T <: Number)
    return b.f(x,p,t)
end

function (b::AbstractBasis)(x::AbstractMatrix{T} where T, p::AbstractVector{T} where T,
    t::AbstractVector{T} where T <: Number, u::AbstractMatrix{T} where T)
    return b.f(x,p,t, u)
end

# IIP
function (b::AbstractBasis)(y::AbstractMatrix{T} where T, x::AbstractMatrix{T} where T)
    t = independent_variable(b)
    return b.f(y,x,parameters(b),[t for i in 1:size(x,2)])
end

function (b::AbstractBasis)(y::AbstractMatrix{T} where T, x::AbstractMatrix{T} where T, p::AbstractVector{T} where T,)
    t = independent_variable(b)
    return b.f(y,x,p,[t for i in 1:size(x,2)])
end

function (b::AbstractBasis)(y::AbstractMatrix{T} where T, x::AbstractMatrix{T} where T, p::AbstractVector{T} where T,
    t::AbstractVector{T} where T <: Number, )
    return b.f(y,x,p,t)
end

function (b::AbstractBasis)(y::AbstractMatrix{T} where T, x::AbstractMatrix{T} where T, p::AbstractVector{T} where T,
    t::AbstractVector{T} where T <: Number, u::AbstractMatrix{T} where T)
    return b.f(y,x,p,t,u)
end

#function (b::AbstractBasis)(x::AbstractArray{T} where T, p::DiffEqBase.NullParameters, args...)
#    return b.f(x,parameters(b),args...)
#end
#
#function (b::AbstractBasis)(y::AbstractArray{T} where T, x::AbstractArray{T} where T, p::DiffEqBase.NullParameters, args...)
#    return b.f(y, x, parameters(b),args...)
#end

## Information and Iteration

Base.length(x::AbstractBasis) = length(x.eqs)
Base.size(x::AbstractBasis) = size(x.eqs)

Base.getindex(x::AbstractBasis, idx) = getindex(equations(x), idx)
Base.firstindex(x::AbstractBasis) = firstindex(equations(x))
Base.lastindex(x::AbstractBasis) = lastindex(equations(x))
Base.iterate(x::AbstractBasis) = iterate(equations(x))
Base.iterate(x::AbstractBasis, id) = iterate(equations(x), id)

## Internal update

function update!(b::AbstractBasis, eval_expression = false)

    ff = _build_ddd_function([bi.rhs for bi in equations(b)],
        states(b), parameters(b), [independent_variable(b)],
        controls(b), eval_expression)
    Core.setfield!(b, :f, ff)

    return
end

function Base.setindex!(x::AbstractBasis, idx, val, eval_expression = false)
    setindex!(equations(x), idx, val)
    update!(x, eval_expression)
    return
end
## Utilities

function unique(b::AbstractArray{Num}, simplify_eqs::Bool)
    b = simplify_eqs ? simplify.(b) : b
    returns = ones(Bool, size(b)...)
    N = maximum(eachindex(b))
    for i ∈ eachindex(b)
        returns[i] = !any([isequal(b[i], b[j]) for j in i+1:N])
    end
    return Num.(b[returns])
end

function Base.unique!(b::AbstractArray{Num}, simplify_eqs = false)
    bs = simplify_eqs ? simplify.(b) : b
    removes = zeros(Bool, size(bs)...)
    N = maximum(eachindex(bs))
    for i ∈ eachindex(bs)
        removes[i] = any([isequal(bs[i], bs[j]) for j in i+1:N])
    end
    deleteat!(b, removes)
end

function unique(b::AbstractArray{Equation}, simplify_eqs::Bool)
    b = simplify_eqs ? simplify.(b) : b
    returns = ones(Bool, size(b)...)
    N = maximum(eachindex(b))
    for i ∈ eachindex(b)
        returns[i] = !any([isequal(b[i].rhs, b[j].rhs) for j in i+1:N])
    end
    return b[returns]
end

function Base.unique!(b::AbstractArray{Equation}, simplify_eqs::Bool)
    bs = [bi.rhs for bi in b]
    bs = simplify_eqs ? simplify.(bs) : bs
    removes = zeros(Bool, size(bs)...)
    N = maximum(eachindex(bs))
    for i ∈ eachindex(bs)
        removes[i] = any([isequal(bs[i], bs[j]) for j in i+1:N])
    end
    deleteat!(b, removes)
end

## Interfacing && merging

function Base.unique!(b::Basis, simplify_eqs = false; eval_expression = false)
    unique!(b.eqs, simplify_eqs)
    update!(b, eval_expression)
end

function Base.unique(b::Basis; kwargs...)
    eqs = unique(equations(b))
    return Basis(eqs, states(b), parameters = parameters(b), iv = independent_variable(b), kwargs...)
end

"""
    deleteat!(basis, inds, eval_expression = false)

    Delete the entries specified by `inds` and update the `Basis` accordingly.
"""
function Base.deleteat!(b::Basis, inds; eval_expression = false)
    deleteat!(b.eqs, inds)
    update!(b, eval_expression)
    return
end

"""
    push!(basis, eq, simplify_eqs = true; eval_expression = false)

    Push the equations(s) in `eq` into the basis and update all internal fields accordingly.
    `eq` can either be a single equation or an array. If `simplify_eq` is true, the equation will be simplified.
"""
function Base.push!(b::Basis, eqs::AbstractArray, simplify_eqs = true; eval_expression = false)
    @inbounds for eq ∈ eqs
        push!(b, eq, false)
    end
    unique!(b, simplify_eqs, eval_expression = eval_expression)
    return
end

function Base.push!(b::Basis, eq, simplify_eqs = true; eval_expression = false)
    push!(equations(b), Variable(:φ, length(b.eqs)+1)~eq)
    unique!(b, simplify_eqs, eval_expression = eval_expression)
    return
end

function Base.push!(b::Basis, eq::Equation, simplify_eqs = true; eval_expression = false)
    push!(equations(b), eq)
    unique!(b, simplify_eqs, eval_expression = eval_expression)
    return
end

"""
    merge(x::Basis, y::Basis; eval_expression = false)

    Return a new `Basis`, which is defined via the union of `x` and `y` .
"""
function Base.merge(x::Basis, y::Basis; eval_expression = false)
    b =  unique(vcat([xi.rhs  for xi ∈ equations(x)], [xi.rhs  for xi ∈ equations(y)]))
    vs = unique(vcat(states(x), states(y)))
    ps = unique(vcat(parameters(x), parameters(y)))

    c = unique(vcat(controls(x), controls(y)))
    observed = unique(vcat(get_observed(x), get_observed(y)))

    return Basis(Num.(b), vs, parameters = ps, controls = c, observed = observed, eval_expression = eval_expression)
end

"""
    merge!(x::Basis, y::Basis; eval_expression = false)

    Updates `x` to include the union of both `x` and `y`.
"""
function Base.merge!(x::Basis, y::Basis; eval_expression = false)
    push!(x, map(x->x.rhs, equations(y)))
    Core.setfield!(x, :states, unique(vcat(states(x),states(y))))
    Core.setfield!(x, :ps, unique(vcat(parameters(x), parameters(y))))
    Core.setfield!(x, :controls, unique(vcat(controls(x), controls(y))))
    Core.setfield!(x, :observed, unique(vcat(get_observed(x), get_observed(y))))
    update!(x, eval_expression)
    return
end

## Additional functionalities

function (==)(x::Basis, y::Basis)
    length(x) == length(y) || return false
    n = zeros(Bool, length(x))
    yrhs = [yi.rhs for yi in equations(y)]
    xrhs = [xi.rhs for xi in equations(x)]
    @inbounds for (i, xi) in enumerate(xrhs)
        n[i] = any(isequal.([xi], yrhs))
        !n[i] && break
    end
    return all(n)
end

free_parameters(b::AbstractBasis; operations = [+]) = count_operation([xi.rhs for xi in b.eqs], operations) + length(b.eqs)

## Callable struct

#(b::Basis)(u, p::DiffEqBase.NullParameters, t) = b(u, [], t)
#(b::Basis)(du, u, p::DiffEqBase.NullParameters, t) = b(du, u, [], t)
#(b::Basis)(u::AbstractVector,  p::AbstractArray = [], t = nothing) = b.f(u, isempty(p) ? parameters(b) : p, isnothing(t) ? zero(eltype(u)) : t)
#(b::Basis)(du::AbstractVector, u::AbstractVector, p::AbstractArray = [], t = nothing) = b.f(du, u, isempty(p) ? parameters(b) : p, isnothing(t) ? zero(eltype(u)) : t)
#
#
#function (b::Basis)(x::AbstractMatrix, p::AbstractArray = [], t::AbstractArray = [])
#    isempty(t) ? nothing : @assert size(x, 2) == length(t)
#
#    if (isempty(p) || eltype(p) <: Num) && !isempty(parameters(b))
#        pi = isempty(p) ? parameters(b) : p
#        res = Array{Any}(undef,length(b), size(x)[2])
#    else
#        pi = p
#        res = zeros(eltype(x), length(b), size(x)[2])
#    end
#
#    @inbounds for i in 1:size(x)[2]
#        res[:, i] .= b.f(x[:, i], isempty(p) ? parameters(b) : p, isempty(t) ? zero(eltype(x)) : t[i])
#    end
#    return res
#end
#
#function (b::Basis)(y::AbstractMatrix, x::AbstractMatrix, p::AbstractArray = [], t::AbstractArray = [])
#    @assert size(x, 2) == size(y, 2) "Provide consistent arrays."
#    @assert size(y, 1) == length(b) "Provide consistent arrays."
#    isempty(t) ? nothing : @assert size(x, 2) == length(t)
#
#    @inbounds for i in 1:size(x, 2)
#        b.f(view(y, :, i), view(x, :, i), isempty(p) ? parameters(b) : p, isempty(t) ? zero(eltype(x)) : t[i])
#    end
#end
#
## Derivatives

"""
    jacobian(basis)

    Returns a function representing the jacobian matrix / gradient of the `Basis` with respect to the
    dependent variables as a function with the common signature `f(u,p,t)` for out of place and `f(du, u, p, t)` for in place computation.
    If control variables are defined, the function can also be called by `f(u,p,t,control)` or `f(du,u,p,t,control)` and assumes `control .= 0` if no control is given.
"""
function jacobian(x::Basis, eval_expression = false)

    j = Symbolics.jacobian([xi.rhs for xi in equations(x)], states(x))

    jac  = _build_ddd_function(expand_derivatives.(j),
        states(x), parameters(x), independent_variable(x),
        controls(x), eval_expression)

    return jac
end

## Get unary and binary functions

function is_unary(f::Function, t::Type = Number)
    f ∈ [+, -, *, /, ^] && return false
    for m in methods(f, (t, ))
        m.nargs - 1 > 1 && return false
    end
    return true
end

function is_binary(f::Function, t::Type = Number)
    f ∈ [+, -, *, /, ^] && return true
    !is_unary(f, t)
end

function ariety(f::Function, t::Type = Number)
    is_unary(f, t) && return 1
    is_binary(f, t) && return 2
    return 0
end

function sort_ops(f::Vector{Function})
    U = Function[]
    B = Function[]
    for fi in f
        is_unary(fi) ? push!(U, fi) : push!(B, fi)
    end
    return U, B
end

## Create linear independent basis

count_operation(x::Number, op::Function, nested::Bool = true) = 0
count_operation(x::Sym, op::Function, nested::Bool = true) = 0
count_operation(x::Num, op::Function, nested::Bool = true) = count_operation(value(x), op, nested)

function count_operation(x, op::Function, nested::Bool = true)
    if operation(x)== op
        if is_unary(op)
            # Handles sin, cos and stuff
            nested && return 1 + count_operation(arguments(x), op)
            return 1
        else
            # Handles +, *
            nested && length(arguments(x))-1 + count_operation(arguments(x), op)
            return length(arguments(x))-1
        end
    elseif nested
        return count_operation(arguments(x), op, nested)
    end
    return 0
end

function count_operation(x, ops::AbstractArray, nested::Bool = true)
    return sum([count_operation(x, op, nested) for op in ops])
end

function count_operation(x::AbstractArray, op::Function, nested::Bool = true)
    sum([count_operation(xi, op, nested) for xi in x])
end

function count_operation(x::AbstractArray, ops::AbstractArray, nested::Bool = true)
    counter = 0
    @inbounds for xi in x, op in ops
        counter += count_operation(xi, op, nested)
    end
    counter
end

function split_term!(x::AbstractArray, o, ops::AbstractArray = [+])
    if istree(o)
        n_ops = count_operation(o, ops, false)
        c_ops = 0
        @views begin
            if n_ops == 0
                x[begin]= o
            else
                counter_ = 1
                for oi in arguments(o)
                    c_ops = count_operation(oi, ops, false)
                    split_term!(x[counter_:counter_+c_ops], oi, ops)
                    counter_ += c_ops + 1
                end
            end
        end
    else
        x[begin] = o
    end
    return
end

split_term!(x::AbstractArray,o::Num, ops::AbstractArray = [+]) = split_term!(x, value(o), ops)

remove_constant_factor(x::Num) = remove_constant_factor(value(x))
remove_constant_factor(x::Number) = one(x)

function remove_constant_factor(x)
    # Return, if the function is nested
    istree(x) || return x
    # Count the number of operations
    n_ops = count_operation(x, [*], false)+1
    # Create a new array
    ops = Array{Any}(undef, n_ops)
    @views split_term!(ops, x, [*])
    filter!(x->!isa(x, Number), ops)
    return Num(prod(ops))
end

function remove_constant_factor(o::AbstractArray)
    oi = Array{Any}(undef, size(o))
    for i in eachindex(o)
        oi[i] = remove_constant_factor(o[i])
    end
    return Num.(oi)
end

function create_linear_independent_eqs(ops::AbstractVector, simplify_eqs::Bool = false)
    o = simplify.(ops)
    o = remove_constant_factor(o)
    n_ops = [count_operation(oi, +, false) for oi in o]
    n_x = sum(n_ops) + length(o)
    u_o = Array{Any}(undef, n_x)
    ind_lo, ind_up = 0, 0
    for i in eachindex(o)
        ind_lo = i > 1 ? sum(n_ops[1:i-1]) + i : 1
        ind_up = sum(n_ops[1:i]) + i

        @views split_term!(u_o[ind_lo:ind_up], o[i], [+])
    end
    u_o = remove_constant_factor(u_o)
    unique!(u_o)
    return simplify_eqs ? simplify.(Num.(u_o)) : Num.(u_o)
end


text:
@asset Cash  AFN  Afghani
@asset Cash  EUR  Euro
@asset Cash  ALL  Lek
@asset Cash  DZD  AlgerianDinar
@asset Cash  USD  USDollar
@asset Cash  AOA  Kwanza
@asset Cash  XCD  EastCaribbeanDollar
@asset Cash  ARS  ArgentinePeso
@asset Cash  AMD  ArmenianDram
@asset Cash  AWG  ArubanFlorin
@asset Cash  AUD  AustralianDollar
@asset Cash  AZN  AzerbaijanManat
@asset Cash  BSD  BahamianDollar
@asset Cash  BHD  BahrainiDinar
@asset Cash  BDT  Taka
@asset Cash  BBD  BarbadosDollar
@asset Cash  BYN  BelarusianRuble
@asset Cash  BZD  BelizeDollar
@asset Cash  XOF  CFAFrancBCEAO
@asset Cash  BMD  BermudianDollar
@asset Cash  INR  IndianRupee
@asset Cash  BTN  Ngultrum
@asset Cash  BOB  Boliviano
@asset Cash  BOV  Mvdol
@asset Cash  BAM  ConvertibleMark
@asset Cash  BWP  Pula
@asset Cash  NOK  NorwegianKrone
@asset Cash  BRL  BrazilianReal
@asset Cash  BND  BruneiDollar
@asset Cash  BGN  BulgarianLev
@asset Cash  BIF  BurundiFranc
@asset Cash  CVE  CaboVerdeEscudo
@asset Cash  KHR  Riel
@asset Cash  XAF  CFAFrancBEAC
@asset Cash  CAD  CanadianDollar
@asset Cash  KYD  CaymanIslandsDollar
@asset Cash  CLP  ChileanPeso
@asset Cash  CLF  UnidadDeFomento
@asset Cash  CNY  YuanRenminbi
@asset Cash  COP  ColombianPeso
@asset Cash  COU  UnidadDeValorReal
@asset Cash  KMF  ComorianFranc
@asset Cash  CDF  CongoleseFranc
@asset Cash  NZD  NewZealandDollar
@asset Cash  CRC  CostaRicanColon
@asset Cash  HRK  Kuna
@asset Cash  CUP  CubanPeso
@asset Cash  CUC  PesoConvertible
@asset Cash  ANG  NetherlandsAntilleanGuilder
@asset Cash  CZK  CzechKoruna
@asset Cash  DKK  DanishKrone
@asset Cash  DJF  DjiboutiFranc
@asset Cash  DOP  DominicanPeso
@asset Cash  EGP  EgyptianPound
@asset Cash  SVC  ElSalvadorColon
@asset Cash  ERN  Nakfa
@asset Cash  SZL  Lilangeni
@asset Cash  ETB  EthiopianBirr
@asset Cash  FKP  FalklandIslandsPound
@asset Cash  FJD  FijiDollar
@asset Cash  XPF  CFPFranc
@asset Cash  GMD  Dalasi
@asset Cash  GEL  Lari
@asset Cash  GHS  GhanaCedi
@asset Cash  GIP  GibraltarPound
@asset Cash  GTQ  Quetzal
@asset Cash  GBP  PoundSterling
@asset Cash  GNF  GuineanFranc
@asset Cash  GYD  GuyanaDollar
@asset Cash  HTG  Gourde
@asset Cash  HNL  Lempira
@asset Cash  HKD  HongKongDollar
@asset Cash  HUF  Forint
@asset Cash  ISK  IcelandKrona
@asset Cash  IDR  Rupiah
@asset Cash  XDR  SDR(SpecialDrawingRight)
@asset Cash  IRR  IranianRial
@asset Cash  IQD  IraqiDinar
@asset Cash  ILS  NewIsraeliSheqel
@asset Cash  JMD  JamaicanDollar
@asset Cash  JPY  Yen
@asset Cash  JOD  JordanianDinar
@asset Cash  KZT  Tenge
@asset Cash  KES  KenyanShilling
@asset Cash  KPW  NorthKoreanWon
@asset Cash  KRW  Won
@asset Cash  KWD  KuwaitiDinar
@asset Cash  KGS  Som
@asset Cash  LAK  LaoKip
@asset Cash  LBP  LebanesePound
@asset Cash  LSL  Loti
@asset Cash  ZAR  Rand
@asset Cash  LRD  LiberianDollar
@asset Cash  LYD  LibyanDinar
@asset Cash  CHF  SwissFranc
@asset Cash  MOP  Pataca
@asset Cash  MKD  Denar
@asset Cash  MGA  MalagasyAriary
@asset Cash  MWK  MalawiKwacha
@asset Cash  MYR  MalaysianRinggit
@asset Cash  MVR  Rufiyaa
@asset Cash  MRU  Ouguiya
@asset Cash  MUR  MauritiusRupee
@asset Cash  XUA  ADBUnitOfAccount
@asset Cash  MXN  MexicanPeso
@asset Cash  MXV  MexicanUnidadDeInversion(UDI)
@asset Cash  MDL  MoldovanLeu
@asset Cash  MNT  Tugrik
@asset Cash  MAD  MoroccanDirham
@asset Cash  MZN  MozambiqueMetical
@asset Cash  MMK  Kyat
@asset Cash  NAD  NamibiaDollar
@asset Cash  NPR  NepaleseRupee
@asset Cash  NIO  CordobaOro
@asset Cash  NGN  Naira
@asset Cash  OMR  RialOmani
@asset Cash  PKR  PakistanRupee
@asset Cash  PAB  Balboa
@asset Cash  PGK  Kina
@asset Cash  PYG  Guarani
@asset Cash  PEN  Sol
@asset Cash  PHP  PhilippinePeso
@asset Cash  PLN  Zloty
@asset Cash  QAR  QatariRial
@asset Cash  RON  RomanianLeu
@asset Cash  RUB  RussianRuble
@asset Cash  RWF  RwandaFranc
@asset Cash  SHP  SaintHelenaPound
@asset Cash  WST  Tala
@asset Cash  STN  Dobra
@asset Cash  SAR  SaudiRiyal
@asset Cash  RSD  SerbianDinar
@asset Cash  SCR  SeychellesRupee
@asset Cash  SLL  Leone
@asset Cash  SGD  SingaporeDollar
@asset Cash  XSU  Sucre
@asset Cash  SBD  SolomonIslandsDollar
@asset Cash  SOS  SomaliShilling
@asset Cash  SSP  SouthSudanesePound
@asset Cash  LKR  SriLankaRupee
@asset Cash  SDG  SudanesePound
@asset Cash  SRD  SurinamDollar
@asset Cash  SEK  SwedishKrona
@asset Cash  CHE  WIREuro
@asset Cash  CHW  WIRFranc
@asset Cash  SYP  SyrianPound
@asset Cash  TWD  NewTaiwanDollar
@asset Cash  TJS  Somoni
@asset Cash  TZS  TanzanianShilling
@asset Cash  THB  Baht
@asset Cash  TOP  PaAnga
@asset Cash  TTD  TrinidadAndTobagoDollar
@asset Cash  TND  TunisianDinar
@asset Cash  TRY  TurkishLira
@asset Cash  TMT  TurkmenistanNewManat
@asset Cash  UGX  UgandaShilling
@asset Cash  UAH  Hryvnia
@asset Cash  AED  UAEDirham
@asset Cash  USN  USDollar(NextDay)
@asset Cash  UYU  PesoUruguayo
@asset Cash  UYI  UruguayPesoEnUnidadesIndexadas(UI)
@asset Cash  UYW  UnidadPrevisional
@asset Cash  UZS  UzbekistanSum
@asset Cash  VUV  Vatu
@asset Cash  VES  BolívarSoberano
@asset Cash  VND  Dong
@asset Cash  YER  YemeniRial
@asset Cash  ZMW  ZambianKwacha
@asset Cash  ZWL  ZimbabweDollar
@asset Cash  XTS  CodesSpecificallyReservedForTestingPurposes
@asset Cash  XXX  TheCodesAssignedForTransactionsWhereNoCurrencyIsInvolved
@asset Bond  XBA  BondMarketsUnitEuropeanCompositeUnit(EURCO)
@asset Bond  XBB  BondMarketsUnitEuropeanMonetaryUnit(E.M.U.-6)
@asset Bond  XBC  BondMarketsUnitEuropeanUnitOfAccount9(E.U.A.-9)
@asset Bond  XBD  BondMarketsUnitEuropeanUnitOfAccount17(E.U.A.-17)
@asset Commodity  XAU  Gold
@asset Commodity  XPD  Palladium
@asset Commodity  XPT  Platinum
@asset Commodity  XAG  Silver


text:
"Tableau for the s-stage Gauss-Lobatto SFIRK method"
function TableauStochasticGLRK(s::Int)
    TableauSIRK(Symbol("StochasticGLRK" * string(s)), CoefficientsGLRK(s), CoefficientsGLRK(s))
end

"""
Tableau for the 2-stage stochastic symplectic DIRK method
  Tableau for the stochastic symplectic DIRK method
  Satisfies the conditions for Lagrange-d'Alembert integrators.
  Satisfies the conditions for strong convergence of order 1.0 for one Wiener process
"""
function TableauStochasticDIRK(c::Number=0.5)

    a_drift = [[c/2  0.0    ]
                [c    (1-c)/2]]

    b_drift = [c, 1-c]

    c_drift = [c/2, (1+c)/2]

    TableauSIRK(:StochasticSymplecticDIRK, 2, a_drift, b_drift, c_drift, 2, a_drift, b_drift, c_drift)
end


text:
module ExperimentalSVM

import StatsBase.predict

export svm, cddual, cddual_shrinking, pegasos, predict
export read_svm_data

typealias SparseOrFullMat Union(Matrix, SparseMatrixCSC)

type SVMFit
	w::Vector{Float64}
	pass::Int
	converged::Bool
end

function Base.show(io::IO, fit::SVMFit)
	@printf io "Fitted linear SVM\n"
	@printf io " * Non-zero weights: %d\n" countnz(fit.w)
	@printf io " * Iterations: %d\n" fit.pass
	@printf io " * Converged: %s\n" string(fit.converged)
end

function predict(fit::SVMFit, X::SparseOrFullMat, ret_class::Bool=true)
	n, l = size(X)
	preds = Array(Float64, l)
	for i in 1:l
		v = inner_prod(fit.w, X[:,i])
		preds[i] = ret_class ? sign(v) : v
	end
	return preds
end

include("calc.jl")

include("io.jl")

include("pegasos.jl")

include("cddual.jl")

include("cddual_shrinking.jl")

function svm(X::SparseOrFullMat,
	         Y::Vector;
	         k::Integer = 5,
	         lambda::Real = 0.1,
	         T::Integer = 100)
	pegasos(X, Y, k = k, lambda = lambda, T = T)
end

end # module


text:
# https://juliagraphics.github.io/Gtk.jl/latest/manual/builder/
using Gtk
# b = GtkBuilder(filename="path/to/myapp.glade")
#
# Alternatively, if we would store above XML definition in a Julia string myapp we can initalize the builder by
#
# b = GtkBuilder(buffer=myapp)
# Now we want to access a widget from the XML file in order to actually display it on the screen. To do so we call
#
# win = b["window1"]
# showall(win)

b = GtkBuilder(filename="myapp.glade")
win = b["window1"] # "window1" came from myapp.glade
showall(win)


text:
using Test, Random, Flux

@testset "DeepONet" begin
    @testset "dimensions" begin
        # Test the proper construction
        # Branch net
        @test size(DeepONet((32,64,72), (24,48,72), σ, tanh).branch_net.layers[end].weight) == (72,64)
        @test size(DeepONet((32,64,72), (24,48,72), σ, tanh).branch_net.layers[end].bias) == (72,)
        # Trunk net
        @test size(DeepONet((32,64,72), (24,48,72), σ, tanh).trunk_net.layers[end].weight) == (72,48)
        @test size(DeepONet((32,64,72), (24,48,72), σ, tanh).trunk_net.layers[end].bias) == (72,)
    end

    # Accept only Int as architecture parameters
    @test_throws MethodError DeepONet((32.5,64,72), (24,48,72), σ, tanh)
    @test_throws MethodError DeepONet((32,64,72), (24.1,48,72))
end


text:
using Zygote
import Base.Broadcast: broadcasted
import Zygote: Context

Zygote._forward(Context(), broadcasted, sin, [0.1, 0.2])


text:
using Trajectory
using Test

x0 = 0.
v0 = 0.
a0 = 0.
t0 = 0.
xf = 10.
vf = 0.
af = 0.
tf = 0.5

@test prod([0.0, 0.0, 0.0, 800.0, -2400.0, 1920.0] .≈ minimumjerk(x0, v0, a0, t0, xf, vf, af, tf))
x, v, a = map(x->x(tf), minimumjerkf(x0, v0, a0, t0, xf, vf, af, tf))
@test xf ≈ round(x)
@test vf ≈ round(v)
@test af ≈ round(a)


text:
using StochasticDiffEq, DiffEqDevTools, Test, Random
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
using DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_additive, prob_sde_additivesystem
Random.seed!(100)

println("Bunch of additive solves")

f_bm(u,p,t) = 0.0
f_analytic_bm(u0,p,t,W) = W
g_bm(u,p,t) = 1.0
prob = SDEProblem(SDEFunction(f_bm,g_bm,analytic=f_analytic_bm),g_bm,0.0,(0.0,1.0))
sol1 = solve(prob,SRA1(),dt=1/2^(3),adaptive=false)
sol2 = solve(prob,SOSRA(),dt=1/2^(3),adaptive=false)
sol3 = solve(prob,SKenCarp(),dt=1/2^(3),adaptive=false)
@test sol1.errors[:l2] ≈ 0.0 atol = 1e-14
@test sol2.errors[:l2] ≈ 0.0 atol = 1e-14
@test sol3.errors[:l2] ≈ 0.0 atol = 1e-14

prob = prob_sde_additive

# Test error in stepping and seeding simultaniously
sol  = solve(prob,SRA(StochasticDiffEq.constructSOSRA()),dt=1/2^(3),seed=1)
sol2 = solve(prob,SOSRA(),dt=1/2^(3),seed=1)
@test sol.t ≈ sol2.t

sol =solve(prob,SRA1(),dt=1/2^(3))
sol =solve(prob,SRA2(),dt=1/2^(3))
sol =solve(prob,SRA3(),dt=1/2^(3))
sol =solve(prob,SOSRA(),dt=1/2^(3))
sol =solve(prob,SOSRA2(),dt=1/2^(3))
sol =solve(prob,SKenCarp(),dt=1/2^(3))
sol =solve(prob,SKenCarp(nlsolve=StochasticDiffEq.NLNewton()),dt=1/2^(3))

prob = prob_sde_additivesystem

sol1 = solve(prob,SKenCarp())
sol2 = solve(prob,SKenCarp(ode_error_est=false))
@test length(sol1) < 50
@test length(sol2) < 100

sol =solve(prob,SRA(),dt=1/2^(3))
sol =solve(prob,SRA1(),dt=1/2^(3))
sol =solve(prob,SRA2(),dt=1/2^(3))
sol =solve(prob,SRA3(),dt=1/2^(3))
sol =solve(prob,SOSRA(),dt=1/2^(3))
sol =solve(prob,SOSRA2(),dt=1/2^(3))
sol =solve(prob,SKenCarp(),dt=1/2^(3))

## Convergence Testing
println("Convergence Test on MultiDimAdditive")
dts = (1/2) .^ (7:-1:4) #14->7 good plot

sim = test_convergence(dts,prob,SRA(),trajectories=100)
@test abs(sim.𝒪est[:l2]-2) < 0.1
sim2 = test_convergence(dts,prob,SRA1(),trajectories=100)
@test abs(sim2.𝒪est[:l∞]-2) <.1 #High tolerance since low dts for testing!
sim2 = test_convergence(dts,prob,SRA2(),trajectories=100)
@test abs(sim2.𝒪est[:l∞]-2) <.1 #High tolerance since low dts for testing!
sim2 = test_convergence(dts,prob,SRA3(),trajectories=100)
@test abs(sim2.𝒪est[:l∞]-2) <.1 #High tolerance since low dts for testing!
sim2 = test_convergence(dts,prob,SOSRA(),trajectories=100)
@test abs(sim2.𝒪est[:l∞]-2) <.1 #High tolerance since low dts for testing!
sim2 = test_convergence(dts,prob,SOSRA2(),trajectories=100)
@test abs(sim2.𝒪est[:l∞]-2) <.1 #High tolerance since low dts for testing!
dts = (1/2) .^ (14:-1:11) #14->7 good plot
Random.seed!(100)
sim2 = test_convergence(dts,prob,SKenCarp(),trajectories=100)
@test abs(sim2.𝒪est[:l∞]-2) <.1 #High tolerance since low dts for testing!
sim2 = test_convergence(dts,prob,SKenCarp(nlsolve=StochasticDiffEq.NLFunctional()),trajectories=100)
@test abs(sim2.𝒪est[:l∞]-2) <.1 #High tolerance since low dts for testing!


text:
import Base: reshape, dropdims, vec

"""
    reshape(x, dims::Tuple)

# Example
```julia
T = Float32
x = Var(rand(T,10,5))
y = reshape(x, 5, 10)
```
"""
reshape(x::Var, dims::Vararg{Int}) = Var(reshape(x.data,dims), ∇reshape!, (x,))
reshape(x::Node, dims) = Node(reshape, (x,dims))

vec(x::Var) = reshape(x, length(x))

function ∇reshape!(y::Var, x::Var)
    isnothing(x.grad) && return
    addto!(x.grad, y.grad)
end

doc"""
    dropdims(x, dims::Int...)

Remove the dimensions of `x` specified by dims.
"""
function dropdims(x::Var, dims::Vararg{Int})
    ydata = dropdims(x.data, dims=dims)
    Var(ydata, ∇dropdims!, (x,))
end
dropdims(x::Node, dims) = Node(dropdims, (x,dims))

function dropdims(x::Var)
    dims = ()
    for i = 1:ndims(x)
        if size(x,i) == 1
            dims = tuple(dims..., i)
        end
    end
    dropdims(x, dims...)
end

function ∇dropdims!(y::Var, x::Var)
    isnothing(x.grad) && return
    addto!(x.grad, y.grad)
end


text:
@reexport module FirstPoincareInvariants

end  # module FirstPoincareInvariants


text:
using Base
using Plots

include("initialconditions.jl")
include("kernels.jl")
include("output.jl")
include("simulationstep.jl")

function runsimulation()
    # Initial quantities
    (x, ρ, u, v, p) = generateQuantities()
    writecsvstep(0, x, ρ, u, v, p)
    println("Step $(0)")
    gr()
    plot(x, ρ)

    # SPH operands
    L = zeros(PARTICLES_COUNT)
    H = zeros(PARTICLES_COUNT)

    # Simulation steps
    for step in 1:TIME_STEPS
        (x, ρ, u, v, p, L, H) = computeSimulationStep(step, x, ρ, u, v, p, L, H)
        writecsvstep(step, x, ρ, u, v, p)
        println("Step $(step)")
    end

    (x, ρ, u, v, p)
end


text:
# This file has been obtained fom the Gridap library
# https://github.com/gridap/Gridap.jl (c) [Santiago
# Badia](mailto:santiago.badia@monash.edu) and [Francesc
# Verdugo](mailto:fverdugo@cimne.upc.edu)

"""
    @abstractmethod

Macro used in generic functions that must be overloaded by derived types.
"""
macro abstractmethod()
  quote
    error("This function belongs to an interface definition and cannot be used.")
  end
end

"""
    @notimplemented
    @notimplemented "Error message"

Macro used to raise an error, when something is not implemented.
"""
macro notimplemented(message="This function is not yet implemented")
  quote
    error($(esc(message)))
  end
end

"""
    @notimplementedif condition
    @notimplementedif condition "Error message"

Macro used to raise an error if the `condition` is true
"""
macro notimplementedif(condition,message="This function is not yet implemented")
  quote
    if $(esc(condition))
      @notimplemented $(esc(message))
    end
  end
end

"""
    @unreachable
    @unreachable "Error message"

Macro used to make sure that a line of code is never reached.
"""
macro unreachable(message="This line of code cannot be reached")
  quote
    error($(esc(message)))
  end
end

max_stars_repo_name: rrgalvan/Jalerkin

max_issues_repo_name: rrgalvan/Jalerkin

max_forks_repo_name: rrgalvan/Jalerkin


text:
using Test
using Statistics
using VectorizedStatistics

using VectorizationBase
#VectorizationBase.vsum(x::Float64) = x

@testset "Vreducibles" begin include("testVreducibles.jl") end
@testset "ArrayStats" begin include("testArrayStats.jl") end
@testset "Correlation and Covariance" begin include("testCovCor.jl") end

@testset "Sorting" begin include("testSorting.jl") end


text:
using LazySets

function bouncingBallNonlinear_specification()
    # initial set in the "down" mode
    X0_down = Hyperrectangle(low=[4.9, -0.2], high=[5.1, 0.0])
    X0 = [(2, X0_down)]

    # time horizon: 6 time units
    time_horizon = 6.0

    # specification
    O = Dict{Symbol, Any}(:T => time_horizon)

    return X0, O
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks