text:
using QuasiArrays, FillArrays, LinearAlgebra, StaticArrays, Random, Base64, Test
import QuasiArrays: AbstractQuasiFill

@testset "QuasiFill" begin
    @testset "fill array constructors and convert" begin
        for (Typ, funcs, func) in ((:QuasiZeros, :zeros, :zero), (:QuasiOnes, :ones, :one))
            @eval begin
                @test $Typ((Base.OneTo(5),)) isa AbstractQuasiVector{Float64}
                @test $Typ((Base.OneTo(5),Base.OneTo(5))) isa AbstractQuasiMatrix{Float64}

                for T in (Int, Float64)
                    ax = [1,3,4]
                    Z = $Typ{T}(ax)
                    qZ = QuasiArray($funcs(T,length(ax)),(ax,))
                    @test eltype(Z) == T
                    @test Z[1] ≡ $func(T)
                    @test Z[3] ≡ $func(T)
                    @test Z[4] ≡ $func(T)
                    @test_throws BoundsError Z[2]

                    @test QuasiArray(Z) == qZ
                    @test QuasiArray{T}(Z) == qZ
                    @test QuasiArray{T,1}(Z) == qZ

                    @test convert(AbstractQuasiArray,Z) ≡ Z
                    @test convert(AbstractQuasiArray{T},Z) ≡ AbstractQuasiArray{T}(Z) ≡ Z
                    @test convert(AbstractQuasiVector{T},Z) ≡ AbstractQuasiVector{T}(Z) ≡ Z

                    Z = $Typ{T}(ax, ax)
                    qZ = QuasiArray($funcs(T,length(ax),length(ax)),(ax,ax))
                    @test eltype(Z) == T
                    @test QuasiArray(Z) == qZ
                    @test QuasiArray{T}(Z) == qZ
                    @test QuasiArray{T,2}(Z) == qZ

                    @test convert(AbstractQuasiArray,Z) ≡ convert(AbstractQuasiFill,Z) ≡ Z
                    @test convert(AbstractQuasiArray{T},Z) ≡ convert(AbstractQuasiFill{T},Z) ≡ AbstractQuasiArray{T}(Z) ≡ Z
                    @test convert(AbstractQuasiMatrix{T},Z) ≡ convert(AbstractQuasiFill{T,2},Z) ≡ AbstractQuasiMatrix{T}(Z) ≡ Z

                    @test AbstractQuasiArray{Float32}(Z) ≡ $Typ{Float32}(ax,ax)
                    @test AbstractQuasiArray{Float32,2}(Z) ≡ $Typ{Float32}(ax,ax)
                end
            end
        end

        @test QuasiFill(1) ≡ QuasiFill{Int}(1) ≡ QuasiFill{Int,0}(1) ≡ QuasiFill{Int,0,Tuple{}}(1,())
        ax = [1,3,4]
        @test QuasiFill(1.0,ax) isa AbstractQuasiVector{Float64}
        @test QuasiFill(1.0,ax,ax) isa AbstractQuasiMatrix{Float64}
        @test QuasiFill(1,ax) ≡ QuasiFill(1,(Inclusion(ax),))
        @test QuasiFill(1,ax,ax) ≡ QuasiFill(1,(Inclusion(ax),Inclusion(ax)))
        @test eltype(QuasiFill(1.0,ax,ax)) == Float64

        for T in (Int, Float64)
            F = QuasiFill{T}(one(T), ax)
            qF = QuasiArray(fill(one(T),length(ax)), (ax,))

            @test eltype(F) == T
            @test QuasiArray(F) == qF
            @test QuasiArray{T}(F) == qF
            @test QuasiArray{T,1}(F) == qF

            F = QuasiFill{T}(one(T), ax, ax)
            qF = QuasiArray(fill(one(T),length(ax),length(ax)), (ax,ax))
            @test eltype(F) == T
            @test QuasiArray(F) == qF
            @test QuasiArray{T}(F) == qF
            @test QuasiArray{T,2}(F) == qF

            @test convert(AbstractQuasiArray,F) ≡ F
            @test convert(AbstractQuasiArray{T},F) ≡ AbstractQuasiArray{T}(F) ≡ F
            @test convert(AbstractQuasiMatrix{T},F) ≡ AbstractQuasiMatrix{T}(F) ≡ F

            @test convert(AbstractQuasiArray{Float32},F) ≡ AbstractQuasiArray{Float32}(F) ≡
                    QuasiFill{Float32}(one(Float32),ax,ax)
            @test convert(AbstractQuasiMatrix{Float32},F) ≡ AbstractQuasiMatrix{Float32}(F) ≡
                    QuasiFill{Float32}(one(Float32),ax,ax)

            @test QuasiFill{T}(F) ≡ QuasiFill{T,2}(F) ≡ typeof(F)(F) ≡ F
        end

        @test QuasiEye(ax) isa QuasiDiagonal{Float64}
        @test QuasiEye(ax) == QuasiEye{Float64}(ax)
        @test eltype(QuasiEye(ax)) == Float64

        @test QuasiEye((Inclusion(ax),)) ≡ QuasiEye(ax)

        for T in (Int, Float64)
            E = QuasiEye{T}(ax)
            M = QuasiMatrix{T}(I, ax, ax)

            @test eltype(E) == T
            @test QuasiArray(E) == M
            @test QuasiArray{T}(E) == M
            @test QuasiArray{T,2}(E) == M

            @test convert(AbstractQuasiArray,E) === E
            @test convert(AbstractQuasiArray{T},E) === E
            @test convert(AbstractQuasiMatrix{T},E) === E

            @test AbstractQuasiArray{Float32}(E) == QuasiEye{Float32}(ax)
            @test QuasiEye{T}((Inclusion(ax),)) ≡ QuasiEye{T}(ax)
        end

        @testset "Bool should change type" begin
            ax = [1,3,4]
            x = QuasiFill(true,ax)
            y = x + x
            @test y isa QuasiFill{Int,1}
            @test y[1] == 2

            x = QuasiOnes{Bool}(ax)
            y = x + x
            @test y isa QuasiFill{Int,1}
            @test y[1] == 2
            @test x + QuasiZeros{Bool}(ax) ≡ x
            @test x - QuasiZeros{Bool}(ax) ≡ x
            @test QuasiZeros{Bool}(ax) + x ≡ x
            @test -x ≡ QuasiFill(-1,ax)
        end

        @testset "copy should return Fill" begin
            ax = [1,3,4]
            x = QuasiFill(1.0,ax)
            @test copy(x) ≡ x
            x = QuasiZeros(10)
            @test copy(x) ≡ x
            x = QuasiFill([1.,2.],ax)
            @test copy(x) == x
            @test copy(x) === x   # because isbits(x)
            @test copy(x) isa QuasiFill
        end

        @testset "in" begin
            ax = [1,3,4]
            for T in [QuasiZeros, QuasiOnes, QuasiFill]
                A = T(4, ax)
                @test FillArrays.getindex_value(A) in A
                @test !(FillArrays.getindex_value(A) + 1 in A)
            end
        end
    end

    @testset "indexing" begin
        ax = [1,3,4]
        A = QuasiFill(3.0,ax)
        @test A[[1,3]] ≡ Fill(3.0,2)
        @test A[Inclusion([1,3])] isa QuasiFill
        @test_throws BoundsError A[1:3]
        @test_throws BoundsError A[1:26]
        A = QuasiFill(3.0, ax, ax)
        @test A[[1,3], [3,4]] ≡ Fill(3.0,2,2)
        @test A[[1,3], Inclusion([3,4])] == QuasiFill(3.0,(Base.OneTo(2),Inclusion([3,4])))

        A = QuasiOnes{Int}(ax,ax)
        @test A[[1,3],1] ≡ Ones{Int}(2)
        @test A[[1,3],[3,4]] ≡ Ones{Int}(2,2)
        @test_throws BoundsError A[1:26]

        A = QuasiZeros{Int}(ax,ax)
        @test A[[1,3],1] ≡ Zeros{Int}(2)
        @test A[[1,3],[3,4]] ≡ Zeros{Int}(2,2)
        @test_throws BoundsError A[1:26]

        @testset "colon" begin
            @test QuasiOnes(ax)[:] == QuasiOnes(ax)
            @test QuasiZeros(ax)[:] == QuasiZeros(ax)
            @test QuasiFill(3.0,ax)[:] == QuasiFill(3.0,ax)

            @test QuasiOnes(ax,ax)[:,:] ≡ QuasiOnes(ax,ax)
        end

        @testset "mixed integer / vector /colon" begin
            ax = [1,3,4]
            a = QuasiFill(2.0,ax)
            z = QuasiZeros(ax)
            @test a[Inclusion(ax)] ≡ a[:] ≡ a
            @test z[Inclusion(ax)] ≡ z[:] ≡ z

            bx = 2:5
            A = QuasiFill(2.0,ax,bx)
            Z = QuasiZeros(ax,bx)
            @test A[:,2] ≡ A[Inclusion(ax),2] ≡ QuasiFill(2.0,ax)
            @test A[ax,2] ≡ Fill(2.0,3)
            @test A[1,:] ≡ A[1,Base.IdentityUnitRange(bx)] ≡ Fill(2.0,(Base.IdentityUnitRange(bx),))

            @test A[:,:] == A[Inclusion(ax),Base.IdentityUnitRange(bx)] == A[Inclusion(ax),:] == A[:,Base.IdentityUnitRange(bx)] == A
            @test Z[:,2] ≡ Z[Inclusion(ax),2] ≡ QuasiZeros(ax)
            @test Z[1,:] ≡ Z[1,Base.IdentityUnitRange(bx)] ≡ Zeros((Base.IdentityUnitRange(bx),))
            @test Z[:,:] == Z[Inclusion(ax),Base.IdentityUnitRange(bx)] == Z[Inclusion(ax),:] == Z[:,Base.IdentityUnitRange(bx)] == Z

            cx = 1:2
            A = QuasiFill(2.0,ax,bx,cx)
            Z = QuasiZeros(ax,bx,cx)
            @test A[:,2,1] ≡ A[Inclusion(ax),2,1] ≡ QuasiFill(2.0,ax)
            @test A[1,:,1] ≡ A[1,Base.IdentityUnitRange(bx),1] ≡ Fill(2.0,(Base.IdentityUnitRange(bx),))
            @test A[:,:,:] ≡ A[Inclusion(ax),Base.IdentityUnitRange(bx),Base.IdentityUnitRange(cx)] ≡ A[Inclusion(ax),:,Base.IdentityUnitRange(cx)] ≡ A[:,Base.IdentityUnitRange(bx),Base.IdentityUnitRange(cx)] ≡ A
        end
    end

    # Check that all pair-wise combinations of + / - elements of As and Bs yield the correct
    # type, and produce numerically correct results.
    function test_addition_and_subtraction(As, Bs, Tout::Type)
        for A in As, B in Bs
            @test A + B isa Tout{promote_type(eltype(A), eltype(B))}
            @test QuasiArray(A + B) == QuasiArray(A) + QuasiArray(B)

            @test A - B isa Tout{promote_type(eltype(A), eltype(B))}
            @test QuasiArray(A - B) == QuasiArray(A) - QuasiArray(B)

            @test B + A isa Tout{promote_type(eltype(B), eltype(A))}
            @test QuasiArray(B + A) == QuasiArray(B) + QuasiArray(A)

            @test B - A isa Tout{promote_type(eltype(B), eltype(A))}
            @test QuasiArray(B - A) == QuasiArray(B) - QuasiArray(A)
        end
    end

    # Check that all permutations of + / - throw a `DimensionMismatch` exception.
    function test_addition_and_subtraction_dim_mismatch(a, b)
        @test_throws DimensionMismatch a + b
        @test_throws DimensionMismatch a - b
        @test_throws DimensionMismatch b + a
        @test_throws DimensionMismatch b - a
    end

    @testset "FillArray addition and subtraction" begin
        ax = [1,3,4]
        bx = 2:5
        test_addition_and_subtraction_dim_mismatch(QuasiZeros(ax), QuasiZeros(bx))
        test_addition_and_subtraction_dim_mismatch(QuasiZeros(ax), QuasiZeros{Int}(bx))
        test_addition_and_subtraction_dim_mismatch(QuasiZeros(ax), QuasiZeros(bx,bx))
        test_addition_and_subtraction_dim_mismatch(QuasiZeros(ax), QuasiZeros{Int}(bx,ax))

        # Construct FillArray for repeated use.
        rng = MersenneTwister(123456)
        A_fill, B_fill = QuasiFill(randn(rng, Float64), ax), QuasiFill(4, ax)

        # Unary +/- constructs a new FillArray.
        @test +A_fill === A_fill
        @test -A_fill === QuasiFill(-A_fill.value, ax)

        # FillArray +/- FillArray should construct a new FillArray.
        test_addition_and_subtraction([A_fill, B_fill], [A_fill, B_fill], QuasiFill)
        test_addition_and_subtraction_dim_mismatch(A_fill, QuasiFill(randn(rng), 5, 2))

        # FillArray + Array (etc) should construct a new Array using `getindex`.
        A_dense, B_dense = QuasiArray(randn(rng, 3),ax), QuasiArray([5, 4, 3],ax)
        test_addition_and_subtraction([A_fill, B_fill], [A_dense, B_dense], QuasiArray)
        test_addition_and_subtraction_dim_mismatch(A_fill, QuasiArray(randn(rng, 3, 4),ax,bx))
    end

    @testset "Other matrix types" begin
        ax = [1,3,4]
        bx = 2:5
        z = QuasiZeros(ax)
        @test QuasiDiagonal(z) == QuasiDiagonal(QuasiArray(z))

        @test QuasiDiagonal(QuasiZeros(ax,bx)) == QuasiDiagonal(z)
        @test convert(QuasiDiagonal, QuasiZeros(ax,ax)) == QuasiDiagonal(z)
        @test_throws BoundsError convert(QuasiDiagonal, QuasiZeros(ax,bx))

        @test convert(QuasiDiagonal{Int}, QuasiZeros(ax,ax)) == QuasiDiagonal(z)
        @test_throws BoundsError convert(QuasiDiagonal{Int}, QuasiZeros(bx,ax))

        @test QuasiDiagonal(QuasiEye(ax)) == QuasiEye(ax)
        @test convert(QuasiDiagonal, QuasiEye(ax)) ==  QuasiEye(ax)
        @test convert(QuasiDiagonal{Int}, QuasiEye(ax)) == QuasiEye(ax)
    end

    @testset "==" begin
        ax,bx = [1,3,4],2:5
        @test QuasiZeros(ax,bx) == QuasiFill(0,ax,bx)
        @test QuasiZeros(ax,bx) ≠ QuasiZeros(3)
        @test QuasiOnes(ax,bx) == QuasiFill(1,ax,bx)
    end

    @testset "Rank" begin
        ax,bx = [1,3,4],2:5
        @test rank(QuasiZeros(ax,bx)) == 0
        @test rank(QuasiOnes(ax,bx)) == 1
        @test rank(QuasiFill(2,ax,bx)) == 1
        @test rank(QuasiFill(0,ax,bx)) == 0
    end

    @testset "Identities" begin
        ax,bx,cx = [1,3,4],2:5,1:5
        A = QuasiArray(randn(4,5),bx,cx)
        B = QuasiArray(randn(3,4),ax,bx)
        @test 1.0 .* QuasiZeros(ax,bx) ≡ QuasiZeros(ax,bx) .* 1.0 ≡ QuasiZeros(ax,bx)
        @test QuasiZeros(ax,bx) * QuasiZeros(bx,cx) == QuasiZeros(ax,bx) * A == B * QuasiZeros(bx,cx) == QuasiZeros(ax,cx)
        @test QuasiZeros(ax,bx) * A isa QuasiZeros
        @test B * QuasiZeros(bx,cx) isa QuasiZeros
        @test_throws DimensionMismatch B * QuasiZeros(ax, ax)

        # Check multiplication by Adjoint vectors works as expected.
        @test QuasiArray(randn(3, 4),ax,bx)' * QuasiZeros(ax) === QuasiZeros(bx)
        @test QuasiArray(randn(4),bx)' * QuasiZeros(bx) == zero(Float64)
        @test QuasiArray([1, 2, 3],ax)' * QuasiZeros{Int}(ax) == zero(Int)
        @test_broken QuasiArray([SVector(1,2)', SVector(2,3)', SVector(3,4)'],ax)' * QuasiZeros{Int}(ax) == SVector(0,0)

        @test +(QuasiZeros{Float64}(ax, bx)) === QuasiZeros{Float64}(ax, bx)
        @test -(QuasiZeros{Float32}(bx, cx)) === QuasiZeros{Float32}(bx, cx)

        # `Zeros` are closed under addition and subtraction (both unary and binary).
        z1, z2 = QuasiZeros{Float64}(ax), QuasiZeros{Int}(ax)
        @test +(z1) === z1
        @test -(z1) === z1

        test_addition_and_subtraction([z1, z2], [z1, z2], QuasiZeros)
        test_addition_and_subtraction_dim_mismatch(z1, QuasiZeros{Float64}(ax, bx))

        # `Zeros` +/- `Fill`s should yield `Fills`.
        fill1, fill2 = QuasiFill(5.0, ax), QuasiFill(5, ax)
        test_addition_and_subtraction([z1, z2], [fill1, fill2], QuasiFill)
        test_addition_and_subtraction_dim_mismatch(z1, QuasiFill(5, bx))

        X = QuasiArray(randn(3, 4), ax, bx)
        for op in [+, -]

            # Addition / subtraction with same eltypes.
            @test op(QuasiZeros(ax, bx), QuasiZeros(ax, bx)) === QuasiZeros(ax, bx)
            @test_throws DimensionMismatch op(X, QuasiZeros(ax, cx))
            @test eltype(op(QuasiZeros(ax, bx), X)) == Float64

            # Different eltypes, the other way around.
            @test op(X, QuasiZeros{Float32}(ax,bx)) isa QuasiMatrix{Float64}
            @test !(op(X, QuasiZeros{Float32}(ax,bx)) === X)
            @test op(X, QuasiZeros{Float32}(ax,bx)) == X
            @test !(op(X, QuasiZeros{ComplexF64}(ax,bx)) === X)
            @test op(X, QuasiZeros{ComplexF64}(ax,bx)) == X

            # Addition / subtraction of Zeros.
            @test eltype(op(QuasiZeros{Float64}(ax,bx), QuasiZeros{Int}(ax,bx))) == Float64
            @test eltype(op(QuasiZeros{Int}(bx,ax), QuasiZeros{Float32}(bx,ax))) == Float32
            @test op(QuasiZeros{Float64}(ax,bx), QuasiZeros{Int}(ax,bx)) isa QuasiZeros{Float64}
            @test op(QuasiZeros{Float64}(ax,bx), QuasiZeros{Int}(ax,bx)) === QuasiZeros{Float64}(ax,bx)
        end

        # Zeros +/- dense where + / - have different results.
        @test +(QuasiZeros(ax,bx), X) == X && +(X, QuasiZeros(ax,bx)) == X
        @test !(QuasiZeros(ax,bx) + X === X) && !(X + QuasiZeros(ax,bx) === X)
        @test -(QuasiZeros(ax,bx), X) == -X

        # Addition with different eltypes.
        @test +(QuasiZeros{Float32}(ax,bx), X) isa QuasiMatrix{Float64}
        @test !(+(QuasiZeros{Float32}(ax,bx), X) === X)
        @test +(QuasiZeros{Float32}(ax,bx), X) == X
        @test !(+(QuasiZeros{ComplexF64}(ax,bx), X) === X)
        @test +(QuasiZeros{ComplexF64}(ax,bx), X) == X

        # Subtraction with different eltypes.
        @test -(QuasiZeros{Float32}(ax,bx), X) isa QuasiMatrix{Float64}
        @test -(QuasiZeros{Float32}(ax,bx), X) == -X
        @test -(QuasiZeros{ComplexF64}(ax,bx), X) == -X

        # Tests for ranges.
        X = QuasiArray(randn(3),ax)
        @test !(QuasiZeros(ax) + X === X)

        # test Base.zero
        @test zero(QuasiZeros(ax)) == QuasiZeros(ax)
        @test zero(QuasiOnes(ax,ax)) == QuasiZeros(ax,ax)
        @test zero(QuasiFill(0.5, ax, ax)) == QuasiZeros(ax,ax)
    end

    @testset "maximum/minimum/svd/sort" begin
        ax = [1,3,4]
        @test maximum(QuasiFill(1, ax)) == minimum(QuasiFill(1, ax)) == 1
        @test sort(QuasiOnes(ax)) == sort!(QuasiOnes(ax))
    end

    @testset "Broadcast" begin
        ax = [1,3,4]
        x = QuasiFill(5,ax)
        @test (.+)(x) ≡ x
        @test (.-)(x) ≡ -x
        @test exp.(x) ≡ QuasiFill(exp(5),ax)
        @test x .+ 1 ≡ QuasiFill(6,ax)
        @test 1 .+ x ≡ QuasiFill(6,ax)
        @test x .+ x ≡ QuasiFill(10,ax)
        @test x .+ QuasiOnes(ax) ≡ QuasiFill(6.0,ax)
        f = (x,y) -> cos(x*y)
        @test f.(x, QuasiOnes(ax)) ≡ QuasiFill(f(5,1.0),ax)

        y = QuasiOnes(ax,ax)
        @test (.+)(y) ≡ QuasiOnes(ax,ax)
        @test (.-)(y) ≡ QuasiFill(-1.0,ax,ax)
        @test exp.(y) ≡ QuasiFill(exp(1),ax,ax)
        @test y .+ 1 ≡ QuasiFill(2.0,ax,ax)
        @test y .+ y ≡ QuasiFill(2.0,ax,ax)
        @test y .* y ≡ y ./ y ≡ y .\ y ≡ y

        @test QuasiZeros{Int}(5) .+ QuasiZeros(5) isa QuasiZeros{Float64}

        # Test for conj, real and imag with complex element types
        @test conj(QuasiZeros{ComplexF64}(ax)) isa QuasiZeros{ComplexF64}
        @test conj(QuasiZeros{ComplexF64}(ax,ax)) isa QuasiZeros{ComplexF64}
        @test conj(QuasiOnes{ComplexF64}(ax)) isa QuasiOnes{ComplexF64}
        @test conj(QuasiOnes{ComplexF64}(ax,ax)) isa QuasiOnes{ComplexF64}
        @test real(QuasiZeros{Float64}(ax)) isa QuasiZeros{Float64}
        @test real(QuasiZeros{Float64}(ax,ax)) isa QuasiZeros{Float64}
        @test real(QuasiZeros{ComplexF64}(ax)) isa QuasiZeros{Float64}
        @test real(QuasiZeros{ComplexF64}(ax,ax)) isa QuasiZeros{Float64}
        @test real(QuasiOnes{Float64}(ax)) isa QuasiOnes{Float64}
        @test real(QuasiOnes{Float64}(ax,ax)) isa QuasiOnes{Float64}
        @test real(QuasiOnes{ComplexF64}(ax)) isa QuasiOnes{Float64}
        @test real(QuasiOnes{ComplexF64}(ax,ax)) isa QuasiOnes{Float64}
        @test imag(QuasiZeros{Float64}(ax)) isa QuasiZeros{Float64}
        @test imag(QuasiZeros{Float64}(ax,ax)) isa QuasiZeros{Float64}
        @test imag(QuasiZeros{ComplexF64}(ax)) isa QuasiZeros{Float64}
        @test imag(QuasiZeros{ComplexF64}(ax,ax)) isa QuasiZeros{Float64}
        @test imag(QuasiOnes{Float64}(ax)) isa QuasiZeros{Float64}
        @test imag(QuasiOnes{Float64}(ax,ax)) isa QuasiZeros{Float64}
        @test imag(QuasiOnes{ComplexF64}(ax)) isa QuasiZeros{Float64}
        @test imag(QuasiOnes{ComplexF64}(ax,ax)) isa QuasiZeros{Float64}

        @testset "Special Zeros/Ones" begin
            @test broadcast(+,QuasiZeros(ax)) ≡ broadcast(-,QuasiZeros(ax)) ≡ QuasiZeros(ax)
            @test broadcast(+,QuasiOnes(ax)) ≡ QuasiOnes(ax)

            @test QuasiZeros(ax) .* QuasiOnes(ax) ≡ QuasiZeros(ax) .* 1 ≡ QuasiZeros(ax)
            @test QuasiZeros(ax) .* QuasiFill(5.0, ax) ≡ QuasiZeros(ax) .* 5.0 ≡ QuasiZeros(ax)
            @test QuasiOnes(ax) .* QuasiZeros(ax) ≡ 1 .* QuasiZeros(ax) ≡ QuasiZeros(ax)
            @test QuasiFill(5.0, ax) .* QuasiZeros(ax) ≡ 5.0 .* QuasiZeros(ax) ≡ QuasiZeros(ax)

            @test QuasiZeros(ax) ./ QuasiOnes(ax) ≡ QuasiZeros(ax) ./ 1 ≡ QuasiZeros(ax)
            @test QuasiZeros(ax) ./ QuasiFill(5.0, ax) ≡ QuasiZeros(ax) ./ 5.0 ≡ QuasiZeros(ax)
            @test QuasiOnes(ax) .\ QuasiZeros(ax) ≡ 1 .\ QuasiZeros(ax) ≡ QuasiZeros(ax)
            @test QuasiFill(5.0, ax) .\ QuasiZeros(ax) ≡ 5.0 .\ QuasiZeros(ax) ≡ QuasiZeros(ax)

            @test conj.(QuasiZeros(5)) ≡ QuasiZeros(5)
            @test conj.(QuasiZeros{ComplexF64}(5)) ≡ QuasiZeros{ComplexF64}(5)

            @test_throws DimensionMismatch broadcast(*, QuasiOnes(ax), 1:6)
            @test_throws DimensionMismatch broadcast(*, 1:6, QuasiOnes(ax))
            @test_throws DimensionMismatch broadcast(*, QuasiFill(1,ax), 1:6)
            @test_throws DimensionMismatch broadcast(*, 1:6, QuasiFill(1,ax))

            @testset "Number" begin
                @test broadcast(*, QuasiZeros(ax), 2) ≡ broadcast(*, 2, QuasiZeros(ax)) ≡ QuasiZeros(ax)
            end

            @testset "Nested" begin
                a = QuasiArray(randn(3), ax)
                b = QuasiArray(randn(3), ax)
                @test a .\ b .* QuasiZeros(ax) ≡ QuasiZeros(ax)
                @test broadcast(*, QuasiZeros(ax), Base.Broadcast.broadcasted(\, a, b)) ≡ QuasiZeros(ax)
            end

            @testset "array-valued" begin
                a = QuasiArray(1:3, ax)
                @test broadcast(*, QuasiFill([1,2],ax), a) == broadcast(*, a, QuasiFill([1,2],ax)) == QuasiArray((1:3) .* Fill([1,2],3), ax)
                @test broadcast(*, QuasiFill([1,2],ax), QuasiZeros(ax)) == broadcast(*, QuasiZeros(ax), QuasiFill([1,2],ax))
                @test broadcast(*, QuasiFill([1,2],ax), QuasiZeros(ax)) isa QuasiFill{Vector{Float64}}
                @test broadcast(*, QuasiArray([[1,2], [3,4,5], [1,2]], ax), QuasiZeros(ax)) == broadcast(*, QuasiZeros(ax), QuasiArray([[1,2], [3,4,5], [1,2]], ax))
            end

            @testset "NaN" begin
                @test QuasiZeros(ax) ./ QuasiZeros(ax) ≡ QuasiZeros(ax) .\ QuasiZeros(ax) ≡ QuasiFill(NaN,ax)
                @test QuasiZeros{Int}(ax,ax) ./ QuasiZeros{Int}(ax) ≡ QuasiZeros{Int}(ax) .\ QuasiZeros{Int}(ax,ax) ≡ QuasiFill(NaN,ax,ax)
            end

            @testset "Addition" begin
                a = QuasiArray(1:3,ax)
                @test QuasiZeros{Int}(ax) .+ a == a .+ QuasiZeros{Int}(ax) == a .- QuasiZeros{Int}(ax) == a
                @test QuasiZeros(ax) .+ a == a .+ QuasiZeros(ax) == a .- QuasiZeros(ax) == a
                @test QuasiZeros{Int}(ax) .+ QuasiFill(1,ax) ≡ QuasiFill(1,ax) .+ QuasiZeros{Int}(ax) ≡ QuasiFill(1,ax) .- QuasiZeros{Int}(ax) ≡ QuasiFill(1,ax)
                @test_throws DimensionMismatch QuasiZeros{Int}([1,2]) .+ a
                @test_throws DimensionMismatch a .+ QuasiZeros{Int}([1,2])
            end
        end

        @testset "support Ref" begin
            @test QuasiFill(1,ax) .- 1 ≡ QuasiFill(1,ax) .- Ref(1) ≡ QuasiFill(1,ax) .- Ref(1I)
            @test QuasiFill([1 2; 3 4],ax) .- Ref(1I) == QuasiFill([0 2; 3 3],ax)
            @test Ref(1I) .+ QuasiFill([1 2; 3 4],ax) == QuasiFill([2 2; 3 5],ax)
        end

        @testset "Special Ones" begin
            a = QuasiArray(1:3,ax)
            @test QuasiOnes{Int}(ax) .* a == a .* QuasiOnes{Int}(ax) == a
            @test QuasiOnes(ax) .* a == a .* QuasiOnes(ax) == a
            @test QuasiOnes{Int}(5) .* QuasiOnes{Int}(5) == QuasiOnes{Int}(5)
            @test_throws DimensionMismatch QuasiOnes{Int}([1,2]) .* a
            @test_throws DimensionMismatch a .* QuasiOnes{Int}([1,2])
            @test_throws DimensionMismatch QuasiOnes{Int}(ax) .* QuasiOnes{Int}([1,2])
        end

        @testset "Zeros -" begin
            @test QuasiZeros(ax) - QuasiZeros(ax) ≡ QuasiZeros(ax)
            @test QuasiOnes(ax) - QuasiZeros(ax) ≡ QuasiOnes(ax)
            @test QuasiOnes(ax) - QuasiOnes(ax) isa QuasiZeros
            @test QuasiFill(1,ax) - QuasiZeros(ax) ≡ QuasiFill(1.0,ax)

            @test QuasiZeros(ax) .- QuasiZeros(ax) ≡ QuasiZeros(ax)
            @test QuasiOnes(ax) .- QuasiZeros(ax) ≡ QuasiOnes(ax)
            @test QuasiOnes(ax) .- QuasiOnes(ax) ≡ QuasiZeros(ax)
            @test QuasiFill(1,ax) .- QuasiZeros(ax) ≡ QuasiFill(1.0,ax)
        end

        @testset "Zero .*" begin
            ax = [1,3,4]
            a = QuasiArray(randn(3),ax)
            b = QuasiArray(1:3,ax)
            @test QuasiZeros{Int}(ax) .* QuasiZeros{Int}(ax) ≡ QuasiZeros{Int}(ax)
            @test a .* QuasiZeros(ax) isa QuasiZeros
            @test QuasiZeros(ax) .* a isa QuasiZeros
            @test b .* QuasiZeros(ax) isa QuasiZeros
            @test QuasiZeros(ax) .* b isa QuasiZeros
            @test_throws DimensionMismatch (1:11) .* QuasiZeros(ax)
        end
    end

    @testset "map" begin
        ax = [1,3,4]
        x = QuasiOnes(ax)
        @test map(exp,x) === QuasiFill(exp(1.0),ax)
        @test map(isone,x) === QuasiFill(true,ax)

        x = QuasiZeros(ax)
        @test map(exp,x) === exp.(x)

        x = QuasiFill(2,ax,ax)
        @test map(exp,x) === QuasiFill(exp(2),ax,ax)
    end

    @testset "0-dimensional" begin
        A = QuasiFill{Int,0,Tuple{}}(3, ())

        @test A[] ≡ 3
        @test A ≡ QuasiFill{Int,0}(3, ()) ≡ QuasiFill(3, ()) ≡ QuasiFill(3)
        @test size(A) == ()
        @test axes(A) == ()

        A = QuasiOnes{Int,0,Tuple{}}(())
        @test A[] ≡ 1
        @test A ≡ QuasiOnes{Int,0}(()) ≡ QuasiOnes{Int}(()) ≡ QuasiOnes{Int}()

        A = QuasiZeros{Int,0,Tuple{}}(())
        @test A[] ≡ 0
        @test A ≡ QuasiZeros{Int,0}(()) ≡ QuasiZeros{Int}(()) ≡ QuasiZeros{Int}()
    end

    @testset "Eye identity ops" begin
        ax = [1,3,4]
        m = QuasiEye(ax)
        D = QuasiDiagonal(QuasiFill(2,ax))

        for op in (permutedims, inv)
            @test op(m) === m
        end
        @test permutedims(D) ≡ D
        @test inv(D) ≡ QuasiDiagonal(QuasiFill(1/2,ax))

        @test copy(m) ≡ m
        @test copy(D) ≡ D
        @test LinearAlgebra.copy_oftype(m, Int) ≡ QuasiEye{Int}(ax)
        @test LinearAlgebra.copy_oftype(D, Float64) ≡ QuasiDiagonal(QuasiFill(2.0,ax))
    end

    @testset "Adjoint/Transpose/permutedims" begin
        ax = [1,3,4]
        bx = 2:5
        @test QuasiOnes{ComplexF64}(ax,bx)' ≡ transpose(QuasiOnes{ComplexF64}(ax,bx)) ≡ QuasiOnes{ComplexF64}(bx,ax)
        @test QuasiZeros{ComplexF64}(ax,bx)' ≡ transpose(QuasiZeros{ComplexF64}(ax,bx)) ≡ QuasiZeros{ComplexF64}(bx,ax)
        @test QuasiFill(1+im, ax, bx)' ≡ QuasiFill(1-im, bx,ax)
        @test transpose(QuasiFill(1+im, ax, bx)) ≡ QuasiFill(1+im, bx,ax)
        @test QuasiOnes(ax)' isa QuasiAdjoint # Vectors still need special dot product
        @test QuasiFill([1+im 2; 3 4; 5 6], ax, bx)' == QuasiFill([1+im 2; 3 4; 5 6]',bx, ax)
        @test transpose(QuasiFill([1+im 2; 3 4; 5 6], ax, bx)) == QuasiFill(transpose([1+im 2; 3 4; 5 6]),bx, ax)

        @test permutedims(QuasiOnes(ax)) ≡ QuasiOnes(1,ax)
        @test permutedims(QuasiZeros(ax)) ≡ QuasiZeros(1,ax)
        @test permutedims(QuasiFill(2.0,ax)) ≡ QuasiFill(2.0,1,ax)
        @test permutedims(QuasiOnes(ax,bx)) ≡ QuasiOnes(bx,ax)
        @test permutedims(QuasiZeros(ax,bx)) ≡ QuasiZeros(bx,ax)
        @test permutedims(QuasiFill(2.0,ax,bx)) ≡ QuasiFill(2.0,bx,ax)
    end

    @testset "setindex!/fill!" begin
        ax = [1,3,4]
        bx = 2:5

        F = QuasiFill(1,ax)
        @test (F[1] = 1) == 1
        @test_throws BoundsError (F[11] = 1)
        @test_throws ArgumentError (F[3] = 2)

        F = QuasiFill(1,ax,bx)
        @test (F[3,3] = 1) == 1
        @test_throws BoundsError (F[51] = 1)
        @test_throws BoundsError (F[1,6] = 1)
        @test_throws ArgumentError (F[1,2] = 2)

        @test (F[:,2] .= 1) == QuasiFill(1,ax)
        @test_throws ArgumentError (F[:,2] .= 2)

        @test fill!(F,1) == F
        @test_throws ArgumentError fill!(F,2)
    end

    @testset "mult" begin
        ax = [1,3,4]
        bx = 2:5
        cx = 1:3

        @test QuasiFill(2,ax)*QuasiFill(3,Base.OneTo(1),bx) ≡ QuasiFill(6,ax,bx)
        @test QuasiFill(2,ax,bx)*QuasiFill(3,bx,cx) ≡ QuasiFill(6*length(bx),ax,cx)
        @test QuasiFill(2,ax,bx)*QuasiFill(3,bx) ≡ QuasiFill(6*length(bx),ax)
        @test_throws DimensionMismatch QuasiFill(2,bx)*QuasiFill(3,2,cx)
        @test_throws DimensionMismatch QuasiFill(2,3,bx)*QuasiFill(3,2,cx)

        @test QuasiOnes(bx)*QuasiFill(3,Base.OneTo(1),cx) ≡ QuasiFill(3.0,bx,cx)
        @test QuasiOnes(ax,bx)*QuasiFill(3,bx,cx) == QuasiFill(3length(bx),ax,cx)
        @test QuasiOnes(ax,bx)*QuasiFill(3,bx) == QuasiFill(3length(bx),ax)

        @test QuasiFill(2,ax)*QuasiOnes(Base.OneTo(1),bx) ≡ QuasiFill(2.0,ax,bx)
        @test QuasiFill(2,ax,bx)*QuasiOnes(bx,cx) == QuasiFill(2length(bx),ax,cx)
        @test QuasiFill(2,ax,bx)*QuasiOnes(bx) == QuasiFill(2length(bx),ax)

        @test QuasiOnes(ax)*QuasiOnes(Base.OneTo(1),bx) == QuasiOnes(ax,bx)
        @test QuasiOnes(ax,bx)*QuasiOnes(bx,cx) == QuasiFill(length(bx),ax,cx)
        @test QuasiOnes(ax,bx)*QuasiOnes(bx) == QuasiFill(length(bx),ax)

        @test QuasiZeros(ax)*QuasiFill(3,Base.OneTo(1),bx) ≡   QuasiZeros(ax,bx)
        @test QuasiZeros(ax,bx)*QuasiFill(3,bx,cx) ≡ QuasiZeros(ax,cx)
        @test QuasiZeros(ax,bx)*QuasiFill(3,bx) ≡    QuasiZeros(ax)

        @test QuasiFill(2,ax)*  QuasiZeros(Base.OneTo(1),bx) ≡  QuasiZeros(ax,bx)
        @test QuasiFill(2,ax,bx)*QuasiZeros(bx,cx) ≡ QuasiZeros(ax,cx)
        @test QuasiFill(2,ax,bx)*QuasiZeros(bx) ≡    QuasiZeros(ax)

        @test QuasiZeros(ax)*QuasiZeros(Base.OneTo(1),bx) ≡ QuasiZeros(ax,bx)
        @test QuasiZeros(ax,bx)*QuasiZeros(bx,cx) ≡ QuasiZeros(ax,cx)
        @test QuasiZeros(ax,bx)*QuasiZeros(bx) ≡ QuasiZeros(ax)

        a = QuasiArray(randn(4),bx)
        A = QuasiArray(randn(1,4),Base.OneTo(1),bx)

        @test QuasiFill(2,ax,Base.OneTo(1))*A ≈ QuasiMatrix(QuasiFill(2,ax,Base.OneTo(1)))*A
        @test QuasiFill(2,ax,bx)*a ≈ QuasiMatrix(QuasiFill(2,ax,bx))*a
        @test QuasiOnes(ax,Base.OneTo(1))*A ≈ QuasiMatrix(QuasiOnes(ax,Base.OneTo(1)))*A
        @test QuasiOnes(ax,bx)*a ≈ QuasiMatrix(QuasiOnes(ax,bx))*a
        @test QuasiZeros(ax,Base.OneTo(1))*A == QuasiZeros(ax,bx)
        @test QuasiZeros(ax,bx)*a == QuasiZeros(ax)

        @test A*QuasiFill(2,bx) ≈ A*QuasiVector(QuasiFill(2,bx))
        @test A*QuasiFill(2,bx,cx) ≈ A*QuasiMatrix(QuasiFill(2,bx,cx))
        @test A*QuasiOnes(bx) ≈ A*QuasiVector(QuasiOnes(bx))
        @test A*QuasiOnes(bx,cx) ≈ A*QuasiMatrix(QuasiOnes(bx,cx))
        @test A*QuasiZeros(bx)  ≡ Zeros((Base.OneTo(1),))
        @test A*QuasiZeros(bx,cx) ≡ QuasiZeros(Base.OneTo(1),cx)

        D = QuasiDiagonal(a)
        @test QuasiZeros(ax,bx)*D ≡ QuasiZeros(ax,bx)
        @test D*QuasiZeros(bx,cx) ≡ QuasiZeros(bx,cx)
        @test D*QuasiZeros(bx) ≡ QuasiZeros(bx)

        D = QuasiDiagonal(QuasiFill(2,bx))
        @test D * QuasiOnes(bx) == QuasiFill(2.0,bx)
        @test D * QuasiOnes(bx,cx) == QuasiFill(2.0,bx,cx)
        @test QuasiOnes(ax,bx) * D == QuasiFill(2.0,ax,bx)
    end

    @testset "zeros/ones/fill" begin
        ax = Inclusion([1,3,4])
        @test zero(ax) ≡ QuasiZeros{Int}((ax,))
        @test one(ax) ≡ QuasiOnes{Int}((ax,))
        @test zeros(ax) ≡ zeros(Float64,ax) ≡ QuasiZeros((ax,))
        @test zeros(Base.OneTo(3), ax) ≡ zeros(Float64,Base.OneTo(3),ax) ≡ QuasiZeros((Base.OneTo(3),ax))
        @test ones(ax) ≡ ones(Float64,ax) ≡ QuasiOnes((ax,))
        @test ones(Base.OneTo(3), ax) ≡ ones(Float64,Base.OneTo(3),ax) ≡ QuasiOnes((Base.OneTo(3),ax))
        @test fill(2,ax) ≡ QuasiFill(2,(ax,))
        @test fill(2,Base.OneTo(3), ax) ≡ QuasiFill(2,(Base.OneTo(3),ax))
    end

    @testset "Static eval" begin
        o = QuasiOnes([1,3,4])
        @test view(o,SVector(1,3)) ≡ o[SVector(1,3)] ≡ Ones((SOneTo(2),))
    end

    @testset "UniformScaling" begin
        A = QuasiArray(randn(3,3), 0:0.5:1,0:0.5:1)
        @test A+2I ≈ 2I+A ≈ QuasiArray(A.parent+2I, 0:0.5:1,0:0.5:1)
        @test A-2I ≈ -(2I-A) ≈ QuasiArray(A.parent-2I, 0:0.5:1,0:0.5:1)
        @test A*(2I) ≈ (2I)*A ≈ 2A
        @test A/(2I) ≈ (2I)\A ≈ A/2

        B = QuasiArray(randn(3,3), 0:0.5:1,1:0.5:2)
        @test_throws DimensionMismatch B+I
        @test_throws DimensionMismatch I+B
        @test B*I ≈ I*B ≈ B
    end

    @testset "show" begin
        @test stringmime("text/plain",ones(Inclusion([1,2,3]))) == "ones(Inclusion([1, 2, 3]))"
        @test stringmime("text/plain",zeros(Inclusion([1,2,3]))) == "zeros(Inclusion([1, 2, 3]))"
        @test stringmime("text/plain",fill(2,Inclusion([1,2,3]))) == "fill(2, Inclusion([1, 2, 3]))"
    end
end


text:
## Module for Introspection
module Introspection

import SymPyCall: Sym, asSymbolic
using PythonCall
import PythonCall: getpy
export args, func, funcname, class, classname, getmembers

# utilities

"""
    funcname(x)

Return name or ""
"""
function funcname(x::Union{Sym, Py})
    y = getpy(x)
    func = pygetattr(y, "func", nothing)
    func == nothing && return ""
    return string(func.__name__)
end

"""
   func(x)

Return function head from an expression

[Invariant:](http://docs.sympy.org/dev/tutorial/manipulation.html)

Every well-formed SymPy expression `ex` must either have `length(args(ex)) == 0` or
`func(ex)(args(ex)...) = ex`.
"""
func(ex::Sym) = return getpy(ex).func

"""
    args(x)

Return arguments of `x`, as a tuple. (Empty if no `:args` property.)
"""
function args(x::Union{Sym, Py})
    asSymbolic(pygetattr(getpy(x), "args", ()))
end

# return class or nothing
class(x::T) where {T <: Union{Sym, Py}} = getattr(x, "__class__", nothing)
classname(x::T) where {T <: Union{Sym, Py}} = (cls = class(x); isnothing(cls) ? "" : cls.__name__)

#function getmembers(x::T) where {T <: Union{Sym, PyObject}}
#    Dict(u=>v for (u,v) in inspect.getmembers(x))
#end

end

ext: jl

max_stars_repo_name: jverzani/SymPyCall

max_issues_repo_name: jverzani/SymPyCall

max_forks_repo_name: jverzani/SymPyCall


text:
using CycUtils: multinomial
using Iterators: chain, partition
using JLD: jldopen

type SampleFactory
    material :: String
    amount :: UInt64
    consumed :: Ref{UInt64}
    neg_prob :: Vector{Float32}
    producer :: Function
    corpus_type :: Union{Type{String}, Type{UInt32}}
    function SampleFactory(material::String, neg_prob::Vector{Float32})
        global batch_producer
        global line_batch_producer
        if endswith(material, ".idx")
            amt = open(material) do f
                seekend(f)
                position(f)
            end
            return new(material, amt, Ref{UInt64}(0), neg_prob, batch_producer,
                       UInt32)
        end

        if !endswith(material, ".bz2")
            error("only accept bzip2 compressed plain text as corpus")
        end
        dn = dirname(material)
        fn = basename(material)[1:end-4]
        mn = joinpath(dn, fn * "-meta.jld")
        if isfile(mn)
            lines = jldopen(mn) do f
                read(f, "lines")
            end
        else
            nl, nw, nb = split(readstring(pipeline(`bzcat $material`, `wc`)))
            nl = parse(Int, nl)
            nw = parse(Int, nw)
            nb = parse(Int, nb)
            lines = nl
            jldopen(mn, "w") do f
                write(f, "lines", lines)
                write(f, "words", nw)
                write(f, "bytes", nb)
            end
        end
        return new(material, lines, Ref{UInt64}(0), neg_prob,
                   line_batch_producer, String)
    end
end

progress(factory::SampleFactory) = factory.consumed[]/factory.amount

const sampleFactory = Ref{Nullable{SampleFactory}}(Nullable{SampleFactory}())
const genBufFac = Ref{UInt32}(0xa)

function line_sampler(factory, m, k, word2idx, window=6)
    wsize = 1:window
    fn = factory.material
    lines = eachline(`bzcat $fn`)
    for line in lines
        sent = [word2idx[w] for w in split(line) if haskey(word2idx, w)]
        lsent = length(sent)
        for i = 1:lsent
            w = rand(wsize)
            st = max(1, i-w)
            en = min(lsent, i+w)

            for j in chain(st:i-1, i+1:en)
                produce((sent[i], sent[j]))
            end
        end
        factory.consumed[] += 1
    end
end

function line_batch_producer(factory, m, k)
    global init_model
    global genBufFac
    const multiplier = get(genBufFac)
    draw! = multinomial(factory.neg_prob)
    const nb_samples = m * multiplier
    idxes = collect(1:nb_samples)

    bigbatch_w = zeros(UInt32, nb_samples)
    bigbatch_c = zeros(UInt32, k+1, nb_samples)

    word2idx = jldopen(init_model) do f
        read(f, "word2idx")
    end
    corpus = Task(()->line_sampler(factory, m, k, word2idx))
    for x in partition(corpus, nb_samples)
        for (i, wc) in enumerate(x)
            w, c = wc
            bigbatch_w[i] = w
            bigbatch_c[1, i] = c
        end
        draw!(@view bigbatch_c[2:end, :])
        shuffle!(idxes)
        bigbatch_w = bigbatch_w[idxes]
        bigbatch_c = bigbatch_c[:, idxes]
        for st in 1:m:nb_samples
            en = st + m - 1
            produce((bigbatch_w[st:en], bigbatch_c[:, st:en]))
        end
    end
end

function batch_producer(factory, m, k)
    global genBufFac
    const multiplier = get(genBufFac)
    draw! = multinomial(factory.neg_prob)
    const nb_samples = m * multiplier
    idxes = collect(1:(2*nb_samples))
    buf = zeros(UInt32, (2, length(idxes)))
    const nbblk = sizeof(UInt32) * 2 * nb_samples

    open(factory.material) do f
        nbtotoal = factory.amount

        for i = 1:div(nbtotoal, nbblk)
            buf_pos = read(f, UInt32, (2, nb_samples))
            factory.consumed[] += nbblk
            buf[:, 1:nb_samples] = buf_pos
            buf[1, nb_samples+1:end] = buf_pos[2, :]
            buf[2, nb_samples+1:end] = buf_pos[1, :]
            shuffle!(idxes)
            buf_o = buf
            buf = reshape(buf_o[:, idxes], (2, m, 2*multiplier))
            for kk = 1:2*multiplier
                c = zeros(UInt32, (k+1, m))
                w = buf[1, :, kk]
                c[1, :] = buf[2, :, kk]
                draw!(@view c[2:end, :])
                produce((w, c))
            end
            buf = buf_o
        end
    end
    nothing
end


text:
# this is a little Monte Carlo exercise that illustrates that
# the OLS estimator is efficient. It compares the OLS estimator
# to an estimator that is the average of the OLS estimators using
# 3 subsamples
using Plots
function main()
pyplot()
reps = 10000 # number of Monte Carlo reps.
n = 21 # sample size
sig = 2.0  # st. dev. of errors

x = [ones(n,1) randn(n,1)]  # x is fixed over repeated samples
beta = [1.0, 2.0] # true beta
PopRegLine = x*beta

betas = zeros(reps,2) # holder for results

for i = 1:reps
	# generate dep var
	e = sig*randn(n)
	y = PopRegLine + e
	# the OLS estimator
	betas[i,1] = (y\x)[2]
	# the average of split sample estimator
	y1 = y[1:7,:]
	y2 = y[8:14,:]
	y3 = y[15:21,:]
	x1 = x[1:7,:]
	x2 = x[8:14,:]
	x3 = x[15:21,:]
	beta_ss1 = y1\x1
	beta_ss2 = y2\x2
	beta_ss3 = y3\x3
	beta_ss = (beta_ss1 + beta_ss3 + beta_ss3)/3.0 # average the 3 estimators
	betas[i,2] = beta_ss[2]
end

histogram(betas[:,1], label="", nbins=30, show=true)
#savefig("efficiency1.svg")
histogram(betas[:,2], label="", nbins=30, reuse=false, show=true)
#savefig("efficiency2.svg")
return
end
main()

max_issues_repo_name: DawievLill/Econometrics

max_forks_repo_name: DawievLill/Econometrics


text:
export ChebPEP

# Returns the chebyshev nodes scaled to interval [a,b]
# computed with type T
function get_chebyshev_nodes(::Type{T},a,b,k) where {T<:Real}
    mypi=T(pi);
    return (T(a)+T(b))/2 .+ (T(b)-T(a))*(cos.((2*Vector(1:k).-1)*mypi/(2*k)))/2
end

import Base.acos
function acos(S::LowerTriangular)
    # Specialized acos-function for a lower triangular matrix
    # Some extra allocations, but preserves realness.
    #
    # This is a work-around julia issue:
    #    https://github.com/JuliaLang/julia/issues/32721
    F=acos(Matrix(S));
    if (all(isreal.(acos.(diag(S)))))
        F .= real(F); # Try to enforce realness
    end
    return F
end
function cheb_f_cosine_formula(a,b,x,k)
    x1=2*(x-a*one(x))/(b-a)-one(x);
    if (istril(x1) && (x1 isa AbstractMatrix))
        x1=LowerTriangular(x1);
    end
    return cos(k*acos(x1));
end
function cheb_f_poly(a,b,x,k)
    x1=2*(x-a*one(x))/(b-a) -one(x);
    if (k==0)
        return one(x1);
    elseif k==1
        return x1
    elseif k==2
        return 2*x1^2-one(x1)
    elseif k==3
        return 4*x1^3-3*x1
    elseif k==4
        return 8*x1^4-8*x1^2+one(x1)
    elseif k==5
        return 16*x1^5 - 20*x1^3 + 5*x1;
    elseif k==6
        return 32*x1^6-48*x1^4+18*x1^2-one(x1);
    elseif k==7
        return 64*x1^7-112*x1^5+56*x1^3-7*x1;
    elseif k==8
        return 128*x1^8-256*x1^6+160*x1^4-32*x1^2+one(x1);
    else
        error("Not implemented")
    end
end

function cheb_f(a,b,x,k)
    return cheb_f_cosine_formula(a,b,x,k);
end

# Evaluate F in the chebyshev nodes
function chebyshev_eval(a,b,k,F::Function)
    x=get_chebyshev_nodes(Float64,a,b,k)
    Fk=F.(x);
    return (Fk,x)
end

function chebyshev_compute_coefficients_naive(a,b,Fk,k)
    xk=get_chebyshev_nodes(Float64,a,b,k)
    @show xk
    if (size(xk) != size(Fk))
        error("Incompatible sizes");
    end

    Tinv=zeros(k,k);
    for i=1:k
        for j=1:k
            Tinv[i,j]=cheb_f(a,b,xk[i],j-1);
        end
    end
    @show Tinv

    Ck=Tinv\Fk

    return Ck

end

# Compute the chebyshev coefficients of the coefficients
# stored in Fk. xk should be the chebyshev points
# Chebyshev Polynomials, 1st Edition, J.C. Mason, David C. Handscomb
# Chapter 8
function chebyshev_compute_coefficients(a,b,Fk,xk)
    # Return coefficient
    k=size(Fk,1);
    Tmat=zeros(k,k);
    for i=1:k
        # for each chebyshev polynomial, compute the k coefficients
        Tmat[i,:]= map(x->cheb_f(a,b,x,i-1)*2/k,xk)
    end
    Tmat[1,:] *= 0.5; # Fix the prime in the sum

    # Compute the coefficients from the Tmat via
    # a "matrix vector multiplication"

    # When Fk are scalars
    #Ck=Tmat*Fk;
    # Generalization when Fk are vectors or matrices:
    Ck=map(i-> sum(Fk .* Tmat[i,:]), 1:k)

    return Ck
end

struct ChebPEP{T<:AbstractMatrix,Ftype} <: AbstractSPMF{T}
    n::Int # size
    # Chebyshev polys are scaled to the interval [a,b]
    a::Ftype
    b::Ftype
    k::Int; # Number of Chebyshev polys
    spmf::SPMF_NEP{T,Ftype}; # The cheb-coefficents are stored directly in the SPMF
end

# Delegate compute functions

function size(nep::ChebPEP)
    return (nep.n,nep.n)
end
function size(nep::ChebPEP,dim)
    return nep.n
end
compute_Mlincomb(nep::ChebPEP,λ::Number,V::AbstractVecOrMat,a::Vector)= compute_Mlincomb(nep.spmf,λ,V,a);
compute_Mlincomb(nep::ChebPEP,λ::Number,V::AbstractVecOrMat)= compute_Mlincomb(nep.spmf,λ,V);
compute_Mder(nep::ChebPEP,λ::Number)=compute_Mder(nep.spmf,λ,0)
compute_Mder(nep::ChebPEP,λ::Number,der::Integer)=compute_Mder(nep.spmf,λ,der)
compute_MM(nep::ChebPEP,par...)=compute_MM(nep.spmf,par...)
get_Av(nep::ChebPEP)=get_Av(nep.spmf)
get_fv(nep::ChebPEP)=get_fv(nep.spmf)

"""
    ChebPEP(orgnep::NEP,k,[a=-1,[b=1]] [,cosine_formula_cutoff=5])

The type `ChebPEP<:AbstractSPMF` represents a polynomial function
where the function is stored using a Chebyshev basis scaled to the
interval `[a,b]`, i.e.,
```math
M(λ)= B_0T_0(λ)+⋯+B_{k-1}T_{k-1}(λ)
```
where ``T_i`` are the scaled and shifted Chebyshev polynomials.

The constructor `ChebPEP` takes `nep::NEP` as an input
and interpolates this NEP in `k` Chebyshev nodes, resulting
in a polynomial of degree `k-1`, represented by its
coefficients in the Chebyshev basis.
Interpolation in Chebyshev nodes and representation
with Chebyshev basis, is known to have
attractive approximation properties, as well
as robustness with respect to round-off errors.

The kwarg `cosine_formula_cutoff` decides how the Chebyshev
polynomials should be computed. For larger degrees, it is better
to use the cosine formula, whereas for low degrees
the explicit monomial expression is more efficient.
The explicit monomial expression will be used for degrees
lower than `cosine_formula_cutoff`.

# Example:

```julia
julia> nep=nep_gallery("dep0");
julia> chebpep=ChebPEP(nep,9);
julia> using LinearAlgebra;
julia> norm(compute_Mder(nep,0.3)-compute_Mder(chebpep,0.3))
1.2881862971045282e-8
julia> chebpep=ChebPEP(nep,19); # Better interpolation
julia> norm(compute_Mder(nep,0.3)-compute_Mder(chebpep,0.3))
2.0312004517316714e-15
```

See also: [`polyeig`](methods.md#NonlinearEigenproblems.NEPSolver.polyeig), [`PEP`](@ref)
"""
function ChebPEP(orgspmf,k,a=-1,b=1;cosine_formula_cutoff=5)
    F=s-> compute_Mder(orgspmf,s);
    (Fk,xk)=chebyshev_eval(a,b,k,F);
    Ck=chebyshev_compute_coefficients(a,b,Fk,xk);

    fv=Array{Function}(undef,0);
    for j=1:k
        if (j>cosine_formula_cutoff)
            push!(fv, S-> cheb_f_cosine_formula(a,b,S,j-1))
        else
            push!(fv, S-> cheb_f_poly(a,b,S,j-1))
        end
    end

    # Determine the types (based on Ck[1])
    MatType=typeof(Ck[1]);
    T=eltype(Ck[1]);

    # Create an SPMF with the correct type
    cheb_spmf=SPMF_NEP(Ck,fv,check_consistency=false,Ftype=T);

    n=size(Fk[1],1);
    # Instantiate the type
    return ChebPEP{MatType,T}(n,a,b,k,cheb_spmf);
end


text:
using PyCall

@pyimport gym

struct GymEnv <: AbstractEnvironment
    name::String
    action_space
    observation_space
    state::PyObject
    pygymenv
end

"""
    GymEnv(name)

Initialize a [OpenAIGym](https://github.com/openai/gym) environment
and wrap it in a `GymEnv` struct.
"""
function GymEnv(name::String)
    pygymenv = nothing
    try
        pygymenv = gym.make(name)
    catch e
        error("Error when initializing [$name]\n$e")
    end
    GymEnv(name,
           gymspace2jlspace(pygymenv[:action_space]),
           gymspace2jlspace(pygymenv[:observation_space]),
           PyNULL(),
           pygymenv)
end

"""
    step!(env::GymEnv, action)

`action` should be a valid sample in `env.action_space`.
"""
function step!(env::GymEnv, action)
    pycall!(env.state, env.pygymenv[:step], PyVector, action)
end

close!(env::GymEnv) = env.pygymenv[:close]()
seed!(env::GymEnv, seed=nothing) = env.pygymenv[:seed](seed)
reset!(env::GymEnv) = env.pygymenv[:reset]()
render(env::GymEnv, mode="human") = env.pygymenv[:render](mode)
action_space(env::GymEnv) = env.action_space
observation_space(env::GymEnv) = env.observation_space
get_state(env::GymEnv) = env.state[1]
get_reward(env::GymEnv) = env.state[2]
is_end(env::GymEnv) = env.state[3]


text:
function linsolve_bicgstab( LF::LF3dGrid, b::Array{Float64,1};
                            x0 = nothing,
                            NiterMax = 1000, TOL=5.e-10,
                            convmsg=false, showprogress=false )
    #
    Npoints = size(b)[1]
    #
    if x0 == nothing
        srand(1234)
        x = randn( Npoints )
    else
        x = copy(x0)
    end
    #
    r = b - apply_Laplacian(LF, x)
    #
    r_hat0 = copy(r)
    #
    rho = 1.0
    α   = 1.0
    ω   = 1.0
    #
    v = zeros(Float64,Npoints)
    p = zeros(Float64,Npoints)

    r_old = copy(r)
    v_old = zeros(Float64,Npoints)
    p_old = zeros(Float64,Npoints)
    x_old = copy(x)
    #
    rho_old = rho
    ω_old = ω
    for iter = 1:NiterMax
        rho = dot( r_hat0, r_old )
        β = (rho/rho_old)*(α/ω_old)
        p = r_old + β*( p_old - ω_old*v_old )
        v = apply_Laplacian(LF,p)
        α = rho/dot(r_hat0,v)
        h = x_old + α*p
        diffV_1 = norm((h-x))
        #@printf("iter, diffV_1 = %8d %e\n", iter, diffV_1)
        s = r_old - α*v
        t = apply_Laplacian(LF,s)
        ω = dot(t,s)/dot(t,t)
        x[:] = h[:] + ω*s[:]
        #
        diffV_2 = norm(x-x_old) #sum( abs(x-x_old) )
        @printf("iter, diffV_2 = %8d %e\n", iter, diffV_2)
        if diffV_1 < TOL || diffV_2 < TOL
            @printf("Convergence achieved in linsolve_bicgstab: %8d\n", iter)
            break
        end
        r = s - ω*t
        #
        rho_old = rho
        r_old = copy(r)
        v_old = copy(v)
        p_old = copy(p)
        x_old = copy(x)
    end
    return x
end


text:
@coeff function K_G(reactor::Reactor)
  0.3183 * N_G / ( epsilon ^ 2 )
end


text:
# Copyright (c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# Dynamic Bounds - pODEs Discrete
# A package for discretize and relax methods for bounding pODEs.
# See https://github.com/PSORLab/DynamicBoundspODEsDiscrete.jl
#############################################################################
# src/DiscretizeRelax/utilities/taylor_functor.jl
# Defines methods to used to compute Taylor coeffients.
#############################################################################
"""
TaylorFunctor!

A function g!(out, y) that perfoms a Taylor coefficient calculation. Provides
preallocated storage. Evaluating this function out is a vector of length nx*(s+1)
where 1:(s+1) are the Taylor coefficients of the first component, (s+2):nx*(s+1)
are the Taylor coefficients of the second component, and so on. This may be
constructed using `TaylorFunctor!(g!, nx::Int, np::Int, k::Val{K}, t::T, q::Q)`
were type `T` should use type `Q` for internal computations. The order of the
TaylorSeries is `k`, the right-hand side function is `g!`, `nx` is the number
of state variables, `np` is the number of parameters.

$(TYPEDFIELDS)
"""
mutable struct TaylorFunctor!{F <: Function, N, T <: Real, S <: Real}
    "Right-hand side function for pODE which operates in place as g!(dx,x,p,t)"
    g!::F
    "Dimensionality of x"
    nx::Int64
    "Dimensionality of p"
    np::Int64
    "Order of TaylorSeries, that is the first k terms are used in the approximation
    and N = k+1 term is bounded"
    k::Int64
    "State variables x"
    x::Vector{S}
    "Decision variables p"
    p::Vector{S}
    "Store temporary STaylor1 vector for calculations"
    xtaylor::Vector{STaylor1{N,S}}
    "Store temporary STaylor1 vector for calculations"
    xaux::Vector{STaylor1{N,S}}
    "Store temporary STaylor1 vector for calculations"
    dx::Vector{STaylor1{N,S}}
    taux::Vector{STaylor1{N,T}}
    vnxt::Vector{Int64}
    fnxt::Vector{Float64}
end

function (d::TaylorFunctor!{F, K, T, S})(out::Vector{Vector{S}}, x::Vector{S},
          p::Vector{S}, t::T) where {F <: Function, K, T <: Real, S <: Real}

    val = Val(K-1)
    for i in eachindex(d.xtaylor)
        d.xtaylor[i] = STaylor1(x[i], val)
    end

    jetcoeffs!(d.g!, t, d.xtaylor, d.xaux, d.dx, K-1, p, d.vnxt, d.fnxt)::Nothing
    for i in eachindex(out)
        for j in eachindex(d.xtaylor)
            out[i][j] = d.xtaylor[j][i-1]
        end
    end

    return nothing
end

function TaylorFunctor!(g!, nx::Int, np::Int, k::Val{K}, t::T, q::Q) where {K, T <: Number, Q <: Number}

    x0 = zeros(T, nx)
    Vⱼ = zeros(T, nx)
    f̃ = Vector{T}[]
    for i = 1:(K+1)
        push!(f̃, zeros(T, nx))
    end
    temp = STaylor1(zeros(T,K+1))
    xtaylor = STaylor1{K+1,T}[]
    xaux = STaylor1{K+1,T}[]
    dx = STaylor1{K+1,T}[]
    taux = STaylor1{K+1,Q}[]
    for i = 1:nx
        push!(xtaylor, temp)
        push!(xaux, temp)
        push!(dx, temp)
        push!(taux, zero(STaylor1{K+1,Q}))
    end
    #push!(xtaylor, fill(temp, nx))
    x = zeros(T, nx)
    p = zeros(T, np)
    vnxt = zeros(Int, nx)
    fnxt = zeros(Float64, nx)

    return TaylorFunctor!{typeof(g!), K+1, Q, T}(g!, nx, np, K, x, p, xtaylor,
                                                 xaux, dx, taux, vnxt, fnxt)
end


text:
function LineArc(line::Geometry, arc::Geometry)
    if line.GType != "line" || arc.GType != "arc"
        error("input type incorrect: line and arc needed")
    end
    Center = ArcCenter(arc)
    SymPy.@vars x y t
    x = line.Sp.x + t*(line.Ep.x - line.Sp.x)
    y = line.Sp.y + t*(line.Ep.y - line.Sp.y)
    Eq = (x - Center.x)^2 + (y - Center.y)^2 - (1/arc.Sp.cur)^2
    Roots = SymPy.solve(Eq)
    #Ips = Array{Tuple{Point2D{Float64},Float64, Float64}, 1}()
    #if !isa(SymPy.N(Roots[1]), Real)  # line and circle do not Intersect
        #println(typeof(Roots[1]))
    #    return (intersect=false, type="none")
    #else
    if isa(SymPy.N(Roots[1]), Real)
        for ti in Roots
            if 0.0 <= ti <= 1.0 + 1.0e-6  # check if intersection point is in line segment
                Ipc = Point2D{Float64}(SymPy.subs(x, t, ti), SymPy.subs(y, t, ti))
                dist1 = SymPy.N(ti*line.Length)
                Result = PointInSector(Ipc, arc)
                if Result[1] # check if intersection point is in arc
                    dist2 = Result[2]
                    return (intersect=true, type="cross", Ip=Ipc, Dist1=dist1, Dist2=dist2) # only the first junction point is returned
                    #push!(Ips, (Ipc, Dist1, Dist2))
                end
            end
        end
    end

    return (intersect=false, type="none")
end


text:
# header.jl

using Test: @test, @inferred

using MIRTio: header_read, header_write
using MIRTio: header_example, header_init, header_size, header_string

hd, hs = header_example()
@test header_size(hd) == hs
ht = header_init(hd) # random data for testing

tname = tempname()
open(tname, "w") do fid
    header_write(fid, ht ; bytes = header_size(hd))
end

open(tname, "r") do fid
    hr = header_read(fid, hd ; seek0 = true)
    @test hr == ht
    hu = header_string(hr)
    @test hu isa NamedTuple
    @test hu.date isa String
    @test all([hu[i] == ht[i] for i=2:length(hu)])
end


text:
"""
    dataset(name::String)::DataFrame

Load the dataset associated with `name`.
"""
function dataset(name::String)::DataFrame

    filename = joinpath(dirname(@__FILE__), "..", "data", string(name, ".csv"))
    if isfile(filename)
        return DataFrame(CSV.File(filename))
    end
    error("There is no dataset with the name '$name'")
end

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
#struct approximator end

mutable struct accumulator end
#<: approximator end

# enum op +, -,*,/

@enum op begin
    Add = 1
    Sub = 2
    Mul = 3
    Div = 4
end

function consume(op, a::Float64, b::Float64)
    if op isa Add()
        addCall(a::Float64, b::Float64)
    elseif op isa Sub
        #subCall()
    elseif op isa Mul #enought- for a banach space operations (covers  rationals  (real )
    elseif op isa Div # unnecessary for this context
    end
end
#TODO: Rethink a generalization about the 4 ops:
"""
we can never be sure about the Input Arguments: a  & b
in this file in particular, we cannot preallocate static types for arguments, then
unless, once we specify one method with particular set of `preallocated` argumnets
then:
      via multi-dispatch, we are constrained by defining each and every possible method there could exist
between the inputs a & b

returns the Sum ∑ #which sum?

"""
function addCall(a::Array{Float64,1}, b::Array{Float64,1}) #vectors a, b
    return Expr(:call, :+, a, b)  # a + b # do vector multiplication  #TODO:
end

function addCall(a::Array{Float64,2}, b::Array{Float64,2}) #matrix  a, b
    return # do matrix multiplication
end

function subCall(a, b)
    return a - b
end

function mulCall(a, b)
    return a * b
end

function divCall(a, b) end

#end

function arrayOp(; op = (:+, :-, :*, :/), a::Array{Float64,1}, b::Array{Float64,1}) #vectors a, b
    return Expr(:call, op, a, b)  # a + b # do vector multiplication  #TODO:
end  # #TODO: matrix Addition
function addCall(a::Array{Float64,1}, b::Array{Float64,1}) #vectors a, b
    return Expr(:call, :+, a, b)  # a + b # do vector multiplication  #TODO:
end
a = rand(10, 10);
b = rand(10, 10);

function arrayOp(;op=:+, a::Matrix, b::Matrix)
end

arrayOp()#:+, a, b)
   # addCall( )

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
resDASSL = Dict{Symbol,Any}()

dassl = IVPTestSuite.Solvers.allsolvers[DASSL.dasslSolve]
for (n,tc) in totest
    suite = TestSuite(tc, dassl, abstols, reltols, [NaN])
    resDASSL[n] = run_ode_testsuite(suite)
end


text:
using PkgUtility
using Test
using WaterPhysics

@testset "WaterPhysics --- functions" begin
    for FT in [Float32, Float64]
        rand_r = (rand(FT) + 20) * FT(1e-6);
        rand_T = rand(FT) + 298;
        rand_α = rand(FT) * 50;
        rand_Ψ = rand(FT) - 3;

        gas_air = TraceGasAir{FT}();
        gas_CO₂ = TraceGasCO₂{FT}();
        gas_O₂  = TraceGasO₂{FT}();
        liq_H₂O = TraceLiquidH₂O{FT}();

        # capillary_pressure
        for result in [ capillary_pressure(rand_r, rand_T),
                        capillary_pressure(rand_r, rand_T, rand_α) ]
            @test FT_test(result, FT);
            @test NaN_test(result);
        end

        # diffusive_coefficient
        for result in [ diffusive_coefficient(rand_T, gas_CO₂, gas_air),
                        diffusive_coefficient(rand_T, gas_CO₂, liq_H₂O),
                        relative_diffusive_coefficient(rand_T),
                        relative_diffusive_coefficient(rand_T, gas_CO₂, liq_H₂O) ]
            @test FT_test(result, FT);
            @test NaN_test(result);
        end
        @test relative_diffusive_coefficient(FT(298.15)) ≈ 1;

        # latent_heat_vapor
        for result in [ latent_heat_vapor(rand_T) ]
            @test FT_test(result, FT);
            @test NaN_test(result);
        end

        # surface_tension and relative_surface_tension
        for result in [ surface_tension(rand_T),
                        relative_surface_tension(rand_T) ]
            @test FT_test(result, FT);
            @test NaN_test(result);
        end
        @test relative_surface_tension(FT(298.15)) ≈ 1;

        # saturation_vapor_pressure
        for result in [ pressure_correction(rand_T, rand_Ψ),
                        saturation_vapor_pressure(rand_T),
                        saturation_vapor_pressure(rand_T, rand_Ψ),
                        saturation_vapor_pressure_slope(rand_T),
                        saturation_vapor_pressure_slope(rand_T, rand_Ψ) ]
            @test FT_test(result, FT);
            @test NaN_test(result);
        end

        # viscosity and relative_viscosity
        for result in [ viscosity(rand_T),
                        relative_viscosity(rand_T) ]
            @test FT_test(result, FT);
            @test NaN_test(result);
        end
        @test relative_viscosity(FT(298.15)) ≈ 1;
    end
end


text:
@testset ExtendedTestSet "Metric function" begin
    @testset ExtendedTestSet "`accuracy`" begin
        @test accuracy([1 0; 0 1], [1 0; 0 1]) == 1
        @test accuracy([0.9 0.1; 0.1 0.9], [1 0; 0 1]) == 1
    end

end


text:
@testset "Print" begin
    buf = IOBuffer()
    print(buf, ReferenceSequence(dna""))
    @test String(take!(buf)) == ""

    print(buf, ReferenceSequence(dna"ACGTN"))
    @test String(take!(buf)) == "ACGTN"

    print(buf, ReferenceSequence(dna"A"^100))
    @test String(take!(buf)) == "A"^100
end

@testset "Show" begin
    buf = IOBuffer()
    show(buf, ReferenceSequence(dna""))
    @test String(take!(buf)) == "< EMPTY SEQUENCE >"

    show(buf, ReferenceSequence(dna"ACGTN"))
    @test String(take!(buf)) == "ACGTN"
end


text:
module RFIOD
export binarysearch, fixedpointiteration, staircasediagram

using Plots

"""
    binarysearch(f::Function, a::Real, b::Real, [ε::Real])

Use the binary search algorithm to find a root of `f` between `a` and `b` with maximum trucation error `ε`.

Also returns the number of iterations that the algorithm takes.

# Examples

```julia-repl
julia> binarysearch(sin,3,4)
(3.141590118408203, 17)

julia> abs(ans[1] - π) < 5e-6
true
```
"""
function binarysearch(f::Function, a::Real, b::Real, ε::Real = 5e-6)

    # Checks that the initial interval is valid
    fa, fb = f(a), f(b)
    iszero(fa) && return (a,0)
    iszero(fb) && return (b,0)
    fa * fb < 0 || error("cannot find a root, function must take opposite signs on each end of the interval")

    # Iteratively finds the halfway point of the interval, then keeps the half where the sign of f changes
    N = 0
    while abs(b-a) > 2ε
        N += 1
        m = (a+b)/2
        fm = f(m)
        iszero(fm) && return (m,N)
        (a, b, fa, fb) = ( fa * fm > 0 ? (m, b, fm, fb) : (a, m, fa, fm) )
    end
    # Returns the midpoint of the interval, which is within ε of the root
    return ((a+b)/2,N)

end

"""
    fixedpointiteration(f::Function, x₀::Real, Nₘₐₓ::Int, [ε::Real])

Use the fixed point iteration algorithm to converge on a fixed point of `f`, with starting value `x₀`.

Terminate if consecutive iterations differ by less than `ε`, or after `Nₘₐₓ` iterations.

Returns the sequence resulting from the iteration.

# Examples

```julia-repl
julia> xs = fixedpointiteration(cos, 0, 100)
92-element Vector{Float64}:
 0.0
 1.0
 [...]
 0.7390851332151607

julia> cos(xs[end])
0.7390851332151607
```
"""
function fixedpointiteration(f::Function, x₀::Real, Nₘₐₓ::Int, ε::Real = eps(Float64))

    # Sets up the initial values for the algorithm
    n = 1
    xₙ₋₁ = convert(Float64,x₀)
    xₙ = f(xₙ₋₁)
    xvalues = [xₙ₋₁, xₙ]

    # Iteratively advances the algorithm by evaluating f at the previous result and storing it
    while abs(xₙ₋₁ - xₙ) ≥ ε && n < Nₘₐₓ
        n += 1
        xₙ₋₁, xₙ = xₙ, f(xₙ)
        push!(xvalues,xₙ)
    end
    # Returns the sequence resulting from the algorithm
    return xvalues

end

"""
staircasediagram(f::Function, xvalues::Vector{Float64})

Create a staircase/cobweb diagram from `xvalues`, the results of the fixed point iteration algorithm on `f`.

# Examples

```julia-repl
julia> staircasediagram(cos, fixedpointiteration(cos, 0, 100))
```
"""
function staircasediagram(f::Function, xvalues::Vector{Float64})

    # Finds the scale at which to draw the diagram
    minx, maxx = min(xvalues...), max(xvalues...)
    δ = (maxx - minx)/100
    xticks = (minx - 10δ):δ:(maxx + 10δ)

    # Plots y = x and y = f(x), then adds in a marker for the starting point and the lines given by the sequence
    p = plot([x -> x, f], xticks, legend = false)
    scatter!(p, (xvalues[1], xvalues[1]), markercolor = :black)
    for n = 1:(length(xvalues)-1)
        plot!(p, [xvalues[n], xvalues[n], xvalues[n+1]],
            [xvalues[n], xvalues[n+1], xvalues[n+1]], linecolor = :black)
    end
    # Returns the plot
    return p

end

end


text:
## PART 1

include("partition.jl");
include("../q1/autocorrelation.jl");

println("\n\tPick ncut = 150, as before");
ncut = 150;

# Partition the set into M/N bins
println("\n\tPartitioning v into M/N bins for N1 = $N1 and N2 = $N2");
vpartsN1 = partition(v, M, N1);
vpartsN2 = partition(v, M, N2);

println("\n\tCounting number of bins");
npartsN1 = size(vpartsN1)[3];
npartsN2 = size(vpartsN2)[3];

println("\n\tCalculating τint for M/N1 bins");
τintN1 = zeros(npartsN1,5);
for npart in [1:npartsN1]
    vpart = vpartsN1[:,:,npart];
    τintN1[npart,:] = τ(N1, vpart, ncut);
end

println("\n\tCalculating στint for M/N1 bins");
μτintN1 = zeros(5);
στintN1 = zeros(5);
for i in [1:5]
    μτintN1[i] = mean(τintN1[:,i]);
    στintN1[i] = std(τintN1[:,i]);
end

println("\n\tCalculating τint for M/N2 bins");
τintN2 = zeros(npartsN2,5);
for npart in [1:npartsN2]
    vpart = vpartsN2[:,:,npart];
    τintN2[npart,:] = τ(N2, vpart, ncut);
end

println("\n\tCalculating στint for M/N2 bins");
μτintN2 = zeros(5);
στintN2 = zeros(5);
for i in [1:5]
    μτintN2[i] = mean(τintN2[:,i]);
    στintN2[i] = std(τintN2[:,i]);
end


text:
infile = "./src/Day11/input.txt"
input = readlines(infile)

seating  = Matrix{String}(undef, length(input), length(input[1]))
for line in 1:length(input)
    s = input[line]
    st = string.(split(s, ""))
    for chr in 1:length(st)
        seating[line, chr] = st[chr]
    end
end

seatmask = [I for I in CartesianIndices((-1:1, -1:1)) if I != CartesianIndex(0,0)]

function update(seating, seatmask)
    newseating = copy(seating)
    neighbours = zeros(Int64, length(input), length(input[1]))
    for I in CartesianIndices(seating)
        for x in seatmask
            if checkbounds(Bool, seating, I + x)
                if seating[I + x] == "#"
                    neighbours[I] += 1
                end
            end
        end
        if seating[I] == "#" && neighbours[I] >= 4
            newseating[I] = "L"
        elseif seating[I] == "L" && neighbours[I] == 0
            newseating[I] = "#"
        end
    end
    return newseating
end

current = copy(seating)
function runtothestables(current, seatmask, maxattempts)
    for i in 1:maxattempts
        current, prev = update(current, seatmask), current
        if all(current .== prev)
            return current, i
        end
    end
    error("Need to increase maxattempts")
end

a,b = runtothestables(current, seatmask, 100)
count(x -> (a[x] == "#"), CartesianIndices(a))

function checkvisibleneighbours(seating, seatmask, pos)
    neighbours = String[]
    for dir in seatmask
        for steps in 1:maximum(size(seating))
            if checkbounds(Bool, seating, pos + steps*dir)
                if seating[pos + steps*dir] != "."
                    push!(neighbours,  seating[pos + steps*dir])
                    break
                end
            else
                break
            end
        end
    end
    return neighbours
end

function part2update(seating, seatmask)
    newseating = copy(seating)
    for I in CartesianIndices(seating)
       n = checkvisibleneighbours(seating, seatmask, I)
       if seating[I] == "#" && count(==("#"), n) >= 5
            newseating[I] = "L"
        elseif seating[I] == "L" && count(==("#"), n) == 0
            newseating[I] = "#"
        end
    end
    return newseating
end

a = part2update(seating, seatmask)

current = copy(seating)
function runtothestables2(current, seatmask, maxattempts)
    for i in 1:maxattempts
        current, prev = part2update(current, seatmask), current
        if all(current .== prev)
            return current, i
        end
    end
    error("Need to increase maxattempts")
end

a, b = runtothestables2(current, seatmask, 100)
count(x -> (a[x] == "#"), CartesianIndices(a))


text:
# This file is a part of DataDeps.jl. License is MIT.

const registry = Dict{String, AbstractDataDep}()

"""
    register(datadep::AbstractDataDep)

Registers the given datadep to be globally available to the program.
this makes `datadep"Name"` work.
`register` should be run within this `__init__` of your module.
"""
function register(datadep::AbstractDataDep)
    name = datadep.name
    if haskey(registry, name)
        @warn("Over-writing registration of the datadep", name)
    end
    if !is_valid_name(name)
        throw(ArgumentError(name, " is not a valid name for a datadep. Valid names must be legal foldernames on Windows."))
    end

    registry[name] = datadep
end

"""
    is_valid_name(name)

This checks if a datadep name is valid.
This basically means it must be a valid folder name on windows.

"""
function is_valid_name(name)
    namechars = collect(name)
    !any( namechars .∈ "\\/:*?<>|") && !any(iscntrl.(namechars))
end


text:
using RDKafka
using Test

a = RDKafka.rd_kafka_version() # 17105151 or 0x010500ff for 1.0.5
# verify the major and minor. Any change to these should be reflected
# by a change in Project.toml
@test (a & 0x01000000) == 0x01000000
@test (a & 0x00050000) == 0x00050000

conf =  RDKafka.kafka_conf_new()
@test RDKafka.kafka_conf_get(conf, "socket.keepalive.enable") == "false"
RDKafka.kafka_conf_set(conf, "socket.keepalive.enable", "true")
@test RDKafka.kafka_conf_get(conf, "socket.keepalive.enable") == "true"

# Verify that uninitalised memory is not exposed to caller
@test RDKafka.kafka_conf_get(conf, "foobar") == ""

RDKafka.kafka_conf_destroy(conf)
@test true # No exceptions


text:
using Plots
using LinearAlgebra

using SparseArrays
using Arpack
#komvos1 komvos2 komvos3
xy = [ 0 0;
       0 1]; #orismos komvon san xyz
nodes = size(xy)[2]; #arithmos komvon

#komvos1_start komvos1_end
numel = [1 2;

]; #orismos ton stixion os grammes me toys akraious komvous
N = size(numel)[1]; #arithmos stoixion
dofs = 3*nodes;
#orismos elastikotitas stixion
a = 10e-2;
E = 220e9.*ones(N,1); I = a^4 .*ones(N,1)./12; A = a^2 .*ones(N,1) #is square with axa
rho = 7.8e3;

#orismos ton desmeumenon komvon
BC = zeros(dofs,1)
#ignored = 1:3:dofs;#ignore turn of axis=ignore freedom
BC[1:3] .= 1; #BC[ignored].=1;
fixedofs = findall(x->x==1,BC[:])
freedofs = findall(x->x==0,BC[:])
#boundary conditions
U = zeros(dofs,1);
F = zeros(dofs,1); F[[5]] .= 10e3;
#euresi olon ton deikton pou tha prospelasoun ton K
alldofs = zeros(N,6);
#orismos olikou mitroou stivarotitas
global K = sparse([],[],Float64[],dofs,dofs);
global M = sparse([],[],Float64[],dofs,dofs);

#global K = zeros(dofs,dofs);
#start main loop to calculate stiffness matrix
for i=1:N
    i1 = numel[i,1]
    i2 = numel[i,2]
    #calclulate directions
    x1 = xy[1,i1]; x2 = xy[1,i2]; dx = x2-x1;
    y1 = xy[2,i1]; y2 = xy[2,i2]; dy = y2-y1;
    L = sqrt(dx^2+dy^2);
    nx = dx/L; ny = dy/L;
    #direction matrix
    c = nx; s = ny;
    Bbend = [-s c 0  0 0 0;
              0 0 1  0 0 0;
              0 0 0 -s c 0;
              0 0 0  0 0 1];
    Btens = [c s 0 0;
             0 0 c s];
    kbend = E[i]*I[i]/L^3;

    #insert tensile
    ktens = E[i]*A[i]/L
    global Ketens = ktens*[1 -1; -1 1]
    Ketens = Btens'*Ketens*Btens
    global Metens = rho*A[i]*L*[2 1; 1 2]./6;
    Metens = Btens'*Metens*Btens;

    Kebend = kbend*[12 6*L   -12   6*L;
                     0 4*L^2 -6*L 2*L^2;
                     0  0     12  -6*L;
                     0  0     0   4*L^2]
    Kebend = Symmetric(Kebend);
    Kebend = Bbend'*Kebend*Bbend;
    Mebend = rho*A[i]*(L/420)*[156   22*L  54     -13*L;
                                0   4*L^2 13*L   -3*L^2;
                                0     0   156     -22*L;
                                0     0    0     4*L^2];
    Mebend = Symmetric(Mebend);
    Mebend = Bbend'*Mebend*Bbend;

    tensiledofs = [1 2 4 5];
    Ke = Kebend .+0;
    Ke[tensiledofs,tensiledofs] = Ke[tensiledofs,tensiledofs][1,:,1,:]  + Ketens;
        #println(size(Ke),"\n")
    edofs = [3*i1-2 3*i1-1 3*i1  3*i2-2 3*i2-1 3*i2]
    K[edofs,edofs] = K[edofs,edofs][1,:,1,:] + Ke;
    #println(K,"\n")
    Me = Mebend .+0;
    Me[tensiledofs,tensiledofs] = Me[tensiledofs,tensiledofs][1,:,1,:]  + Metens;
    M[edofs,edofs] = M[edofs,edofs][1,:,1,:] + Me;
end
Kfinal = K[freedofs,freedofs]
Mfinal = M[freedofs,freedofs]
rhs = F
rhs = rhs[freedofs]
rhs = rhs - K[freedofs,fixedofs]*U[fixedofs]

Ufree = Kfinal\rhs
U[freedofs] = Ufree
U = reshape(U,(3,nodes))
show(U[:])
#println(varinfo())
omega = sqrt.(eigvals(inv(Matrix(Mfinal))*Matrix(Kfinal)))/(2*pi)
scatter(omega)

#calculate

max_stars_repo_path: axnai julia/test2.jl

max_stars_repo_name: billyon/first

max_issues_repo_path: axnai julia/test2.jl

max_issues_repo_name: billyon/first

max_forks_repo_path: axnai julia/test2.jl

max_forks_repo_name: billyon/first


text:
using SparseArrays
using GenericLinearAlgebra
export presolve

"""
    presolve(optimizer::SDPA_GMP.Optimizer{T}) where T

Identifies linearly dependent constraints in the problem. This is done by a naive Gaussian elimination.

Returns a vector with the indices of redundant constraints, which should be removed from the formulation. The rest of the constraints form a maximal linearly independent subset of the original set of constraints.
"""
function presolve(optimizer::SDPA_GMP.Optimizer{T}) where T
    start = time()
    totaldim = sum(abs.(optimizer.blockdims))
    F = spzeros(T, length(optimizer.b), totaldim^2)
    cVec = optimizer.b
    for entry in optimizer.elemdata
        constr_index, blk, i, j, value = entry
        if constr_index != 0
            offset = sum(abs.(optimizer.blockdims)[1:blk-1])
            linear_index = (j + offset - 1) * totaldim + (i + offset)
            F[constr_index, linear_index] = value
        end
    end
    aug_mat = hcat(F, cVec)
    redundant_F = collect(setdiff!(
        Set(1:length(optimizer.b)),
        Set(rowvals(reduce!(aug_mat)))
    ))
    for i in redundant_F
        if aug_mat[i, end] != 0
            error("Inconsistency at $i th constraint. Problem is dual infeasible. ")
        end
    end
    finish = time() - start
    n = length(redundant_F)
    @info "Presolve finished in $finish seconds. $n constraint(s) eliminated."
    return sort!(redundant_F)
end

"""
    function reduce!(A::SparseMatrixCSC{T}, ɛ = T <: Union{Rational,Integer} ? 0 : eps(norm(A, Inf))) where T

Identifies linearly dependent constraints in the problem. This is done by a naive Gaussian elimination.

Returns a vector with the indices of redundant constraints, which should be removed from the formulation. The rest of the constraints form a maximal linearly independent subset of the original set of constraints.
"""
function reduce!(A::SparseMatrixCSC{T}, ɛ = T <: Union{Rational,Integer} ? 0 : eps(norm(A, Inf))) where T
    nr, nc = size(A)
    i = j = 1
    visited_rows = Int[]
    while i <= nr && j <= nc - 1 # avoid touching the cVec
        rows = setdiff!(rowvals(A)[nzrange(A, j)], visited_rows)
        if isempty(rows)
            j += 1
            continue
        end
        (m, mi) = findmax(abs.([A[p, j] for p in rows]))
        mi = rows[mi]
        if m <= ɛ
            if ɛ > 0
                A[:, j] .= zero(T)
                dropzeros!(A)
            end
            j += 1
        else
            push!(visited_rows, mi)
            d = A[mi, j]
            A[mi, j:nc] = A[mi, j:nc] / d

            for k in setdiff!(rows, mi)
                d = A[k, j]
                A[k, j:nc] -= d * A[mi, j:nc]
            end
            dropzeros!(A)
            i += 1
            j += 1
        end
    end
    return droptol!(A, ɛ)
end


text:
using aoc2019.Nanofactory

input = joinpath(@__DIR__, "input")
eqs = readlines(input)

formulae = Nanofactory.parse_equations(eqs)

inventory = Dict("FUEL" => -1)
Nanofactory.produce(inventory, formulae)
p1 = abs(inventory["ORE"])
@assert p1 == 201324

low = 1
high = 1000000000000
while low < high
    global low
    global high
    guess = div(low + high + 1, 2)
    inventory = Dict("FUEL" => -1 * guess)
    Nanofactory.produce(inventory, formulae)
    ore = abs(inventory["ORE"])
    if ore < 1000000000000
        low = guess
    else
        high = guess - 1
    end
end

p2 = low
@assert p2 == 6326857

print("-----------------------------------------------------------------------\n")
print("space stoichiometry -- part one\n    ore  : $p1\n")
print("space stoichiometry -- part two\n    fuel : $p2\n")
print("-----------------------------------------------------------------------\n")


text:
function create_index(client::HODClient, flavor = "standard", index_fields = "",
                      parametric_fields = "")

end

function parse_index(client::HODClient, index::HODIndex)

end

function has_index(client::HODClient, name::UTF8String)

end

function get_index(client::HODClient, name::UTF8String)

end

function delete_index(client::HODClient, name)

end

function list_indexes(client::HODClient, field_type = "", flavor = "standard")

end

function post(client::HODClient, handler, data = Dict(), files = Dict(),
              async = false; kwargs...)

end

function post_async(client::HODClient, handler, data = Dict(), files = Dict();
                    kwargs...)

end

function call_API(client::HODClient, url, data = Dict(), files = Dict(),
                 async = false)

end

function get_async_result(client::HODClient, job_id)

end


text:
export
    CUDNN_BATCHNORM_PER_ACTIVATION, CUDNN_BATCHNORM_SPATIAL # mode

function batchnorm_inference(mode, x, scale, bias, estimated_mean, estimated_var, epsilon::Float64)
    T = eltype(x)
    h = handle(x)
    y = similar(x)
    xdesc = tensor_desc(x)
    ydesc = tensor_desc(y)
    sbmvdesc = tensor_desc(scale)

    cudnnBatchNormalizationForwardInference(h, mode, T[1], T[0], xdesc, x,
    ydesc, y, sbmvdesc, scale, bias, estimated_mean, estimated_var, Cdouble(epsilon))

    cudnnDestroyTensorDescriptor(xdesc)
    cudnnDestroyTensorDescriptor(ydesc)
    cudnnDestroyTensorDescriptor(sbmvdesc)
    y
end

function batchnorm_training(mode, x, scale, bias, averagefactor::Float64, epsilon::Float64)
    T = eltype(x)
    h = handle(x)
    y = similar(x)
    xdesc = tensor_desc(x)
    ydesc = tensor_desc(y)
    sbmvdesc = tensor_desc(scale)
    running_mean = zeros(scale)
    running_var = zeros(scale)
    save_mean = similar(scale)
    save_invvar = similar(scale)

    cudnnBatchNormalizationForwardTraining(h, mode, T[1], T[0], xdesc, x, ydesc, y,
    sbmvdesc, scale, bias, Cdouble(averagefactor), running_mean, running_var,
	Cdouble(epsilon), save_mean, save_invvar)

    cudnnDestroyTensorDescriptor(xdesc)
    cudnnDestroyTensorDescriptor(ydesc)
    cudnnDestroyTensorDescriptor(sbmvdesc)
    y, running_mean, running_var, save_mean, save_invvar
end

function ∇batchnorm!(mode, x, dy, dx, scale, dscale, dbias, epsilon::Float64, saved_mean, saved_invvar)
    T = eltype(x)
    h = handle(x)
    xdesc = tensor_desc(x)
    dydesc = tensor_desc(dy)
    dxdesc = tensor_desc(dx)
    dsbdesc = tensor_desc(scale)

    cudnnBatchNormalizationBackward(h, mode, T[1], T[0], T[1], T[0],
    xdesc, x, dydesc, dy, dxdesc, dx, dsbdesc, scale, dscale, dbias,
    Cdouble(epsilon), saved_mean, saved_invvar)

    cudnnDestroyTensorDescriptor(xdesc)
    cudnnDestroyTensorDescriptor(dydesc)
    cudnnDestroyTensorDescriptor(dxdesc)
    cudnnDestroyTensorDescriptor(sbdesc)
end


text:
"""
    MultiDimEquations

The package provides handy functions to work with NDSparse (from IndexedTables.jl package)
or standard arrays in order to write equations in a concise and readable way.
The formers can be accessed by name, but are somehow slower, the latters are faster
but need to be accessed by position.

`defVars()` defines either empty NDSparse with the dimensions (axis) required or
arrays with the required size and prefilled with `missing` values. More variables
can be defined at once.

`defLoadVar()` and `defLoadVars()` allow to define and import the variables from
a DataFrame in long format (dim1|dim2|...|value), with the second one allowing to
import more varaibles at once, given the presence of a column in the database with
the variable names.

`@meq` is a macro that allows to write your "model equations" more concisely, e.g.
`@meq par1[d1 in DIM1, d2 in DIM2, dfix3] =  par2[d1,d2]+par3[d1,d2]` would be
expanded to use the list comprehension expression above:
`[par1[d1,d2,dfix3] =  par2[d1,d2]+par3[d1,d2] for d1 in DIM1, d2 in DIM2]`

"""
module MultiDimEquations

using DataFrames, IndexedTables

export defLoadVars, defLoadVar, defVars, @meq, getSafe, toDataFrame

##############################################################################
#
# defVars()
#
##############################################################################

"""
    defVars(dimNames, dimTypes; <keyword arguments>)

Define empty NDSparse IndexedTable(s) with the specific dimension(s) and type.

# Arguments
* `dimNames`: Array of names of the dimensions to define
* `dimTypes`: Array of types of the dimensions
* `valueType` (def: `Float64`):  Type of the value column of the table
* `n` (def=`1`): Number of copies of the specified tables to return (useful to define multiple variables at once. In such cases a tuple is returned)

# # Examples
# ```julia
# julia> price,demand,supply = defVars(["region","item","class"],[String,String,Int64],valueType=Float64,n=3 )
# julia> waterContent = defVars(["region","item"],[String,String])
# julia> price["US","apple",1] = 3.2
# julia> waterContent["US","apple"] = 0.2
# ```
# """
function defVars(dimNames, dimTypes; valueType=Float64,n=1)
      values = [Array{T,1}() for T in vcat(dimTypes,valueType)]
      if n==1
          return NDSparse(values...,names=Symbol.(dimNames))
      else
          return [deepcopy(NDSparse(values...,names=Symbol.(dimNames))) for i in 1:n]
      end
end

"""
    defVars(size; <keyword arguments>)

Define multidimensional array(s) with the specific dimension(s) and type filled all with `missing` values.

# Arguments
* `size`: Tuple of the dimensions required
* `valueType` (def: `Float64`):  Inner type of the array required
* `n` (def=`1`): Number of copies of the specified tables to return (useful to define multiple variables at once. In such cases a tuple is returned)
* `missingValue` (def=`missing`): How to fill the matrix with

# # Examples
# ```julia
# julia> price,demand,supply = defVars((3,4,5),valueType=Float64,n=3 )
# julia> waterContent = defVars((3,4))
# julia> price[2,3,1] = 3.2
# julia> waterContent[2,3] = 0.2
# ```
# """
function defVars(size; valueType=Float64,n=1,missingValue=missing)
    missingValueType = typeof(missingValue)
    val = fill(missingValue,size...)
    val = convert(Array{Union{missingValueType,valueType},length(size)},val)
    if n==1
        return val
    else
        return [deepcopy(val) for i in 1:n]
    end
end

##############################################################################
#
# defLoadVars()
#
##############################################################################

"""
    defLoadVar(df, dimsNameCols; <keyword arguments>)

Define the required IndexedTables or Arrays and load the data from a DataFrame in long format while specifing the dimensional columns.

# Arguments
* `df`: The source of the dataframe, that must be in the format dim1|dim2|...|value
* `dimsNameCols`: The names of the columns corresponding to the dimensions over which the variables are defined (the keys)
* `valueCol (def: "value")`: The name of the column in the df containing the values
* `sparse (def: "true")`: Wheter to return `NDSparse` elements (from IndexedTable) or standard arrays
* `missingValue` (def=`missing`): How to fill the matrix with (relevant only for arrays)

# Notes
* Sparse indexed tables can be accessed by element but are slower, standard arrays need to be accessed by position but are faster

# Examples
```julia
julia> vol  = defVars(volumeData,["region","treeSpecie","year"], valueCol="value")
```
"""
function defLoadVar(df, dimsNameCols; valueCol="value",sparse=true,missingValue=missing)
    valueType = eltype(df[!,valueCol])
    nDims     = length(dimsNameCols)
    if sparse
        sDimensions = [Symbol(d) for d in dimsNameCols]
        dimValues  = [df[:,dim] for dim in dimsNameCols]
        values     = df[:,valueCol]
        t = IndexedTables.NDSparse(dimValues..., names=sDimensions, values)
        return t
    else
        toReturn = Array[]
        dimItems = [unique(df[!,col]) for col in dimsNameCols]
        size = [length(dimItem_i) for dimItem_i in dimItems]
        varArray = defVars(size, valueType=valueType,n=1,missingValue=missingValue)

        for i in CartesianIndices(varArray)
           cIdx = Tuple(i)
           particularDims = [map(x->dimItems[d][x], cIdx[d]) for d in 1:nDims]
           selectionArray = fill(false,Base.size(df,1))
           for (i,row) in enumerate(eachrow(df))
               dimFilter = true
               for (d,dim) in enumerate(dimsNameCols)
                   if (row[dim] != particularDims[d])
                       dimFilter = false
                       break
                   end
               end
               selectionArray[i] = dimFilter
           end
           particularValueArray = df[selectionArray,valueCol]
           if(length(particularValueArray)>1)
               @error "In converting a long dataframe to an array, I found more than one record with the same keys."
               particularValue = missingValue
           elseif length(particularValueArray)==0 # this combination is not found
               particularValue = missingValue
           else
               particularValue = particularValueArray[1]
           end
           varArray[i]  = particularValue
        end
        varArray
    end
end

##############################################################################
#
# defLoadVars()
#
##############################################################################

"""
    defLoadVars(vars, df, dimsNameCols; <keyword arguments>)

Define the required IndexedTables or Arrays and load the data from a DataFrame in long format while specifing the dimensional columns and the column containing the variable names.

Like `varLoadVar` but here we are extracting multiple variables at once, with one column of the input dataframe storing the variable name.

# Arguments
* `vars`: The array of variables to lookup
* `df`: The source of the dataframe, that must be in the format dim1|dim2|...|varName|value
* `dimsNameCols`: The name of the column containing the dimensions over which the variables are defined
* `varNameCol (def: "varName")`: The name of the column in the df containing the variables names
* `valueCol (def: "value")`: The name of the column in the df containing the values
* `sparse (def: "true")`: Wheter to return `NDSparse` elements (from IndexedTable) or standard arrays
* `missingValue` (def=`missing`): How to fill the matrix with (relevant only for arrays)
* `fullKeys` (def=`true`): Wheter all the outputed variables should be based on the full set of keys,
   or look for each of them on the keys found specifically for it in the filtered input dataframe

# Notes
* Sparse indexed tables can be accessed by element but are slower, standard arrays need to be accessed by position but are faster
* At the moment `sameKeys=true` is implemented only for the non-sparse case

# Examples
```julia
julia> (vol,numberOfTrees)  = defVars(["vol","numberOfTrees"], forestData,["region","treeSpecie","year"], varNameCol="parName", valueCol="value")
```
"""
function defLoadVars(vars, df, dimsNameCols; varNameCol="varName", valueCol="value",sparse=true, missingValue=missing,fullKeys=true)
    if (!fullKeys) # each variable its own keys set, I can use defLoadVar
        if sparse
            toReturn = NDSparse[]
        else
            toReturn = Array[]
        end
        for var in vars
            filteredDf = df[df[!,varNameCol] .== var,:]
            varData = defLoadVar(filteredDf, dimsNameCols; valueCol=valueCol,sparse=sparse,missingValue=missingValue)
            if length(vars) > 1
                push!(toReturn,varData)
            else
                return varData
            end
        end
        return (toReturn...,)
    else # all output variables share the same set of keys, I can't reuse defLoadVar
        valueType = eltype(df[!,valueCol])
        nDims     = length(dimsNameCols)
        if sparse
            toReturn = NDSparse[]
            sDimensions = [Symbol(d) for d in dimsNameCols]
            for var in vars
                filteredDf = df[df[!,varNameCol] .== var,:]
                dimValues  = [filteredDf[:,dim] for dim in dimsNameCols]
                values     = filteredDf[:,valueCol]
                t = IndexedTables.NDSparse(dimValues..., names=sDimensions, values)
                if length(vars) > 1
                    push!(toReturn,t)
                else
                    return t
                end
            end
            return (toReturn...,)
        else
            toReturn = Array[]
            dimItems = [unique(df[!,col]) for col in dimsNameCols] # this should be general, i.e. independent on the specific variable we are looking at
            size = [length(dimItem_i) for dimItem_i in dimItems]
            for var in vars
                dfVar = df[df[!,varNameCol] .== var,:]
                varArray = defVars(size, valueType=valueType,n=1,missingValue=missingValue)
                for i in CartesianIndices(varArray)
                   cIdx = Tuple(i)
                   particularDims = [map(x->dimItems[d][x], cIdx[d]) for d in 1:nDims]

                   selectionArray = fill(false,Base.size(dfVar,1))
                   for (i,row) in enumerate(eachrow(dfVar))
                       dimFilter = true
                       for (d,dim) in enumerate(dimsNameCols)
                           if (row[dim] != particularDims[d])
                               dimFilter = false
                               break
                           end
                       end
                       selectionArray[i] = dimFilter
                   end
                   particularValueArray = dfVar[selectionArray,valueCol]
                   if(length(particularValueArray)>1)
                       @error "In converting a long dataframe to an array, I found more than one record with the same keys."
                       particularValue = missingValue
                   elseif length(particularValueArray)==0
                       particularValue = missingValue
                   else
                       particularValue = particularValueArray[1]
                   end
                   varArray[i]  = particularValue
                end
                push!(toReturn,varArray)
            end
            return (toReturn...,)
        end

    end
end

##############################################################################
#
# meq()
#
##############################################################################

"""
    meq(exp)

Macro to expand functions like `t[d1 in dim1, d2 in dim2, dfix,..] = value`

With this macro it is possible to write:

```
@meq par1[d1 in DIM1, d2 in DIM2, dfix3] =  par2[d1,d2]+par3[d1,d2]
```

and obtain
```
[par1[d1,d2,dfix3] =  par2[d1,d2]+par3[d1,d2] for d1 in DIM1, d2 in DIM2]
```

That is, it is possible to write "model equations" in a concise and readable way.

"""
macro meq(eq) # works without MacroTools
    lhs                   = eq.args[1]
    rhs                   = eq.args[2]
    lhs_par               = lhs.args[1]
    lhs_dims              = lhs.args[2:end]
    loop_counters         = [d.args[2] for d in lhs_dims if typeof(d) == Expr]
    loop_sets             = [d.args[2] for d in lhs_dims if typeof(d) == Expr]
    loop_wholeElements    = []
    lhs_dims_placeholders = []
    for d in lhs_dims
        if typeof(d) == Expr && d.args[1] == :in
            push!(lhs_dims_placeholders,d.args[2])
            push!(loop_wholeElements, :($(d.args[2]) = $(d.args[3])))
        else
            push!(lhs_dims_placeholders,d)
        end
    end
    if (lhs.head == :ref)      # lhs is an array
        ret = Expr(:comprehension, :($lhs_par[$(lhs_dims_placeholders...)] = $(rhs)),loop_wholeElements...)
#    elseif (lhs.head == :call) # lhs is a function call
#        ret = Expr(:comprehension, :($lhs_par($(rhs),$(lhs_dims_placeholders...))),loop_wholeElements...)
    else
        error("Didn't understand the Left Hand Side.")
    end
    #show(ret)
    return esc(ret)
end

##############################################################################
#
# getSafe()
#
##############################################################################

"""
    getSafe(idxtable,indices,missingValue=missing)

Return the value stored in a NDSParse table or missingValue if the specified keys are not present.

# Arguments
* `idxtable`: The NDSParse table to lookup
* `indices`: A tuple with the indices to use (`:` is supported)
* `missingValue`: The value to return if the specified keys are not found

# Examples
```julia
julia> volBlackForest2014  = getSafe(forestVolumes,("BlackForest",2014),0.0)
```
"""
function getSafe(idxtable,indices,missingValue=missing)
    try
        return idxtable[indices...]
    catch  e
        if isa(e, KeyError)
            return missingValue
        else
            rethrow(e)
        end
    end
end

##############################################################################
#
# toDataFrame()
#
##############################################################################

"""
    toDataFrame(idxtable::NDSparse)

Convert a single-value-column NDSParse table to `DataFrame`.

# Arguments
* `idxtable`: The NDSParse table to convert

# Examples
```julia
julia> content = [["banana","banana","apple","apple","orange"],["us",missing,"us","eu","us"],[1.1,2.2,3.3,missing,5.5]]
julia> dimNames = ["item","region"]
julia> t = NDSparse(content...,names=Symbol.(dimNames))
julia> df = toDataFrame(t)
```
"""
function toDataFrame(t::NDSparse)
    names = vcat(keys(keys(t)[1])...,:value)
    cols  = columns(t)
    return DataFrame(map((n,v) -> Pair(n,v), names, cols))
end

end # end module

max_forks_repo_name: sylvaticus/MultiDimEquations


text:
using Oceananigans
using Oceananigans.Units
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom
using Oceananigans.Models.HydrostaticFreeSurfaceModels: FFTImplicitFreeSurfaceSolver
using Printf

grid = RectilinearGrid(CPU(),
                       topology = (Periodic, Bounded, Bounded),
                       size = (128, 128, 16),
                       x = (-500kilometers, 500kilometers),
                       y = (-500kilometers, 500kilometers),
                       z = (-4kilometers, 0),
                       halo = (3, 3, 3))

name = @sprintf("baroclinic_adjustment_Nx%d_Nz%d", grid.Nx, grid.Nz)

Lz = grid.Lz
width = 50kilometers
bump(x, y) = - Lz * (1 - 0.5 * exp(-(x^2 + y^2) / 2width^2))
grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bump))

#free_surface = ImplicitFreeSurface()
fft_preconditioner = FFTImplicitFreeSurfaceSolver(grid)
free_surface = ImplicitFreeSurface(solver_method=:PreconditionedConjugateGradient, preconditioner=fft_preconditioner)
#free_surface = ImplicitFreeSurface(solver_method=:HeptadiagonalIterativeSolver)

# Physics
Δx = grid.Lx / grid.Nx
κ₄h = Δx^4 / 1day
κz = 1e-2

diffusive_closure = VerticalScalarDiffusivity(VerticallyImplicitTimeDiscretization(), ν=κz, κ=κz)
horizontal_closure = HorizontalScalarBiharmonicDiffusivity(ν=κ₄h, κ=κ₄h)

model = HydrostaticFreeSurfaceModel(; grid, free_surface,
                                    coriolis = BetaPlane(latitude = -45),
                                    buoyancy = BuoyancyTracer(),
                                    closure = (diffusive_closure, horizontal_closure),
                                    tracers = (:b, :c),
                                    momentum_advection = WENO5(),
                                    tracer_advection = WENO5())

# Initial condition: a baroclinically unstable situation!
ramp(y, δy) = min(max(0, y/δy + 1/2), 1)

# Parameters
N² = 4e-6 # [s⁻²] buoyancy frequency / stratification
M² = 8e-8 # [s⁻²] horizontal buoyancy gradient

δy = 50kilometers
δz = 400
Lz = grid.Lz

δc = 2δy
δb = δy * M²
ϵb = 1e-2 * δb # noise amplitude

bᵢ(x, y, z) = N² * z + δb * ramp(y, δy) + ϵb * randn()
cᵢ(x, y, z) = exp(-y^2 / 2δc^2) * exp(-(z + Lz/4)^2 / 2δz^2)

set!(model, b=bᵢ, c=cᵢ)

Δt = 10minutes
simulation = Simulation(model; Δt, stop_time=30days)

#=
for i = 1:10
    @time [time_step!(simulation) for j = 1:10]
    try
        @show simulation.model.free_surface.implicit_step_solver.preconditioned_conjugate_gradient_solver.iteration
    catch; end
end
=#

wizard = TimeStepWizard(cfl=0.2, max_change=1.1, max_Δt=simulation.Δt)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(1))

wall_clock = Ref(time_ns())

function print_progress(sim)

    elapsed = 1e-9 * (time_ns() - wall_clock[])

    msg = @sprintf("Iter: %d, time: %s, wall time: %s, max|w|: %6.3e, m s⁻¹, next Δt: %s\n",
                   iteration(sim), prettytime(sim), prettytime(elapsed),
                   maximum(abs, sim.model.velocities.w), prettytime(sim.Δt))

    wall_clock[] = time_ns()

    try
        solver_iterations = sim.model.free_surface.implicit_step_solver.preconditioned_conjugate_gradient_solver.iteration
        msg *= @sprintf("solver iterations: %d", solver_iterations)
    catch
    end

    @info msg

    return nothing
end

simulation.callbacks[:print_progress] = Callback(print_progress, IterationInterval(10))

b, c = model.tracers
u, v, w = model.velocities
ζ = Field(∂x(v) - ∂y(u))

simulation.output_writers[:surface] = JLD2OutputWriter(model, (; ζ, b, c),
                                                       schedule = TimeInterval(1hour),
                                                       indices = (:, :, grid.Nz),
                                                       prefix = name * "_slices",
                                                       force = true)

simulation.output_writers[:fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers),
                                                      schedule = TimeInterval(10days),
                                                      with_halos = false,
                                                      prefix = name * "_fields",
                                                      force = true)

run!(simulation)

using GLMakie

set_theme!(Theme(font="Arial"))

filepath = name * "_slices.jld2"

ζt = FieldTimeSeries(filepath, "ζ")
bt = FieldTimeSeries(filepath, "b")
ct = FieldTimeSeries(filepath, "c")

t = ζt.times
Nt = length(t)

fig = Figure(resolution=(1800, 600))

axζ = Axis(fig[1, 1], xlabel="x (km)", ylabel="y (km)") #, title="Vorticity at surface")
axb = Axis(fig[1, 2], xlabel="x (km)", ylabel="y (km)") #, title="Buoyancy at surface")
axc = Axis(fig[1, 3], xlabel="x (km)", ylabel="y (km)") #, title="Tracer at surface")

slider = Slider(fig[3, :], range=1:Nt, startvalue=Nt)
n = slider.value

xζ, yζ, zζ = 1e-3 .* nodes((Face, Face, Center), grid)
xc, yc, zc = 1e-3 .* nodes((Center, Center, Center), grid)

ζⁿ = @lift interior(ζt[$n], :, :, 1)
bⁿ = @lift interior(bt[$n], :, :, 1)
cⁿ = @lift interior(ct[$n], :, :, 1)

ζlim = maximum(abs, parent(ζt[Nt]))
bmax = maximum(parent(bt[1]))
bmin = minimum(parent(bt[1]))

hmζ = heatmap!(axζ, xζ, yζ, ζⁿ, colorrange=(-ζlim, ζlim), colormap=:redblue)
hmb = heatmap!(axb, xc, yc, bⁿ, colorrange=(bmin, bmax), colormap=:thermal)
hmc = heatmap!(axc, xc, yc, cⁿ, colorrange=(0, 0.2), colormap=:deep)

Colorbar(fig[2, 1], hmζ, vertical=false, flipaxis=true, label="Vertical vorticity (s⁻¹)")
Colorbar(fig[2, 2], hmb, vertical=false, flipaxis=true, label="Buoyancy (m s⁻²)")
Colorbar(fig[2, 3], hmc, vertical=false, flipaxis=true, label="Tracer concentration")

title = @lift "Baroclinic adjustment at t = " * prettytime(t[$n])
Label(fig[0, :], title)

display(fig)

record(fig, "bumpy_baroclinic_adjustment.mp4", 1:Nt, framerate=12) do nn
    @info "Rendering frame $nn of $Nt..."
    n[] = nn
end


text:
# Translating
# https://proplot.readthedocs.io/en/latest/colorbars_legends.html#Axes-colorbars-and-legends
# to Julia

# Import proplot using PyCall
using PyCall
proplot = pyimport("proplot")

f, axs = proplot.subplots(nrows=2, share=0, axwidth="55mm", panelpad="1em")
[ax.format(suptitle="Stacked colorbars demo") for ax in axs]
N = 10
# Repeat for both axes
for (j, ax) in enumerate(axs)
    ax.format(
        xlabel="data", xlocator=range(0, 0.8, length=5),
        title="Subplot #$(j)"
    )
    for (i, (x0, y0, x1, y1, cmap, scale)) in enumerate((
         (0, 0.5, 1, 1, "grays", 0.5),
         (0, 0, 0.5, 0.5, "reds", 1),
         (0.5, 0, 1, 0.5, "blues", 2)
    ))
        j == 2 && i == 1 && continue
        data = rand(N, N) * scale
        x, y = range(x0, x1, length=N+1), range(y0, y1, length=N+1)
        m = ax.pcolormesh(
            x, y, data, cmap=cmap,
            levels=range(0, scale, length=11)
        )
        ax.colorbar(m, loc="l", label="dataset #$(i)")
    end
end

f.savefig("figs/Axes-colorbars-and-legends-2.svg", transparent=false)

max_stars_repo_name: natgeo-wong/JuliaProPlotExamples

max_issues_repo_name: natgeo-wong/JuliaProPlotExamples

max_forks_repo_name: natgeo-wong/JuliaProPlotExamples