text:
using ApproxFunFourier, ApproxFunBase, Test, SpecialFunctions, LinearAlgebra
import ApproxFunBase: testspace, testtransforms, testmultiplication, testraggedbelowoperator,
                    testbandedoperator, testblockbandedoperator, testbandedblockbandedoperator,
                    testcalculus, Block, Vec, testfunctional, EmptyDomain, UnionDomain
import SpecialFunctions: factorial

@testset "Periodic Domains" begin
    @test 0.1 ∈ PeriodicSegment(2π,0)
    @test 100.0 ∈ PeriodicSegment(0,2π)
    @test -100.0 ∈ PeriodicSegment(0,2π)

    @test 10.0 ∈ PeriodicLine()
    @test -10.0 ∈ PeriodicLine()
    @test -10.0+im ∉ PeriodicLine()

    @test ApproxFunBase.Vec(0,0.5) ∈ PeriodicSegment(ApproxFunBase.Vec(0.0,0), ApproxFunBase.Vec(0,1))

    @test ApproxFunBase.Vec(1,0) ∈ Circle((0.,0.),1.)

    @test ∂(PeriodicSegment() × ChebyshevInterval()) isa UnionDomain
    @test ∂(ChebyshevInterval() × PeriodicSegment()) isa UnionDomain
    @test ∂(PeriodicSegment() × PeriodicSegment()) isa EmptyDomain
end

@testset "Cos/SinSpace" begin
    for d in (PeriodicSegment(0.1,0.5),Circle(1.0+im,2.0))
        testtransforms(CosSpace(d);minpoints=2)
        testtransforms(SinSpace(d))
    end
    @test sum(Fun(1,CosSpace())) ≈ 2π
    @test sum(Fun(SinSpace(),[1])) == 0

    f=Fun(t->cos(t)+cos(3t),CosSpace)
    @test f(0.1) ≈ cos(0.1)+cos(3*0.1)
    @test (f*f-Fun(t->(cos(t)+cos(3t))^2,CosSpace)).coefficients|>norm <100eps()

    ## Calculus
    f=Fun(t->cos(t),CosSpace)
    D=Derivative(space(f))
    @test (D*f)(.1) ≈ -sin(.1)
    @test f'(.1) ≈ -sin(.1)

    f=Fun(t->sin(t),SinSpace)
    D=Derivative(space(f))
    @test (D*f)(.1) ≈ cos(.1)
    @test f'(.1) ≈ cos(.1)

    ## Multiplication
    s=Fun(t->(sin(t)+sin(2t))*cos(sin(t)),SinSpace)
    b=Fun(t->(sin(t)+sin(3t)),SinSpace)

    @test (s*s)(.1) ≈ s(.1)^2
    @test (s*b)(.1) ≈ s(.1)*b(.1)

    s=Fun(t->(cos(t)+cos(2t))*cos(cos(t)),CosSpace)
    b=Fun(t->(1+cos(t)+cos(3t)),CosSpace)

    @test (s*s)(.1) ≈ s(.1)^2
    @test (s*b)(.1) ≈ s(.1)*b(.1)

    s=Fun(t->(cos(t)+cos(2t))*cos(cos(t)),CosSpace)
    b=Fun(t->(sin(t)+sin(3t)),SinSpace)

    @test (s*b)(.1) ≈ s(.1)*b(.1)

    s=Fun(t->(sin(t)+sin(2t))*cos(sin(t)),SinSpace)
    b=Fun(t->(1+cos(t)+cos(3t)),CosSpace)

    @test (s*b)(.1) ≈ s(.1)*b(.1)

    ##  Norms
    @test sum(Fun(CosSpace(),[1.]))/(2π) ≈ 1.
    @test sum(Fun(CosSpace(),[0.,1.])^2)/(2π) ≈ 0.5
    @test sum(Fun(CosSpace(),[0.,0.,1.])^2)/(2π) ≈ 0.5
    @test sum(Fun(CosSpace(),[0.,0.,0.,1.])^2)/(2π) ≈ 0.5

    @test sum(Fun(SinSpace(),[0.,1.])^2)/(2π) ≈ 0.5
    @test sum(Fun(SinSpace(),[0.,0.,1.])^2)/(2π) ≈ 0.5
    @test sum(Fun(SinSpace(),[0.,0.,0.,1.])^2)/(2π) ≈ 0.5

    ## Bug in multiplicaiton
    @test Fun(SinSpace(),Float64[])^2 == Fun(SinSpace(),Float64[])
end

@testset "Taylor/Hardy" begin
    for d in (PeriodicSegment(0.1,0.5),Circle(1.0+im,2.0))
        testtransforms(Taylor(d))
        testtransforms(Hardy{false}(d))
    end
    f=Fun(exp,Taylor(Circle()))
    @test f(exp(0.1im)) ≈ exp(exp(0.1im))
    @test f(1.0) ≈ exp(1.0)
    g=Fun(z->1/(z-0.1),Hardy{false}(Circle()))
    @test (f(1.)+g(1.)) ≈ (exp(1.) + 1/(1-.1))

    @test Fun(Taylor())  == Fun(Taylor(),[0.,1.])

    @test Fun(Taylor())(1.0) ≈ 1.0
    @test Fun(Taylor(Circle(0.1,2.2)))(1.0) ≈ 1.0
    @test Fun(Taylor(Circle(0.1+0.1im,2.2)))(1.0) ≈ 1.0

    @test Multiplication(Fun(Taylor()),Taylor())[1:3,1:3] == [0. 0. 0.; 1. 0. 0.; 0. 1. 0.]

    # check's Derivative constructor works
    D=Derivative(Taylor(PeriodicSegment()))
end

@testset "Fourier" begin
    for d in (PeriodicSegment(0.1,0.5),Circle(1.0+im,2.0))
        testspace(Laurent(d);hasintegral=false)
        testspace(Fourier(d);hasintegral=false)
    end

    P = ApproxFunBase.plan_transform(Fourier(),4)
    v = randn(4) .+ im.*randn(4)
    @test P*v == P*real(v) + im*(P*imag(v))

    @test norm(Fun(cos,Fourier)'+Fun(sin,Fourier)) < 100eps()

    f=Fun(x->exp(-10sin((x-.1)/2)^2),Fourier)
    @test real(f)(.1) ≈ f(.1)

    f=Fun(cos,Fourier)
    @test norm((Derivative(space(f))^2)*f+f)<100eps()

    a=Fun(t->exp(cos(t)+sin(t)),Fourier)
    b=Fun(t->sin(t)+cos(3t)+1,Fourier)

    @test (a*b)(.1) ≈ a(.1)*b(.1)

    a=Fun(t->exp(cos(t)),CosSpace)
    b=Fun(t->sin(t)+cos(3t)+1,Fourier)

    @test (a*b)(.1) ≈ a(.1)*b(.1)

    a=Fun(t->sin(sin(t)),SinSpace)
    b=Fun(t->sin(t)+cos(3t)+1,Fourier)

    @test (a*b)(.1) ≈ a(.1)*b(.1)

    @test Fun(Fourier(),[1.])^2 ≈ Fun(Fourier(),[1.])

    ## Conversion between reverse
    C = Conversion(SinSpace()⊕CosSpace(),Fourier())
    @test C[Block(1), Block(1)] ≈ [0 1; 1 0]
    @test ApproxFunBase.defaultgetindex(C, Block.(1:2), Block.(1:2)) isa AbstractMatrix
    testbandedoperator(C)

    ## Test Multiplication
    mySin = Fun(Fourier(),[0,1.0])
    A = Multiplication(mySin,Fourier())
    @test A.op[1,1] == 0

    mySin = Fun(Fourier(),[0,1])
    A = Multiplication(mySin,Fourier())
    @test A.op[1,1] == 0

    @test norm(ApproxFunBase.Reverse(Fourier())*Fun(t->cos(cos(t-0.2)-0.1),Fourier()) - Fun(t->cos(cos(-t-0.2)-0.1),Fourier())) < 10eps()
    @test norm(ApproxFunBase.ReverseOrientation(Fourier())*Fun(t->cos(cos(t-0.2)-0.1),Fourier()) - Fun(t->cos(cos(t-0.2)-0.1),Fourier(PeriodicSegment(2π,0)))) < 10eps()
end

@testset "Laurent" begin
    f=Fun(x->exp(-10sin((x-.1)/2)^2),Laurent)
    @test f(.5) ≈ (Conversion(space(f),Fourier(domain(f)))*f)(.5)
    @test f(.5) ≈ Fun(f,Fourier)(.5)

    @test Fun(Laurent(0..2π),[1,1.,1.])(0.1) ≈ 1+2cos(0.1)
    @test Fun(Laurent(-1..1),[1,1.,1.])(0.1) ≈ 1+2cos(π*(0.1+1))
    @test Fun(Laurent(0..1),[1,1.,1.])(0.1) ≈ 1+2cos(2π*0.1)

    @test norm(Fun(cos,Laurent)'+Fun(sin,Laurent)) < 100eps()

    B=Evaluation(Laurent(0..2π),0,1)
    @test B*Fun(sin,domainspace(B)) ≈ 1.0

    ## Diagonal Derivative
    D = Derivative(Laurent())
    @test isdiag(D)
end

@testset "Circle" begin
    Γ=Circle(1.1,2.2)
    z=Fun(Fourier(Γ))
    @test space(z)==Fourier(Γ)
    @test z(1.1+2.2exp(0.1im)) ≈ 1.1+2.2exp(0.1im)

    @test abs(Fun(cos,Circle())(exp(0.1im))-cos(exp(0.1im)))<100eps()
    @test abs(Fun(cos,Circle())'(exp(0.1im))+sin(exp(0.1im)))<100eps()
    @test abs(Fun(cos,Circle())'(exp(0.1im))+Fun(sin,Circle())(exp(0.1im)))<100eps()

    @test norm(Fun(cos,Circle())'+Fun(sin,Circle()))<100eps()

    f=Fun(exp,Circle())
    @test component(f,1)(exp(0.1im)) ≈ exp(exp(0.1im))
    @test f(exp(0.1im)) ≈ exp(exp(0.1im))
    @test norm(f'-f)<100eps()
    @test norm(integrate(f)+1-f)<100eps()

    @test (Fun(z->sin(z)*cos(1/z),Circle())*Fun(z->exp(z)*airyai(1/z),Circle()))(exp(.1im)) ≈
                (z->sin(z)*cos(1/z)*exp(z)*airyai(1/z))(exp(.1im))

    for d in (Circle(),Circle(0.5),Circle(-0.1,2.))
        S=Taylor(d)
        D=Derivative(S)
        ef=Fun(exp,S)
        @test norm((D*ef-ef).coefficients)<4000eps()
        @test norm((D^2*ef-ef).coefficients)<200000eps()
        u=[Evaluation(S,0.),D-I]\[1.;0.]
        @test norm((u-ef).coefficients)<200eps()
        @test norm((Integral(S)*Fun(exp,S)+ef.coefficients[1]-ef).coefficients)<100eps()

        f=Fun(z->exp(1/z)-1,Hardy{false}(d))
        df=Fun(z->-1/z^2*exp(1/z),Hardy{false}(d))
        @test norm((Derivative()*f-df).coefficients)<1000eps()
        @test norm((Derivative()^2*f-df').coefficients)<100000eps()
        @test norm((f'-df).coefficients)<1000eps()
    end

    d=Circle()
    S=Taylor(d)
    D=Derivative(S)
    D-I
    ef=Fun(exp,S)
    @test norm((D*ef-ef).coefficients)<1000eps()
    @test norm((D^2*ef-ef).coefficients)<100000eps()
    u=[Evaluation(S,0.),D-I]\[1.;0.]

    # Check bug in off centre Circle
    c2=-0.1+.2im;r2=0.3;
    d2=Circle(c2,r2)
    z=Fun(identity,d2)

    @test z(-0.1+.2im+0.3*exp(0.1im)) ≈ (-0.1+.2im+0.3*exp(0.1im))

    # false Circle
    @test Fun(exp,Fourier(Circle(0.,1.,false)))(exp(0.1im)) ≈ exp(exp(.1im))
    @test Fun(exp,Laurent(Circle(0.,1.,false)))(exp(0.1im)) ≈ exp(exp(.1im))

    ## Reverse orientation
    f=Fun(z->1/z,Taylor(1/Circle()))
    @test f(exp(0.1im)) ≈ exp(-0.1im)

    ## exp(z)
    z=Fun(identity,Circle())
    cfs=exp(z).coefficients[1:2:end]
    for k=1:length(cfs)
        @test abs(cfs[k]-1/factorial(1.0(k-1))) ≤ 1E-10
    end

    ## Test bug in multiplication
    y = Fun(Circle())
    @test (y^2) ≈ Fun(z->z^2,domain(y))
end

@testset "Calculus" begin
    for f in (Fun(θ->sin(sin(θ)),SinSpace()),Fun(θ->cos(θ)+cos(3θ),CosSpace()),
                Fun(θ->sin(sin(θ)),Fourier()),Fun(θ->cos(θ)+cos(3θ),CosSpace()))
        @test norm(integrate(f)'-f)<10eps()
    end
end

@testset "Negatively oriented circles" begin
    f1 = Fun(z -> exp(1/z), Circle(0.0,0.2))
    f̃1 = Fun(z -> exp(1/z), Circle(0.0,0.2,false))
    f̃2 = Fun(z -> exp(1/z), Circle(0.0,0.3,false))

    @test f1(0.2exp(0.1im)) ≈ exp(1/(0.2exp(0.1im)))
    @test f̃1(0.2exp(0.1im)) ≈ exp(1/(0.2exp(0.1im)))
    @test f̃2(0.3exp(0.1im)) ≈ exp(1/(0.3exp(0.1im)))

    @test sum(f1) ≈ -sum(f̃1)
    @test sum(f̃1) ≈ sum(f̃2)
end

@testset "Fourier inplace" begin
    S = Fourier()

    x = [1.,2,3,4,5]
    y = similar(x)
    z = similar(x)
    P = ApproxFunBase.plan_transform(S, x)
    P! = ApproxFunBase.plan_transform!(S, x)
    mul!(y, P, x)
    @test x ≈ [1.,2,3,4,5]
    mul!(z, P!, x)
    @test x ≈ [1.,2,3,4,5]
    @test y ≈ z ≈ P*x ≈ P!*copy(x)

    P = ApproxFunBase.plan_itransform(S, x)
    P! = ApproxFunBase.plan_itransform!(S, x)
    mul!(y, P, x)
    @test x ≈ [1.,2,3,4,5]
    mul!(z, P!, x)
    @test x ≈ [1.,2,3,4,5]
    @test y ≈ z ≈ P*x ≈ P!*copy(x)
end

@testset "Vec circle" begin
    d=Circle((0.,0.),1.)
    f=Fun(xy->exp(-xy[1]-2cos(xy[2])),Fourier(d),40)
    @test f(cos(0.1),sin(0.1)) ≈ exp(-cos(0.1)-2cos(sin(0.1)))
    @test f(Vec(cos(0.1),sin(0.1))) ≈ exp(-cos(0.1)-2cos(sin(0.1)))

    f=Fun((x,y)->exp(-x-2cos(y)),Fourier(d),40)
    @test f(cos(0.1),sin(0.1)) ≈ exp(-cos(0.1)-2cos(sin(0.1)))

    f=Fun((x,y)->exp(-x-2cos(y)),Fourier(d))
    @test f(cos(0.1),sin(0.1)) ≈ exp(-cos(0.1)-2cos(sin(0.1)))
end

@testset "Operators" begin
    d=PeriodicSegment(0.,2π)
    a=Fun(t-> 1+sin(cos(10t)),d)
    D=Derivative(d)
    L=D+a

    @time testbandedoperator(D)
    @time testbandedoperator(Multiplication(a,Space(d)))

    f=Fun(t->exp(sin(t)),d)
    u=L\f

    @test norm(L*u-f) < 100eps()

    d=PeriodicSegment(0.,2π)
    a1=Fun(t->sin(cos(t/2)^2),d)
    a0=Fun(t->cos(12sin(t)),d)
    D=Derivative(d)
    L=D^2+a1*D+a0

    @time testbandedoperator(L)

    f=Fun(space(a1),[1,2,3,4,5])

    testbandedoperator(Multiplication(a0,Fourier(0..2π)))

    @test (Multiplication(a0,Fourier(0..2π))*f)(0.1)  ≈ (a0(0.1)*f(0.1))
    @test ((Multiplication(a1,Fourier(0..2π))*D)*f)(0.1)  ≈ (a1(0.1)*f'(0.1))
    @test (L.ops[1]*f)(0.1) ≈ f''(0.1)
    @test (L.ops[2]*f)(0.1) ≈ a1(0.1)*f'(0.1)
    @test (L.ops[3]*f)(0.1) ≈ a0(0.1)*f(0.1)
    @test (L*f)(0.1) ≈ f''(0.1)+a1(0.1)*f'(0.1)+a0(0.1)*f(0.1)

    f=Fun(t->exp(cos(2t)),d)
    u=L\f

    @test norm(L*u-f) < 1000eps()

    @time for M in (Multiplication(Fun(CosSpace(),[1.]),CosSpace()),
            Multiplication(Fun(CosSpace(),[1.]),SinSpace()),
            Multiplication(Fun(SinSpace(),[1.]),SinSpace()),
            Multiplication(Fun(SinSpace(),[1.]),CosSpace()),
            Derivative(SinSpace()),Derivative(CosSpace()))
        testbandedoperator(M)
    end
end

@testset "Integral equations" begin
    @time for S in (Fourier(Circle()),Laurent(Circle()),Taylor(Circle()),CosSpace(Circle()))
        testfunctional(DefiniteLineIntegral(S))
    end

    Σ = DefiniteIntegral()
    f1 = Fun(t->cos(cos(t))/t,Laurent(Circle()))
    f2 = Fun(t->cos(cos(t))/t,Fourier(Circle()))
    @test Σ*f1 ≈ Σ*f2

    f1=Fun(t->cos(cos(t)),Laurent(-π..π))
    f2=Fun(t->cos(cos(t)),Fourier(-π..π))
    @test Σ*f1 ≈ Σ*f2
end

@testset "tensor of mult for Fourier #507" begin
    mySin = Fun(Fourier(),[0,1])
    A = Multiplication(mySin,Fourier())
    L = A ⊗ A
    @test L[1,1] == 0
end

@testset "Piecewise + Constant" begin
    Γ=Circle() ∪ Circle(0.0,0.4)
    o=ones(Γ)
    @test o(1.) ≈ 1.0
    @test o(0.4) ≈ 1.0

    @time G=Fun(z->in(z,component(Γ,2)) ? [1 0; -1/z 1] : [z 0; 0 1/z],Γ)
    @test (G-I)(exp(0.1im)) ≈ (G(exp(0.1im))-I)
end

@testset "Array" begin
    @testset "Fourier Vector" begin
        a = [1 2; 3 4]
        f = Fun(θ->[sin(θ),sin(2θ)],Fourier())
        @test (a*f)(0.1) ≈ a*f(0.1)
        @test Fun(a)*f ≈ a*f
        @test Fun(a*Array(f)) ≈ a*f
        @test norm(f) ≈ sqrt(2π)
        @test norm(f,2) ≈ sqrt(2π)
    end

    @testset "CosSpace Vector" begin
        a = [1 2; 3 4]
        f = Fun(θ->[1,cos(θ)],CosSpace())
        @test (a*f)(0.1) ≈ [1+2cos(0.1); 3+4cos(0.1)]
        @test (a*f)(0.1) ≈ a*f(0.1)
        @test Fun(a)*f ≈ a*f
        @test Fun(a*Array(f)) ≈ a*f
    end

    @testset "CosSpace Matrix" begin
        a = [1 2; 3 4]
        m = Fun(θ->[1 cos(θ); cos(2θ) cos(cos(θ))],CosSpace())
        @test (a*m)(0.1) ≈ a*m(0.1)
        @test (m*a)(0.1) ≈ m(0.1)*a
        @test Fun(a)*m   ≈ a*m
        @test Fun(a*Array(m))   ≈ a*m

        @test (a+m)(0.1) ≈ a+m(0.1)
        @test (m+a)(0.1) ≈ m(0.1)+a

        @test (m+I)(0.1) ≈ m(0.1)+I
    end

    @testset "SinSpace Vector" begin
        a = [1 2; 3 4]
        f = Fun(θ->[sin(θ),sin(2θ)],SinSpace())
        @test (a*f)(0.1) ≈ a*f(0.1)
        @test Fun(a)*f ≈ a*f
        @test Fun(a*Array(f)) ≈ a*f

        @test all(sp -> sp isa SinSpace, space(a*f).spaces)
    end

    @testset "CosSpace Matrix" begin
        a = [1 2; 3 4]
        m = Fun(θ->[sin(3θ) sin(θ); sin(2θ) sin(sin(θ))],SinSpace())
        @test (a*m)(0.1) ≈ a*m(0.1)
        @test (m*a)(0.1) ≈ m(0.1)*a
        @test Fun(a)*m   ≈ a*m
        @test Fun(a*Array(m))   ≈ a*m

        @test all(sp -> sp isa SinSpace, space(a*m).spaces)

        @test (a+m)(0.1) ≈ a+m(0.1)
        @test (m+a)(0.1) ≈ m(0.1)+a

        @test (m+I)(0.1) ≈ m(0.1)+I
    end

    @testset "Two circles" begin
        Γ = Circle() ∪ Circle(0.5)

        f = Fun(z -> in(z,component(Γ,2)) ? 1 : z,Γ)
        @test f(exp(0.1im)) ≈ exp(0.1im)
        @test f(0.5exp(0.1im)) ≈ 1

        G = Fun(z -> in(z,component(Γ,2)) ? [1 -z^(-1); 0 1] :
                                            [z 0; 0 z^(-1)], Γ);

        @test G(exp(0.1im)) ≈ [exp(0.1im) 0 ; 0 exp(-0.1im)]
        @test G(0.5exp(0.1im)) ≈ [1 -2exp(-0.1im) ; 0 1]

        G1=Fun(Array(G)[:,1])

        @test G1(exp(0.1im)) ≈ [exp(0.1im),0.]
        @test G1(0.5exp(0.1im)) ≈ [1,0.]

        let f = G[1,1], sp = space(G1)[1,1] # test diagonal eltype bug
            D =  Diagonal([map(Multiplication,components(f),sp.spaces)...])
            @test eltype(Matrix(D)) == Operator{ComplexF64}
            testraggedbelowoperator(Multiplication(f, sp))
        end

        M = Multiplication(G, space(G1))
        testblockbandedoperator(M)

        for z in (0.5exp(0.1im),exp(0.2im))
            @test G[1,1](z) ≈ G[1](z)
            @test (M.op.ops[1,1]*G1[1])(z) ≈ M.f[1,1](z)*G1[1](z)
            @test (M.op.ops[2,1]*G1[1])(z) ≈ M.f[2,1](z)*G1[1](z)
            @test (M.op.ops[1,2]*G1[2])(z) ≈ M.f[1,2](z)*G1[2](z)
            @test (M.op.ops[2,2]*G1[2])(z) ≈ M.f[2,2](z)*G1[2](z)
        end

        u = M*G1
        @test norm(u(exp(.1im))-[exp(.2im),0])<100eps()
        @test norm(u(.5exp(.1im))-[1,0])<100eps()
    end

    @testset "Circle" begin
        G = Fun(z->[-1 -3; -3 -1]/z +
                   [ 2  2;  1 -3] +
                   [ 2 -1;  1  2]*z, Circle())

        @test G[1,1](exp(0.1im)) == G(exp(0.1im))[1,1]

        F̃ = Array((G-I)[:,1])
        F = (G-I)[:,1]

        @test Fun(F) ≡ F

        @test F(exp(0.1im)) ≈ [-exp(-0.1im)+1+2exp(0.1im);-3exp(-0.1im)+1+1exp(0.1im)]
        @test Fun(F̃,space(F))(exp(0.1im)) ≈ [-exp(-0.1im)+1+2exp(0.1im);-3exp(-0.1im)+1+1exp(0.1im)]

        @test coefficients(F̃,space(F)) == F.coefficients
        @test Fun(F̃,space(F)) == F

        @test F == Fun(vec(F),space(F))

        @test inv(G(exp(0.1im))) ≈ inv(G)(exp(0.1im))

        @test Fun(Matrix(I,2,2),space(G))(exp(0.1im)) ≈ Matrix(I,2,2)
        @test Fun(I,space(G))(exp(0.1im)) ≈ Matrix(I,2,2)
    end
end

@testset "Taylor()^2, checks bug in type of plan_transform" begin
    f = Fun((x,y)->exp((x-0.1)*cos(y-0.2)),Taylor()^2)
    @test f(0.2,0.3) ≈ exp(0.1*cos(0.1))
end

@testset "Periodic Poisson" begin
    d=PeriodicSegment()^2
    S=Space(d)

    f=Fun((x,y)->exp(-10(sin(x/2)^2+sin(y/2)^2)),d)
    A=Laplacian(d)+0.1I
    testbandedblockbandedoperator(A)
    @time u=A\f
    @test u(.1,.2) ≈ u(.2,.1)
    @test (lap(u)+.1u-f)|>coefficients|>norm < 1000000eps()
end

@testset "Low Rank" begin
    ## Periodic
    f=LowRankFun((x,y)->cos(x)*sin(y),PeriodicSegment(),PeriodicSegment())
    @test f(.1,.2) ≈ cos(.1)*sin(.2)

    f=LowRankFun((x,y)->cos(cos(x)+sin(y)),PeriodicSegment(),PeriodicSegment())
    @test f(.1,.2) ≈ cos(cos(.1)+sin(.2))
    @test norm(Float64[cos(cos(x)+sin(y)) for x=ApproxFunBase.vecpoints(f,1),y=ApproxFunBase.vecpoints(f,2)]-values(f))<10000eps()

    f=ProductFun((x,y)->cos(cos(x)+sin(y)),PeriodicSegment()^2)
    @test f(.1,.2) ≈ cos(cos(.1)+sin(.2))
    x,y=points(f)
    @test norm(Float64[cos(cos(x[k,j])+sin(y[k,j])) for k=1:size(f,1),j=1:size(f,2)]-values(f))<10000eps()

    d=PeriodicSegment()^2
    f=ProductFun((x,y)->exp(-10(sin(x/2)^2+sin(y/2)^2)),d)
    @test (transpose(f)-f|>coefficients|>norm)< 1000eps()
end

@testset "off domain evaluate" begin
    g = Fun(1, PeriodicSegment(Vec(0,-1) , Vec(π,-1)))
    @test g(0.1,-1) ≈ 1
    @test g(0.1,1) ≈ 0
end

@testset "PeriodicLine" begin
    d=PeriodicLine()
    D=Derivative(d)

    f = Fun(x->sech(x-0.1),d,200)
    @test f(1.) ≈ sech(1-0.1)

    f=Fun(x->sech(x-0.1),d)
    @test f(1.) ≈ sech(1-0.1)

    @test ≈((D*f)(.2),-0.0991717226583897;atol=100000eps())
    @test ≈((D^2*f)(.2),-0.9752522555114987;atol=1000000eps())

    f=Fun(z->2exp(z^2),PeriodicLine(0.,π/2))
    @test f(1.1im) ≈ 2exp(-1.1^2)
end


text:
# add more julia processes
using Distributed
nprocs() < 4 ? addprocs(5-nprocs()) : nothing

# tell Julia that you want to use these modules/namespaces
using RoME
## Inter-operating visualization packages for Caesar/RoME/IncrementalInference exist
using RoMEPlotting

# start with an empty factor graph object
fg = initfg()

# Add the first pose :x0
addVariable!(fg, :x0, Pose2)

# Add at a fixed location PriorPose2 to pin :x0 to a starting location (10,10, pi/4)
addFactor!(fg, [:x0], PriorPose2( MvNormal([10; 10; pi/6.0], Matrix(Diagonal([0.1;0.1;0.05].^2))) ) )

# Drive around in a hexagon
for i in 0:5
    psym = Symbol("x$i")
    nsym = Symbol("x$(i+1)")
    addVariable!(fg, nsym, Pose2)
    pp = Pose2Pose2(MvNormal([10.0;0;pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
    addFactor!(fg, [psym;nsym], pp )
end

# perform inference, and remember first runs are slower owing to Julia's just-in-time compiling
getSolverParams(fg).drawtree=true
getSolverParams(fg).showtree=true
tree, smt, hist = solveTree!(fg)

pl = drawPoses(fg) #, meanmax=:mean)
# For scripting use-cases you can export the image
pl |> Gadfly.PDF("/tmp/test1.pdf", 20cm, 10cm)  # or PNG(...)

# another type of plot
plotPose(fg, :x6)

# Add landmarks with Bearing range measurements
addVariable!(fg, :l1, Point2, tags=[:LANDMARK])
p2br = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))
addFactor!(fg, [:x0; :l1], p2br )

# Initialize :l1 numerical values but do not rerun solver
ensureAllInitialized!(fg)
pl = drawPosesLandms(fg)
Gadfly.draw(Gadfly.PDF("/tmp/test2.pdf", 20cm, 10cm),pl)  # or PNG(...)

# Add landmarks with Bearing range measurements
p2br2 = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))
addFactor!(fg, [:x6; :l1], p2br2 )

# solve
batchSolve!(fg, drawpdf=true)

# redraw
pl = drawPosesLandms(fg, meanmax=:mean)
Gadfly.draw(Gadfly.PDF("/tmp/test3.pdf", 20cm, 10cm),pl)  # or PNG(...)

#

ext: jl


text:
using GalerkinSparseGrids
using Base.Test
using Cubature

#---------------------------------------
# Testing regular hier DG reconstruction
#---------------------------------------

print("Testing full discontinuous Galerkin (DG) reconstruction 1-D... ")

for k in 1:5
    for l in 1:6
        dict = coeffs_DG(1, k, l, x->sin(4*x[1]); scheme="full")
        err = x->(reconstruct_DG(dict, [x])-sin(4*x))^2
        @test hquadrature(err, 0, 1; abstol=1.0e-10)[1] < 1/(1<<(l+k-1))
    end
end

for k in 1:5
    for l in 1:6
        dict = coeffs_DG(1, k, l, x->sin(4*x[1]); scheme="sparse")
        err = x->(reconstruct_DG(dict, [x])-sin(4*x))^2
        @test hquadrature(err, 0, 1; abstol=1.0e-10)[1] < 1/(1<<(l+k-1))
    end
end

println("Test Passed.")

# Sparsification should be trivial in 1D

print("Testing that both reconstructions are equivalent... ")

for k in 1:5
    for l in 1:6
        full_coeffs = coeffs_DG(1, k, l, x->sin(4*x[1]); scheme="full")
        sparse_coeffs = coeffs_DG(1, k, l, x->sin(4*x[1]); scheme="full")
        diff_full = x->(reconstruct_DG(full_coeffs, [x])-sin(4*x))^2
        diff_sparse = x->(reconstruct_DG(sparse_coeffs, [x])-sin(4*x))^2
        diff_both = x->(reconstruct_DG(full_coeffs, [x])-reconstruct_DG(sparse_coeffs, [x]))^2

        @test abs(hquadrature(diff_full,0,1;abstol=1.0e-9)[1]-
                  hquadrature(diff_sparse,0,1;abstol=1.0e-9)[1]) < 1.0e-14
        @test hquadrature(diff_both,0,1;abstol=1.0e-9)[1] < 1.0e-14
    end
end

println("Test Passed.")

# Testing sparse DG reconstruction in multidimensional space:

print("Testing sparse DG basis reconstruction 2-D... ")

for k in 1:5
    for l in 1:6
        dict = coeffs_DG(2, k, l, x->sin(4*x[1]+x[2]))
        err = x->(reconstruct_DG(dict, [x[1],x[2]])-sin(4*x[1]+x[2]))^2
        @test hcubature(err, [0,0], [1,1]; abstol=1.0e-10, maxevals=500)[1] < 1/(1<<(l+k-2))
    end
end

println("Test Passed.")


text:
# =============================================================================
# Draw christmas tree
# Luis Lopez Diaz - llopeza@uninorte.edu.co
# =============================================================================
# Input data
l = 28;                         # Number of lines for drawing
# =============================================================================
# cheking data
if l < 0
    error("Number of line needs to be positive.");
    return;
end
if ~(typeof(l) == Int)
    error("Number of line needs to be integer");
    return;
end
# =============================================================================
# main code
endw = "*";
form = "*o";
for i in 1:l
    for j in i+1:l print(" "); end  # Print blank space
    for j in 1:i-1 print(form); end # Print form
    print(endw * "\n")              # Print ending
end


text:
# install.jl collects scripts to install many third party libraries

export install_adept, install_blas, install_openmpi, install_hypre,
    install_had, install_mfem, install_matplotlib, install_eigen

function install_blas(blas_binary::Bool = true)
    if Sys.iswindows()
        if isfile(joinpath(ADCME.LIBDIR, "openblas.lib"))
            @info "openblas.lib exists"
            return
        end
        if blas_binary
            http_file("https://github.com/kailaix/tensorflow-1.15-include/releases/download/v0.1.0/openblas.lib", joinpath(ADCME.LIBDIR, "openblas.lib"))
            return
        end
        @info "You are building openblas from source on Windows, and this process may take a long time.
Alternatively, you can place your precompiled binary to $(joinpath(ADCME.LIBDIR, "openblas.lib"))"
        PWD = pwd()
        change_directory()
        http_file("https://github.com/xianyi/OpenBLAS/archive/v0.3.9.zip", "OpenBlas.zip")
        uncompress("OpenBLAS.zip", "OpenBlas-0.3.9")
        change_directory("OpenBlas-0.3.9/build")
        require_cmakecache() do
            ADCME.cmake(CMAKE_ARGS="-DCMAKE_Fortran_COMPILER=flang -DBUILD_WITHOUT_LAPACK=no -DNOFORTRAN=0 -DDYNAMIC_ARCH=ON")
        end
        require_library("lib/Release/openblas") do
            ADCME.make()
        end
        change_directory("lib/Release")
        copy_file("openblas.lib", joinpath(ADCME.LIBDIR, "openblas.lib"))
        cd(PWD)
    else
        required_file = get_library(joinpath(ADCME.LIBDIR, "openblas"))
        if !isfile(required_file)
            files = readdir(ADCME.LIBDIR)
            files = filter(x->!isnothing(x), match.(r"(libopenblas\S*.\S*)", files))[1]
            target = joinpath(ADCME.LIBDIR, files[1])
            symlink(target, required_file)
            @info "Symlink $(required_file) --> $(files[1])"
        end
    end
end

"""
    install_eigen()

Install eigen library
"""
function install_eigen()
    require_file("$(ADCME.PREFIXDIR)/eigen.zip") do
        download("https://gitlab.com/libeigen/eigen/-/archive/3.3.7/eigen-3.3.7.zip","$(ADCME.PREFIXDIR)/eigen.zip")
    end

    require_file("$(ADCME.PREFIXDIR)/eigen3") do
        UNZIP =  joinpath(ADCME.BINDIR, "unzip")
        if Sys.iswindows()
            require_file("$(ADCME.PREFIXDIR)/unzip.exe") do
                download("http://stahlworks.com/dev/unzip.exe", joinpath(ADCME.PREFIXDIR, "unzip.exe"))
            end
            UNZIP =  joinpath(ADCME.PREFIXDIR, "unzip.exe")
        end
        run(`$UNZIP -qq $(ADCME.PREFIXDIR)/eigen.zip -d $(ADCME.PREFIXDIR)`)
        mv("$(ADCME.PREFIXDIR)/eigen-3.3.7", "$(ADCME.PREFIXDIR)/eigen3", force=true)
    end

end

"""
    install_adept()

Install adept-2 library: https://github.com/rjhogan/Adept-2
"""
function install_adept(; blas_binary::Bool = true)
    change_directory()
    git_repository("https://github.com/ADCMEMarket/Adept-2", "Adept-2")
    cd("Adept-2/adept")
    install_blas(blas_binary)
    cp("$(@__DIR__)/../deps/AdeptCMakeLists.txt", "./CMakeLists.txt", force=true)
    change_directory("build")
    require_cmakecache() do
        ADCME.cmake()
    end
    if Sys.iswindows()
        require_file("$(ADCME.LIBDIR)/adept.lib") do
            ADCME.make()
        end
    else
        require_library("$(ADCME.LIBDIR)/adept") do
            ADCME.make()
        end
    end
end

function install_openmpi()
    filepath = joinpath(@__DIR__, "..", "deps", "install_openmpi.jl")
    include(filepath)
end

function install_hypre()
    filepath = joinpath(@__DIR__, "..", "deps", "install_hypre.jl")
    include(filepath)
end

function install_mfem()
    PWD = pwd()
    change_directory()
    http_file("https://github.com/kailaix/mfem/archive/shared-msvc-dev.zip", "mfem.zip")
    uncompress("mfem.zip", "mfem-shared-msvc-dev")
    change_directory("mfem-shared-msvc-dev/build")
    require_file("CMakeCache.txt") do
        ADCME.cmake(CMAKE_ARGS = ["-DCMAKE_INSTALL_PREFIX=$(joinpath(ADCME.LIBDIR, ".."))", "-SHARED=YES", "-STATIC=NO"])
    end
    require_library("mfem") do
        ADCME.make()
    end
    require_file(joinpath(ADCME.LIBDIR, get_library_name("mfem"))) do
        if Sys.iswindows()
            run_with_env(`cmd /c $(ADCME.CMAKE) --install .`)
        else
            run_with_env(`$(ADCME.NINJA) install`)
        end
    end
    if Sys.iswindows()
mfem_h = """
// Auto-generated file.
#undef NO_ERROR
#undef READ_ERROR
#undef WRITE_ERROR
#undef ALIAS
#undef REGISTERED
#include "mfem/mfem.hpp"
"""
        open(joinpath(ADCME.LIBDIR, "..",  "include", "mfem.hpp"), "w") do io
            write(io, mfem_h)
            @info "Fixed mfem.hpp definitions"
        end
    end
    cd(PWD)
end

function install_had()
    change_directory()
    git_repository("https://github.com/kailaix/had", "had")
end

function install_matplotlib()
    PIP = get_pip()
    run(`$PIP install matplotlib`)
    if Sys.isapple()
        CONDA = get_conda()
        pkgs = run(`$CONDA list`)
        if occursin("pyqt", pkgs)
            return
        end
        if !isdefined(Main, :Pkg)
            error("Package Pkg must be imported in the main module using `import Pkg` or `using Pkg`")
        end
        run(`$CONDA install -y pyqt`)
        Main.Pkg.build("PyPlot")
    end
end


text:
println("Hello world !")

# something on main branch


text:
include("app/devtools.jl")
include("app/supporttypes.jl")
include("app/config.jl")
include("app/dashapp.jl")
include("app/callbacks.jl")


text:
module ProtoBufTestCodec
using ProtoBuf
using Test
using Random
import ProtoBuf.meta

macro _rand_int(T,mx,a)
   esc(quote
       $(T)(round(rand() * $(mx)) + $(a))
   end)
end

const TestTypeJType = Ref{Type}(Int64)
const TestTypeWType = Ref{Symbol}(:int64)
const TestTypeFldNum = Ref{Int}(1)
const TestTypePack = Ref{Vector{Symbol}}(ProtoBuf.DEF_PACK)
mutable struct TestType <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function TestType(; kwargs...)
        obj = new(meta(TestType), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end
function meta(::Type{TestType})
    allflds = Pair{Symbol,Union{Type,String}}[:val => TestTypeJType[]]
    wtypes = Dict{Symbol,Symbol}(:val => TestTypeWType[])
    meta(ProtoMeta(TestType), TestType, allflds, ProtoBuf.DEF_REQ, [TestTypeFldNum[]], ProtoBuf.DEF_VAL, TestTypePack[], wtypes, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
end
function Base.getproperty(obj::TestType, name::Symbol)
    if name === :val
        return (obj.__protobuf_jl_internal_values[name])::Any
    else
        getfield(obj, name)
    end
end

mutable struct TestStr <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function TestStr(; kwargs...)
        obj = new(meta(TestStr), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end
function meta(::Type{TestStr})
    allflds = Pair{Symbol,Union{Type,String}}[:val => AbstractString]
    meta(ProtoMeta(TestStr), TestStr, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
end
function Base.getproperty(obj::TestStr, name::Symbol)
    if name === :val
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct TestOptional <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function TestOptional(; kwargs...)
        obj = new(meta(TestOptional), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end
const TestOptionalReq = Symbol[]
function meta(::Type{TestOptional})
    allflds = Pair{Symbol,Union{Type,String}}[:sVal1 => TestStr, :sVal2 => TestStr, :iVal2 => Array{Int64,1}]
    meta(ProtoMeta(TestOptional), TestOptional, allflds, TestOptionalReq, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
end
function Base.getproperty(obj::TestOptional, name::Symbol)
    if name === :sVal1
        return (obj.__protobuf_jl_internal_values[name])::TestStr
    elseif name === :sVal2
        return (obj.__protobuf_jl_internal_values[name])::TestStr
    elseif name === :iVal2
        return (obj.__protobuf_jl_internal_values[name])::Array{Int64,1}
    else
        getfield(obj, name)
    end
end

mutable struct TestNested <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function TestNested(; kwargs...)
        obj = new(meta(TestNested), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end

end
const TestNestedReq = Symbol[]
function meta(::Type{TestNested})
    allflds = Pair{Symbol,Union{Type,String}}[:fld1 => TestType, :fld2 => TestOptional, :fld3 => Array{TestStr,1}]
    meta(ProtoMeta(TestNested), TestNested, allflds, TestNestedReq, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
end
function Base.getproperty(obj::TestNested, name::Symbol)
    if name === :fld1
        return (obj.__protobuf_jl_internal_values[name])::TestType
    elseif name === :fld2
        return (obj.__protobuf_jl_internal_values[name])::TestOptional
    elseif name === :fld3
        return (obj.__protobuf_jl_internal_values[name])::Array{TestStr,1}
    else
        getfield(obj, name)
    end
end

mutable struct TestDefaults <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function TestDefaults(; kwargs...)
        obj = new(meta(TestDefaults), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end
function meta(::Type{TestDefaults})
    allflds = Pair{Symbol,Union{Type,String}}[:iVal1 => Int64, :sVal2 => AbstractString, :iVal2 => Array{Int64,1}]
    defaults = Dict{Symbol,Any}(:iVal1 => 10, :iVal2 => [1,2,3])
    meta(ProtoMeta(TestDefaults), TestDefaults, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, defaults, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
end
function Base.getproperty(obj::TestDefaults, name::Symbol)
    if name === :iVal1
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :sVal2
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :iVal2
        return (obj.__protobuf_jl_internal_values[name])::Array{Int64,1}
    else
        getfield(obj, name)
    end
end

mutable struct TestOneofs <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function TestOneofs(; kwargs...)
        obj = new(meta(TestOneofs), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end
function meta(::Type{TestOneofs})
    allflds = Pair{Symbol,Union{Type,String}}[:iVal1 => Int64, :iVal2 => Int64, :iVal3 => Int64]
    oneofs = Int[0,1,1]
    oneof_names = [:optval]
    meta(ProtoMeta(TestOneofs), TestOneofs, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
end
function Base.getproperty(obj::TestOneofs, name::Symbol)
    if name === :iVal1
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :iVal2
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :iVal3
        return (obj.__protobuf_jl_internal_values[name])::Int64
    else
        getfield(obj, name)
    end
end

mutable struct TestMaps <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function TestMaps(; kwargs...)
        obj = new(meta(TestMaps), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end
function meta(::Type{TestMaps})
    allflds = Pair{Symbol,Union{Type,String}}[:d1 => Dict{Int,Int}, :d2 => Dict{Int32,String}, :d3 => Dict{String,String}]
    meta(ProtoMeta(TestMaps), TestMaps, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
end
function Base.getproperty(obj::TestMaps, name::Symbol)
    if name === :d1
        return (obj.__protobuf_jl_internal_values[name])::Dict{Int,Int}
    elseif name === :d2
        return (obj.__protobuf_jl_internal_values[name])::Dict{Int32,String}
    elseif name === :d3
        return (obj.__protobuf_jl_internal_values[name])::Dict{String,String}
    else
        getfield(obj, name)
    end
end

mutable struct TestFilled <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function TestFilled(; kwargs...)
        obj = new(meta(TestFilled), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end
function meta(::Type{TestFilled})
    allflds = Pair{Symbol,Union{Type,String}}[:fld1 => TestType, :fld2 => TestType]
    meta(ProtoMeta(TestFilled), TestFilled, allflds, [:fld1], ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
end
function Base.getproperty(obj::TestFilled, name::Symbol)
    if name === :fld1
        return (obj.__protobuf_jl_internal_values[name])::TestType
    elseif name === :fld2
        return (obj.__protobuf_jl_internal_values[name])::TestType
    else
        getfield(obj, name)
    end
end

const TestEnum = (;[
    Symbol("UNIVERSAL") => Int32(0),
    Symbol("WEB") => Int32(1),
    Symbol("IMAGES") => Int32(2),
    Symbol("LOCAL") => Int32(3),
    Symbol("NEWS") => Int32(4),
    Symbol("PRODUCTS") => Int32(5),
    Symbol("VIDEO") => Int32(6),
]...)

assert_equal(::Type{Array{T,1}}, ::Type{Array{U,1}}) where {T,U} = @test (T <: U) || (U <: T)
assert_equal(T::Type, U::Type) = @test (T <: U) || (U <: T)
assert_equal(val1::Bool, val2::Bool) = @test val1 == val2
assert_equal(val1::AbstractString, val2::AbstractString) = @test val1 == val2
assert_equal(val1::Number, val2::Number) = @test val1 == val2
assert_equal(val1::Array{T,1}, val2::Array{T,1}) where {T<:Number} = @test val1 == val2
assert_equal(val1::Array{T,1}, val2::Array{T,1}) where {T<:AbstractString} = @test val1 == val2
assert_equal(val1::Dict{K,V}, val2::Dict{K,V}) where {K,V} = @test val1 == val2
function assert_equal(val1, val2)
    typ1 = typeof(val1)
    typ2 = typeof(val2)
    assert_equal(typ1, typ2)

    for fld in propertynames(typ1)
        fldfill1 = hasproperty(val1, fld)
        fldfill2 = hasproperty(val2, fld)
        @test fldfill1 == fldfill2
        fldfill1 && assert_equal(getproperty(val1, fld), getproperty(val2, fld))
    end
end

function test_types()
    pb = PipeBuffer()

    @testset "types" begin
        @testset "enum" begin
            TestTypeJType[] = Int32
            TestTypeWType[] = :int32
            TestTypePack[] = ProtoBuf.DEF_PACK
            TestTypeFldNum[] = @_rand_int(Int, 100, 1)
            testmeta = meta(TestType)
            testval = TestType(; val=@_rand_int(Int32, 10^9, 0))
            readval = TestType()
            writeproto(pb, testval, testmeta)
            readproto(pb, readval, testmeta)
            assert_equal(testval, readval)
        end

        @testset "integers" begin
            let typs = [Int32,Int64,UInt32,UInt64,Int32,Int64,UInt64,Int64,UInt32,Int32], ptyps=[:int32,:int64,:uint32,:uint64,:sint32,:sint64,:fixed64,:sfixed64,:fixed32,:sfixed32]
                for (typ,ptyp) in zip(typs,ptyps)
                    TestTypeJType[] = typ
                    TestTypeWType[] = ptyp
                    for idx in 1:100
                        TestTypeFldNum[] = @_rand_int(Int, 100, 1)
                        testmeta = meta(TestType)
                        testval = TestType(; val=convert(typ, @_rand_int(UInt32, 10^9, 0)))
                        readval = TestType()
                        writeproto(pb, testval, testmeta)
                        readproto(pb, readval, testmeta)
                        assert_equal(testval, readval)
                    end
                end
            end

            let typs = [Int32,Int64,Int32,Int64], ptyps=[:int32,:int64,:sint32,:sint64]
                for (typ,ptyp) in zip(typs,ptyps)
                    TestTypeJType[] = typ
                    TestTypeWType[] = ptyp
                    for idx in 1:100
                        TestTypeFldNum[] = convert(typ, -1 * @_rand_int(Int32, 10^9, 0))
                        testmeta = meta(TestType)
                        testval = TestType(; val=convert(typ, @_rand_int(UInt32, 10^9, 0)))
                        readval = TestType()
                        writeproto(pb, testval, testmeta)
                        readproto(pb, readval, testmeta)
                        assert_equal(testval, readval)
                    end
                end
            end
        end

        @testset "varint overflow" begin
            ProtoBuf._write_uleb(pb, Int64(-1))
            @test ProtoBuf._read_uleb(pb, Int8) == 0
            ProtoBuf._write_uleb(pb, Int64(1))
            @test ProtoBuf._read_uleb(pb, Int8) == 1
            write(pb, 0xff)
            ProtoBuf._write_uleb(pb, Int64(-1))
            @test ProtoBuf._read_uleb(pb, Int32) == 0
        end

        @testset "bool" begin
            let typs = [Bool], ptyps=[:bool]
                for (typ,ptyp) in zip(typs,ptyps)
                    TestTypeJType[] = typ
                    TestTypeWType[] = ptyp
                    for idx in 1:100
                        TestTypeFldNum[] = @_rand_int(Int, 100, 1)
                        testmeta = meta(TestType)
                        testval = TestType(; val=convert(typ, @_rand_int(UInt32, 1, 0)))
                        readval = TestType()
                        writeproto(pb, testval, testmeta)
                        readproto(pb, readval, testmeta)
                        assert_equal(testval, readval)
                    end
                end
            end
        end

        @testset "double, float" begin
            let typs = [Float64,Float32], ptyps=[:double,:float]
                for (typ,ptyp) in zip(typs,ptyps)
                    TestTypeJType[] = typ
                    TestTypeWType[] = ptyp
                    for idx in 1:100
                        TestTypeFldNum[] = @_rand_int(Int, 100, 1)
                        testmeta = meta(TestType)
                        testval = TestType(; val=convert(typ, @_rand_int(UInt32, 10^9, 0)))
                        readval = TestType()
                        writeproto(pb, testval, testmeta)
                        readproto(pb, readval, testmeta)
                        assert_equal(testval, readval)
                    end
                end
            end
        end

        @testset "string" begin
            let typs = [AbstractString], ptyps=[:string]
                for (typ,ptyp) in zip(typs,ptyps)
                    TestTypeJType[] = typ
                    TestTypeWType[] = ptyp
                    for idx in 1:100
                        TestTypeFldNum[] = @_rand_int(Int, 100, 1)
                        testmeta = meta(TestType)
                        testval = TestType(; val=randstring(50))
                        readval = TestType()
                        writeproto(pb, testval, testmeta)
                        readproto(pb, readval, testmeta)
                        assert_equal(testval, readval)
                    end
                end
            end
        end
    end
end

function test_repeats()
    pb = PipeBuffer()

    @testset "Repeated" begin
        @testset "Empty repeated int64" begin
            TestTypeJType[] = Vector{Int64}
            TestTypeWType[] = :int64
            TestTypePack[] = ProtoBuf.DEF_PACK
            testval = TestType()
            readval = TestType()
            testmeta = meta(TestType)
            writeproto(pb, testval, testmeta)
            readproto(pb, readval, testmeta)
            assert_equal(testval, readval)
            @test hasproperty(readval, :val)
        end

        @testset "Repeated int64" begin
            TestTypeJType[] = Vector{Int64}
            TestTypeWType[] = :int64
            TestTypePack[] = ProtoBuf.DEF_PACK
            for idx in 1:100
                TestTypeFldNum[] = @_rand_int(Int, 100, 1)
                testval = TestType(; val=collect(Int64, randstring(50)))
                readval = TestType()
                testmeta = meta(TestType)
                writeproto(pb, testval, testmeta)
                readproto(pb, readval, testmeta)
                assert_equal(testval, readval)
            end
        end

        @testset "Repeated and packed int64" begin
            TestTypePack[] = Symbol[:val]
            for idx in 1:100
                TestTypeFldNum[] = @_rand_int(Int, 100, 1)
                testval = TestType(; val=collect(Int64, randstring(50)))
                readval = TestType()
                testmeta = meta(TestType)
                writeproto(pb, testval, testmeta)
                readproto(pb, readval, testmeta)
                assert_equal(testval, readval)
            end
        end

        @testset "Repeated string" begin
            TestTypeJType[] = Vector{AbstractString}
            TestTypeWType[] = :string
            TestTypePack[] = ProtoBuf.DEF_PACK
            for idx in 1:100
                testval = TestType(; val=AbstractString[randstring(5) for i in 1:10])
                readval = TestType()
                testmeta = meta(TestType)
                writeproto(pb, testval, testmeta)
                readproto(pb, readval, testmeta)
                assert_equal(testval, readval)
            end
        end
    end
end

function test_optional()
    @testset "Optional fields" begin
        pb = PipeBuffer()
        testval = TestOptional(; sVal1=TestStr(; val=""), sVal2=TestStr(; val=""), iVal2=Int64[1,2,3])
        readval = TestOptional(; sVal1=TestStr(; val=""), sVal2=TestStr(; val=""), iVal2=Int64[])

        for idx in 1:100
            testval.sVal1 = TestStr(; val=string(@_rand_int(Int, 100, 0)))
            testval.sVal2 = TestStr(; val=randstring(5))
            testval.iVal2 = Int64[@_rand_int(Int,100,0) for i in 1:10]

            empty!(TestOptionalReq)
            kwargs = Dict{Symbol,Any}(:iVal2 => Int64[@_rand_int(Int,100,0) for i in 1:10])
            if rand(Bool)
                push!(TestOptionalReq, :sVal1)
                kwargs[:sVal1] = TestStr(; val=string(@_rand_int(Int, 100, 0)))
            end
            if rand(Bool)
                push!(TestOptionalReq, :sVal2)
                kwargs[:sVal2] = TestStr(; val=randstring(5))
            end
            testval = TestOptional(; kwargs...)
            readval = TestOptional()
            testmeta = meta(TestOptional)
            writeproto(pb, testval, testmeta)
            readproto(pb, readval, testmeta)
            assert_equal(testval, readval)
        end
    end
end

function test_nested()
    @testset "Nested types" begin
        pb = PipeBuffer()

        TestTypeJType[] = Int64
        TestTypeWType[] = :int64
        TestTypePack[] = ProtoBuf.DEF_PACK
        TestTypeFldNum[] = 1

        for idx in 1:100
            o1 = rand(Bool)
            o2 = rand(Bool)
            o21 = rand(Bool)
            o22 = rand(Bool)

            empty!(TestNestedReq)
            testnestedkwargs = Dict{Symbol,Any}()
            if o1
                push!(TestNestedReq, :fld1)
                testnestedkwargs[:fld1] = TestType(; val=@_rand_int(Int64, 10^9, 0))
            end
            if o2
                push!(TestNestedReq, :fld2)
                empty!(TestOptionalReq)
                testfld2kwargs = Dict{Symbol,Any}(:iVal2=>Int64[@_rand_int(Int, 100, 0) for i in 1:10])
                if o21
                    push!(TestOptionalReq, :sVal1)
                    testfld2kwargs[:sVal1] = TestStr(; val=string(@_rand_int(Int, 100, 0)))
                end
                if o22
                    push!(TestOptionalReq, :sVal2)
                    testfld2kwargs[:sVal2] = TestStr(; val=randstring(5))
                end
                testnestedkwargs[:fld2] = TestOptional(; testfld2kwargs...)
            end

            testval = TestNested(; testnestedkwargs...)
            readval = TestNested()
            testmeta = meta(TestNested)

            writeproto(pb, testval, testmeta)
            readproto(pb, readval, testmeta)

            assert_equal(testval, readval)
        end
    end
end

function test_defaults()
    @testset "Default values" begin
        pb = PipeBuffer()

        testval = TestDefaults()
        readval = TestDefaults()
        testval.iVal1 = @_rand_int(Int, 100, 0)
        writeproto(pb, testval)
        readproto(pb, readval)

        assert_equal(TestDefaults(; iVal1=testval.iVal1, sVal2="", iVal2=[1,2,3]), readval)
    end
end

function test_oneofs()
    @testset "Oneofs" begin
        testval = TestOneofs(; iVal1=1, iVal3=3)
        @test isfilled(testval)
        @test hasproperty(testval, :iVal1)
        @test !hasproperty(testval, :iVal2)
        @test hasproperty(testval, :iVal3)
        @test which_oneof(testval, :optval) === :iVal3

        testval.iVal2 = 10
        @test hasproperty(testval, :iVal1)
        @test hasproperty(testval, :iVal2)
        @test !hasproperty(testval, :iVal3)
        @test which_oneof(testval, :optval) === :iVal2

        testval.iVal1 = 10
        @test hasproperty(testval, :iVal1)
        @test hasproperty(testval, :iVal2)
        @test !hasproperty(testval, :iVal3)
        @test which_oneof(testval, :optval) === :iVal2

        testval.iVal3 = 10
        @test hasproperty(testval, :iVal1)
        @test !hasproperty(testval, :iVal2)
        @test hasproperty(testval, :iVal3)
        @test which_oneof(testval, :optval) === :iVal3
    end
end

function test_maps()
    @testset "Maps" begin
        pb = PipeBuffer()

        testval = TestMaps()
        readval = TestMaps()
        writeproto(pb, testval)
        readproto(pb, readval)
        assert_equal(testval, readval)

        testval = TestMaps()
        readval = TestMaps()
        testval.d1 = Dict{Int,Int}()
        writeproto(pb, testval)
        readproto(pb, readval)
        @test hasproperty(readval, :d1)
        assert_equal(testval.d1, readval.d1)

        testval = TestMaps()
        readval = TestMaps()
        testval.d2 = Dict{Int32,String}()
        writeproto(pb, testval)
        readproto(pb, readval)
        @test hasproperty(readval, :d2)
        assert_equal(testval.d2, readval.d2)

        testval = TestMaps()
        readval = TestMaps()
        testval.d3 = Dict{String,String}()
        writeproto(pb, testval)
        readproto(pb, readval)
        @test hasproperty(readval, :d3)
        assert_equal(testval.d3, readval.d3)

        testval = TestMaps()
        readval = TestMaps()
        testval.d1 = Dict{Int,Int}()
        testval.d1[1] = 1
        testval.d1[2] = 2
        writeproto(pb, testval)
        readproto(pb, readval)
        @test hasproperty(readval, :d1)
        assert_equal(testval, readval)

        testval = TestMaps()
        readval = TestMaps()
        testval.d2 = Dict{Int32,String}()
        testval.d2[Int32(1)] = convert(String, "One")
        testval.d2[Int32(2)] = convert(String, "Two")
        writeproto(pb, testval)
        readproto(pb, readval)
        @test hasproperty(readval, :d2)
        assert_equal(testval, readval)

        testval = TestMaps()
        readval = TestMaps()
        testval.d3 = Dict{String,String}()
        testval.d3["1"] = "One"
        testval.d3["2"] = "Two"
        writeproto(pb, testval)
        readproto(pb, readval)
        @test hasproperty(readval, :d3)
        assert_equal(testval, readval)
    end
end

function test_misc()
    @testset "Miscellaneous functionality" begin
        testfld = TestOptional(; sVal1=TestStr(; val="1"), sVal2=TestStr(; val=""), iVal2=Int64[1,2,3])
        readfld = TestOptional(; sVal1=TestStr(; val=""), sVal2=TestStr(; val="1"), iVal2=Int64[])
        copy!(readfld, testfld)
        assert_equal(readfld, testfld)

        tf = TestFilled()
        @test !isfilled(tf)
        TestTypeJType[] = AbstractString
        TestTypeWType[] = :string
        TestTypeFldNum[] = 1
        TestTypePack[] = ProtoBuf.DEF_PACK
        tf.fld1 = TestType(; val="")
        @test isfilled(tf)

        iob = IOBuffer()
        show(iob, meta(TestOptional))
        @test !isempty(take!(iob))
    end
end

function test_enums()
    @testset "Enums" begin
        @test getproperty(TestEnum, lookup(TestEnum, 0)) == TestEnum.UNIVERSAL
        @test getproperty(TestEnum, lookup(TestEnum, 1)) == TestEnum.WEB
        @test getproperty(TestEnum, lookup(TestEnum, 2)) == TestEnum.IMAGES
        @test getproperty(TestEnum, lookup(TestEnum, 3)) == TestEnum.LOCAL
        @test getproperty(TestEnum, lookup(TestEnum, 4)) == TestEnum.NEWS
        @test getproperty(TestEnum, lookup(TestEnum, 5)) == TestEnum.PRODUCTS
        @test getproperty(TestEnum, lookup(TestEnum, 6)) == TestEnum.VIDEO

        @test enumstr(TestEnum, TestEnum.LOCAL) == "LOCAL"
        @test_throws ErrorException enumstr(TestEnum, Int32(12))
    end
end

end # module ProtoBufTestCodec

@testset "Codec" begin
    ProtoBufTestCodec.test_types()
    ProtoBufTestCodec.test_enums()
    ProtoBufTestCodec.test_oneofs()
    ProtoBufTestCodec.test_maps()
    ProtoBufTestCodec.test_repeats()
    ProtoBufTestCodec.test_optional()
    ProtoBufTestCodec.test_nested()
    ProtoBufTestCodec.test_defaults()
    ProtoBufTestCodec.test_misc()
end


text:
# Alignment Anchor
# ================
#
# Sequence alignment anchor type and alignment data data structures.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"""
A type to store the operation enocded in an alignment (CIGAR Operations).
Also stores the position in the alignment view of the sequences, and the
corresponding position in the unaltered source sequence (nucleotide or protein).

It stores the position values as Int types and the alignment operation is stored
as a type `Operation`, these are defined in the file `operations.jl` and loaded
into Bio.Align namespace as a series of global constants.
"""
immutable AlignmentAnchor
    seqpos::Int
    refpos::Int
    op::Operation
end

function AlignmentAnchor(pos::Tuple{Int,Int}, op)
    return AlignmentAnchor(pos[1], pos[2], op)
end

function Base.show(io::IO, anc::AlignmentAnchor)
    print(io, "AlignmentAnchor(", anc.seqpos, ", ", anc.refpos, ", '", anc.op, "')")
end

# Alignment
# ---------

immutable Alignment
    anchors::Vector{AlignmentAnchor}
    firstref::Int
    lastref::Int

    function Alignment(anchors::Vector{AlignmentAnchor}, check::Bool=true)
        # optionally check coherence of the anchors
        if check
            # empty alignments are valid, representing an unaligned sequence
            if !isempty(anchors)
                if anchors[1].op != OP_START
                    error("Alignments must begin with on OP_START anchor.")
                end

                for i in 2:length(anchors)
                    if anchors[i].refpos < anchors[i-1].refpos ||
                       anchors[i].seqpos < anchors[i-1].seqpos
                        error("Alignment anchors must be sorted.")
                    end

                    op = anchors[i].op
                    if convert(UInt8, op) > convert(UInt8, OP_MAX_VALID)
                        error("Anchor at index $(i) has an invalid operation.")
                    end

                    # reference skip/delete operations
                    if isdeleteop(op)
                        if anchors[i].seqpos != anchors[i-1].seqpos
                            error("Invalid anchor positions for reference deletion.")
                        end
                    # reference insertion operations
                    elseif isinsertop(op)
                        if anchors[i].refpos != anchors[i-1].refpos
                            error("Invalid anchor positions for reference insertion.")
                        end
                    # match operations
                    elseif ismatchop(op)
                        if anchors[i].refpos - anchors[i-1].refpos !=
                           anchors[i].seqpos - anchors[i-1].seqpos
                            error("Invalid anchor positions for match operation.")
                        end
                    end
                end
            end
        end

        # compute first and last aligned reference positions
        firstref = 0
        for i in 1:length(anchors)
            if ismatchop(anchors[i].op)
                firstref = anchors[i-1].refpos + 1
                break
            end
        end

        lastref = 0
        for i in length(anchors):-1:1
            if ismatchop(anchors[i].op)
                lastref = anchors[i].refpos
                break
            end
        end

        return new(anchors, firstref, lastref)
    end
end

"""
Construct an `Alignment`
"""
function Alignment(cigar::AbstractString, seqpos::Int=1, refpos::Int=1)
    # path starts prior to the first aligned position pair
    seqpos -= 1
    refpos -= 1

    n = 0
    anchors = AlignmentAnchor[AlignmentAnchor(seqpos, refpos, OP_START)]
    for c in cigar
        if isdigit(c)
            n = n * 10 + convert(Int, c - '0')
        else
            if n == 0
                error("CIGAR operations must be prefixed by a positive integer.")
            end
            op = Operation(c)
            if ismatchop(op)
                seqpos += n
                refpos += n
            elseif isinsertop(op)
                seqpos += n
            elseif isdeleteop(op)
                refpos += n
            else
                error("The $(op) CIGAR operation is not yet supported.")
            end

            push!(anchors, AlignmentAnchor(seqpos, refpos, op))
            n = 0
        end
    end

    return Alignment(anchors)
end

@compat function Base.:(==)(a::Alignment, b::Alignment)
    return a.anchors == b.anchors && a.firstref == b.firstref && a.lastref == b.lastref
end

function Base.show(io::IO, aln::Alignment)
    # print a representation of the reference sequence
    anchors = aln.anchors
    for i in 2:length(anchors)
        if ismatchop(anchors[i].op)
            for _ in anchors[i-1].refpos+1:anchors[i].refpos
                write(io, '·')
            end
        elseif isinsertop(anchors[i].op)
            for _ in anchors[i-1].seqpos+1:anchors[i].seqpos
                write(io, '-')
            end
        elseif isdeleteop(anchors[i].op)
            for _ in anchors[i-1].refpos+1:anchors[i].refpos
                write(io, '·')
            end
        end
    end
    write(io, '\n')

    # print a representation of the aligned sequence
    for i in 2:length(anchors)
        if ismatchop(anchors[i].op)
            for i in anchors[i-1].seqpos+1:anchors[i].seqpos
                write(io, '·')
            end
        elseif isinsertop(anchors[i].op)
            for i in anchors[i-1].seqpos+1:anchors[i].seqpos
                write(io, '·')
            end
        elseif isdeleteop(anchors[i].op)
            for _ in anchors[i-1].refpos+1:anchors[i].refpos
                write(io, '-')
            end
        end
    end
end

# find the index of the first anchor that satisfies `i ≤ pos`
@generated function findanchor{isseq}(aln::Alignment, i::Integer, ::Type{Val{isseq}})
    pos = isseq ? :seqpos : :refpos
    quote
        anchors = aln.anchors
        lo = 1
        hi = endof(anchors)
        if !(anchors[lo].$pos < i ≤ anchors[hi].$pos)
            return 0
        end
        # binary search
        @inbounds while hi - lo > 2
            m = (lo + hi) >> 1
            if anchors[m].$pos < i
                lo = m
            else  # i ≤ anchors[m].$pos
                hi = m
            end
            # invariant (activate this for debugging)
            #@assert anchors[lo].$pos < i ≤ anchors[hi].$pos
        end
        # linear search
        @inbounds for j in lo+1:hi
            if i ≤ aln.anchors[j].$pos
                return j
            end
        end
        # do not reach here
        @assert false
        return 0
    end
end

"""
    seq2ref(i, aln)

Map a position from sequence to reference.
"""
function seq2ref(i::Integer, aln::Alignment)
    idx = findanchor(aln, i, Val{true})
    if idx == 0
        throw(ArgumentError("invalid sequence position: $i"))
    end
    anchor = aln.anchors[idx]
    refpos = anchor.refpos
    if ismatchop(anchor.op)
        refpos += i - anchor.seqpos
    end
    return refpos, anchor.op
end

"""
    ref2seq(i, aln)

Map a position from reference to sequence.
"""
function ref2seq(i::Integer, aln::Alignment)
    idx = findanchor(aln, i, Val{false})
    if idx == 0
        throw(ArgumentError("invalid reference position: $i"))
    end
    anchor = aln.anchors[idx]
    seqpos = anchor.seqpos
    if ismatchop(anchor.op)
        seqpos += i - anchor.refpos
    end
    return seqpos, anchor.op
end

# AlignedSequence
# ---------------

immutable AlignedSequence{S}
    seq::S
    aln::Alignment
end

function AlignedSequence{S}(seq::S, anchors::Vector{AlignmentAnchor},
                            check::Bool=true)
    return AlignedSequence(seq, Alignment(anchors, check))
end

function AlignedSequence(seq::Sequence, ref::Sequence)
    return AlignedSequence(seq, 1, ref, 1)
end

function AlignedSequence(seq::Sequence, seqpos::Integer,
                         ref::Sequence, refpos::Integer)
    if length(seq) != length(ref)
        throw(ArgumentError("two sequences must be the same length"))
    end
    seqpos -= 1
    refpos -= 1
    op = OP_START
    newseq = similar(seq, 0)  # sequence without gap symbols
    anchors = AlignmentAnchor[]
    for (x, y) in zip(seq, ref)
        if x == gap(eltype(seq)) && y == gap(eltype(ref))
            throw(ArgumentError("two sequences must not have gaps at the same position"))
        elseif x == gap(eltype(seq))
            op′ = OP_DELETE
        elseif y == gap(eltype(ref))
            op′ = OP_INSERT
        elseif x == y
            op′ = OP_SEQ_MATCH
        else
            op′ = OP_SEQ_MISMATCH
        end

        if op′ != op
            push!(anchors, AlignmentAnchor(seqpos, refpos, op))
            op = op′
        end

        if x != gap(eltype(seq))
            seqpos += 1
            push!(newseq, x)
        end
        if y != gap(eltype(ref))
            refpos += 1
        end
    end
    push!(anchors, AlignmentAnchor(seqpos, refpos, op))
    return AlignedSequence(newseq, anchors)
end

"""
First position in the reference sequence.
"""
function IntervalTrees.first(alnseq::AlignedSequence)
    return alnseq.aln.firstref
end

"""
Last position in the reference sequence.
"""
function IntervalTrees.last(alnseq::AlignedSequence)
    return alnseq.aln.lastref
end

function seq2ref(i::Integer, alnseq::AlignedSequence)
    return seq2ref(i, alnseq.aln)
end

function ref2seq(i::Integer, alnseq::AlignedSequence)
    return ref2seq(i, alnseq.aln)
end

# simple letters and dashes representation of an alignment
function Base.show(io::IO, alnseq::AlignedSequence)
    # print a representation of the reference sequence
    anchors = alnseq.aln.anchors
    for i in 2:length(anchors)
        if ismatchop(anchors[i].op)
            for _ in anchors[i-1].refpos+1:anchors[i].refpos
                write(io, '·')
            end
        elseif isinsertop(anchors[i].op)
            for _ in anchors[i-1].seqpos+1:anchors[i].seqpos
                write(io, '-')
            end
        elseif isdeleteop(anchors[i].op)
            for _ in anchors[i-1].refpos+1:anchors[i].refpos
                write(io, '·')
            end
        end
    end
    write(io, '\n')

    for i in 2:length(anchors)
        if ismatchop(anchors[i].op)
            for i in anchors[i-1].seqpos+1:anchors[i].seqpos
                print(io, alnseq.seq[i])
            end
        elseif isinsertop(anchors[i].op)
            for i in anchors[i-1].seqpos+1:anchors[i].seqpos
                print(io, alnseq.seq[i])
            end
        elseif isdeleteop(anchors[i].op)
            for _ in anchors[i-1].refpos+1:anchors[i].refpos
                write(io, '-')
            end
        end
    end
end

"""
Output a CIGAR string encoding of an `Alignment`. This is not entirely lossless as
it discards the alignments start positions.
"""
function cigar(aln::Alignment)
    anchors = aln.anchors
    out = IOBuffer()
    seqpos = anchors[1].seqpos
    refpos = anchors[1].refpos

    for i in 2:length(anchors)
        n = max(anchors[i].seqpos - anchors[i-1].seqpos,
                anchors[i].refpos - anchors[i-1].refpos)
        print(out, n, Char(anchors[i].op))
    end
    return takebuf_string(out)
end


text:
## labelmap

function labelmap!(dict::Dict{T,Vector{Int}}, idx::Int, elem) where {T}
    labelmap!(dict, idx, T(elem))
end

function labelmap!(dict::Dict{T,Vector{Int}}, idx::Int, elem::T) where {T}
    if !haskey(dict, elem)
        push!(dict, elem => [idx])
    else
        push!(dict[elem], idx)
    end
    dict
end

function labelmap!(dict::Dict{T,Vector{Int}}, idx::AbstractVector{Int}, itr::S) where {T,S}
    for (i, elem) in zip(idx, itr)
        labelmap!(dict, i, elem)
    end
    dict
end

function labelmap(itr::T) where {T}
    dict = Dict{eltype(T),Vector{Int}}()
    for (idx, elem) in enumerate(itr)
        labelmap!(dict, idx, elem)
    end
    dict
end

## labelfreq

function labelfreq!(dict::Dict{T,Int}, elem) where {T}
    labelfreq!(dict, T(elem))
end

function labelfreq!(dict::Dict{T,Int}, elem::T) where {T}
    cnt = get(dict, elem, 0)
    dict[elem] = cnt + 1
    dict
end

function labelfreq!(dict::Dict{T,Int}, iter::AbstractVector) where {T}
    for elem in iter
        labelfreq!(dict, T(elem))
    end
    dict
end

labelfreq(itr::T) where {T} = countmap(itr, alg = :dict)::Dict{eltype(T),Int}
labelfreq(dict::Dict{T,Vector{Int}})  where {T} = Dict(k => length(v) for (k,v) in dict)::Dict{T,Int}

## Errors

labelmap(A::AbstractMatrix)  = throw(ArgumentError("labelmap not supported for matrices"))
labelfreq(A::AbstractMatrix) = throw(ArgumentError("labelfreq not supported for matrices"))

## General Dict

labelenc(x::Dict) = labelenc(label(x))
nlabel(x::Dict) = length(keys(x))
label(x::Dict) = _arrange_label(collect(keys(x)))

## Convert label map to label vector

function labelmap2vec(lm::Dict{T, Vector{Int}}) where T
    isempty(lm) && return Vector{T}(undef, 0)
    labelvec = Vector{T}(undef, mapreduce(length, +, values(lm)))
    @inbounds for (k, v) in lm
        labelvec[v] .= k
    end
    return labelvec
end


text:
using SafeTestsets

##
@time @safetestset "Miscelaneous" begin
    include("./MiscTest.jl")
end

##
@time @safetestset "IO" begin
    include("./ioTest.jl")
end

##
@time @safetestset "Construct Dln" begin
    include("./ConstructNetworkTest.jl")
end

#
@time @safetestset "Segment Force" begin
    include("./SegmentForcesTest.jl")
end

##
@time @safetestset "Mobility" begin
    include("./MobilityTest.jl")
end

##
@time @safetestset "Topology" begin
    include("./TopologyTest.jl")
end

##
@time @safetestset "FEM" begin
    include("./FEMTest.jl")
end

##
@time @safetestset "Dislocation-FEM" begin
    include("./DislocationFEMTest.jl")
end

##
@time @safetestset "Post-process" begin
    include("./PlotTest.jl")
end

##
@time @safetestset "Integral tests" begin
    include("./IntegralTest.jl")
end

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
function compute_ERI!(I::IntegralHelper{T, Chonky, O}, aoints::IntegralHelper{T, Chonky, AtomicOrbitals}) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}
    AOERI = aoints["ERI"]
    C = I.orbitals.C
    if eltype(I.orbitals.C) !== T
        C = T.(C)
    end
    @tensoropt MOERI[i,j,k,l] :=  AOERI[μ, ν, ρ, σ]*C[μ, i]*C[ν, j]*C[ρ, k]*C[σ, l]
    I["ERI"] = MOERI
end

function compute_OOOO!(I::IntegralHelper{T,Chonky,O}, aoints::IntegralHelper{T,Chonky, AtomicOrbitals}) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}

    AOERI = aoints["ERI"]

    core = Options.get("drop_occ")
    ndocc = I.molecule.Nα
    o = (1+core):ndocc
    Co = I.orbitals.C[:,o]
    if eltype(I.orbitals.C) !== T
        Co = T.(Co)
    end
    @tensoropt (μ=>100x, ν=>100x, ρ=>100x, σ=>100x, i=>10x, j=>10x, k=>10x, l=>10x, a=>80x, b=>80x, c=>80x, d=>80) begin
        OOOO[i,j,k,l] :=  AOERI[μ, ν, ρ, σ]*Co[μ, i]*Co[ν, j]*Co[ρ, k]*Co[σ, l]
    end
    I["OOOO"] = OOOO
end

function compute_OOOV!(I::IntegralHelper{T,Chonky,O}, aoints::IntegralHelper{T,Chonky, AtomicOrbitals}) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}

    AOERI = aoints["ERI"]

    inac = Options.get("drop_vir")
    core = Options.get("drop_occ")
    ndocc = I.molecule.Nα
    nbf = size(I.orbitals.C,1)
    o = (1+core):ndocc
    v = (ndocc+1):(nbf - inac)
    Cv = I.orbitals.C[:,v]
    Co = I.orbitals.C[:,o]
    if eltype(I.orbitals.C) !== T
        Cv = T.(Cv)
        Co = T.(Co)
    end
    @tensoropt (μ=>100x, ν=>100x, ρ=>100x, σ=>100x, i=>10x, j=>10x, k=>10x, l=>10x, a=>80x, b=>80x, c=>80x, d=>80) begin
        OOOV[i,j,k,a] :=  AOERI[μ, ν, ρ, σ]*Co[μ, i]*Co[ν, j]*Co[ρ, k]*Cv[σ, a]
    end
    I["OOOV"] = OOOV
end

function compute_OOVV!(I::IntegralHelper{T,Chonky,O}, aoints::IntegralHelper{T,Chonky, AtomicOrbitals}) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}

    AOERI = aoints["ERI"]

    inac = Options.get("drop_vir")
    core = Options.get("drop_occ")
    ndocc = I.molecule.Nα
    nbf = size(I.orbitals.C,1)
    o = (1+core):ndocc
    v = (ndocc+1):(nbf - inac)
    Cv = I.orbitals.C[:,v]
    Co = I.orbitals.C[:,o]
    if eltype(I.orbitals.C) !== T
        Cv = T.(Cv)
        Co = T.(Co)
    end
    @tensoropt (μ=>100x, ν=>100x, ρ=>100x, σ=>100x, i=>10x, j=>10x, k=>10x, l=>10x, a=>80x, b=>80x, c=>80x, d=>80) begin
        OOVV[i,j,a,b] :=  AOERI[μ, ν, ρ, σ]*Co[μ, i]*Co[ν, j]*Cv[ρ, a]*Cv[σ, b]
    end
    I["OOVV"] = OOVV
end

function compute_OVOV!(I::IntegralHelper{T,Chonky,O}, aoints::IntegralHelper{T,Chonky, AtomicOrbitals}) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}

    AOERI = aoints["ERI"]

    inac = Options.get("drop_vir")
    core = Options.get("drop_occ")
    ndocc = I.molecule.Nα
    nbf = size(I.orbitals.C,1)
    o = (1+core):ndocc
    v = (ndocc+1):(nbf - inac)
    Cv = I.orbitals.C[:,v]
    Co = I.orbitals.C[:,o]
    if eltype(I.orbitals.C) !== T
        Cv = T.(Cv)
        Co = T.(Co)
    end
    @tensoropt (μ=>100x, ν=>100x, ρ=>100x, σ=>100x, i=>10x, j=>10x, k=>10x, l=>10x, a=>80x, b=>80x, c=>80x, d=>80) begin
        OVOV[i,a,j,b] :=  AOERI[μ, ν, ρ, σ]*Co[μ, i]*Cv[ν, a]*Co[ρ, j]*Cv[σ, b]
    end
    I["OVOV"] = OVOV
end

function compute_OVVV!(I::IntegralHelper{T,Chonky,O}, aoints::IntegralHelper{T,Chonky, AtomicOrbitals}) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}

    AOERI = aoints["ERI"]

    inac = Options.get("drop_vir")
    core = Options.get("drop_occ")
    ndocc = I.molecule.Nα
    nbf = size(I.orbitals.C,1)
    o = (1+core):ndocc
    v = (ndocc+1):(nbf - inac)
    Cv = I.orbitals.C[:,v]
    Co = I.orbitals.C[:,o]
    if eltype(I.orbitals.C) !== T
        Cv = T.(Cv)
        Co = T.(Co)
    end
    @tensoropt (μ=>100x, ν=>100x, ρ=>100x, σ=>100x, i=>10x, j=>10x, k=>10x, l=>10x, a=>80x, b=>80x, c=>80x, d=>80) begin
        OVVV[i,a,b,c] :=  AOERI[μ, ν, ρ, σ]*Co[μ, i]*Cv[ν, a]*Cv[ρ, b]*Cv[σ, c]
    end
    I["OVVV"] = OVVV
end

function compute_VVVV!(I::IntegralHelper{T,Chonky,O}, aoints::IntegralHelper{T,Chonky, AtomicOrbitals}) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}

    AOERI = aoints["ERI"]

    inac = Options.get("drop_vir")
    ndocc = I.molecule.Nα
    nbf = size(I.orbitals.C,1)
    v = (ndocc+1):(nbf - inac)
    Cv = I.orbitals.C[:,v]
    if eltype(I.orbitals.C) !== T
        Cv = T.(Cv)
    end
    @tensoropt (μ=>100x, ν=>100x, ρ=>100x, σ=>100x, i=>10x, j=>10x, k=>10x, l=>10x, a=>80x, b=>80x, c=>80x, d=>80) begin
        VVVV[a,b,c,d] :=  AOERI[μ, ν, ρ, σ]*Cv[μ, a]*Cv[ν, b]*Cv[ρ, c]*Cv[σ, d]
    end
    I["VVVV"] = VVVV
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
using ViscosityDrop

### Setup simulation and run
filename = "example"
dir = "."
Ny, Nz = 64, 64
do_all(filename, dir, Ny=Ny, Nz=Nz, yLength=2, zLength=2,
                      ν=1e-1,
                      z_position=0.6, width=0.2, amplitude=-1e-1,
                      Δt=0.002, stop_time=120.0, iteration_interval=1000,
                      schedule_interval=0.5, output_fields = :all)
### Post-process (also saves to txt file)
out_results = get_results(dir, dir, "example.jld2"; Ncourse=11, verbose=false)

### Plot results
using Plots

p1 = plot(out_results[:,1],out_results[:,2],
          label="test1", title="Circle radius");
p2 = plot(out_results[:,1],out_results[:,3],
          label="test1", title="y-position");
p3 = plot(out_results[:,1],out_results[:,4],
          label="test1", title="z-position");
p4 = plot(out_results[:,1],out_results[:,5],
          label="test1", title="Circle misfit");
p_all = plot(p1, p2, p3, p4, layout = (2, 2), legend = false)
savefig(p_all, "results.pdf")


text:
@testset "Primes" begin
  @testset "Next prime" begin
    for T in [UInt32, UInt64, UInt128, Int32, Int64, Int128, BigInt, fmpz]
      @test @inferred next_prime(T(0)) == T(2)
      @test @inferred next_prime(T(2)) == T(3)
      @test @inferred next_prime(T(3)) == T(5)
    end

    if Int == Int64
      @test_throws InexactError next_prime(Int(9223372036854775783))
    elseif Int == Int32
      @test_throws InexactError next_prime(Int32(2147483647))
    end

    for B in 1:100
      for i in 1:10
        z = rand(fmpz(10)^B:fmpz(10)^(B + 1))
        p = @inferred next_prime(z)
        @test @inferred isprime(p)
      end
    end
  end

  @testset "PrimesSet" begin
    for T in [Int, BigInt, fmpz]
      P = @inferred PrimesSet(T(0), T(100))
      PP = collect(P)
      @test length(PP) == 25
      @test all(isprime(p) for p in P)
      @test PP[1] == T(2)
      @test PP[end] == T(97)

      P = @inferred PrimesSet(T(2), T(100))
      PP = collect(P)
      @test length(PP) == 25
      @test all(isprime(p) for p in P)
      @test PP[1] == T(2)
      @test PP[end] == T(97)

      P = @inferred collect(PrimesSet(T(100), T(1000), T(5), T(1)))
      PP = collect(P)
      @test length(PP) == 35
      @test PP[1] == T(101)
      @test PP[end] == T(991)
      @test all(isprime(p) && iszero(mod(p - T(1), T(5))) for p in P)

      P = @inferred PrimesSet(T(100), T(1000))
      PP = collect(P)
      @test length(PP) == 143
      @test PP[1] == T(101)
      @test PP[end] == T(997)

      for i in 1:10
        B = rand(100:200)
        modd = rand(2:20)
        a = rand(0:modd)
        while gcd(modd, a) != 1
          a = rand(0:modd)
        end
        P = @inferred PrimesSet(T(B), 2 * T(B), T(modd), T(a))
        @test all(isprime(p) && iszero(mod(p - T(a), T(modd))) for p in P)
      end
    end

    for B in 20:25
      P = @inferred PrimesSet(fmpz(10)^B, fmpz(10)^B + fmpz(10)^5)
      @test all(isprime(p) for p in P)
    end
  end

  @testset "Factorization" begin
    z = fmpz(35406006584838250779859646288524045541522102721903852419351442034218745921487618213430210114517205290360200098727562436996683618008616544114512403804879033305568279270240314426184896768288124359592373197652123675245898751520012787598855558789031381024625387676970030570435823658379242400351150075732901033905965864786338976200877211640473166089014964716251974918808271240247178842792601062285288192064508461909781716550249247324239984258941121466497704357904580222495174166761944293422767709291135213610316677755849561319492135401953603941166886125617969544902674910709468700828636641323083845265496724211468938267526059318324623810616578510536586846196125620883078603703124017829344050945482120041075683588528213410512026977220472164498442361107233172928990996889702611643333700128616817177297106110731760700226123826960770570905079841104576624248095584085121616172787911816790002442329981176247943435649)
    fac = factor(z)
    @test unit(fac) * prod(p^e for (p, e) in fac) == z
  end
end


text:
#TODO: add "expand" method to go from shorter to greater frequencies
#TODO: include functionality for filling the missing values (approx, locf, etc.)

# Date vector
mondays(t::AbstractArray{Date,1}) = dayofweek(t) .== Monday
tuesdays(t::AbstractArray{Date,1}) = dayofweek(t) .== Tuesday
wednesdays(t::AbstractArray{Date,1}) = dayofweek(t) .== Wednesday
thursdays(t::AbstractArray{Date,1}) = dayofweek(t) .== Thursday
fridays(t::AbstractArray{Date,1}) = dayofweek(t) .== Friday
saturdays(t::AbstractArray{Date,1}) = dayofweek(t) .== Saturday
sundays(t::AbstractArray{Date,1}) = dayofweek(t) .== Sunday
# DateTime vector
mondays(t::AbstractArray{DateTime,1}) = dayofweek(t) .== Monday
tuesdays(t::AbstractArray{DateTime,1}) = dayofweek(t) .== Tuesday
wednesdays(t::AbstractArray{DateTime,1}) = dayofweek(t) .== Wednesday
thursdays(t::AbstractArray{DateTime,1}) = dayofweek(t) .== Thursday
fridays(t::AbstractArray{DateTime,1}) = dayofweek(t) .== Friday
saturdays(t::AbstractArray{DateTime,1}) = dayofweek(t) .== Saturday
sundays(t::AbstractArray{DateTime,1}) = dayofweek(t) .== Sunday
# Time series input
mondays(x::TS) = x[mondays(x.index)]
tuesdays(x::TS) = x[tuesdays(x.index)]
wednesdays(x::TS) = x[wednesdays(x.index)]
thursdays(x::TS) = x[thursdays(x.index)]
fridays(x::TS) = x[fridays(x.index)]
saturdays(x::TS) = x[saturdays(x.index)]
sundays(x::TS) = x[sundays(x.index)]

# If `cal` is true, only returns observations on the last calendar day of the month.
# If `cal` is false, returns all observations for which the next index is a new month.
function bow(t::AbstractArray{Date,1}; cal::Bool=false)
    if cal
        return BitArray{1}(t .== firstdayofweek(t))
    else
        return BitArray{1}([0;diff(week(t)).!=0])
    end
end
function eow(t::AbstractArray{Date,1}; cal::Bool=false)
    if cal
        return BitArray{1}(t .== lastdayofweek(t))
    else
        return BitArray{1}([diff(week(t)).!=0;0])
    end
end
function bom(t::AbstractArray{Date,1}; cal::Bool=false)
    if cal
        return BitArray{1}(t .== firstdayofmonth(t))
    else
        return BitArray{1}([0;diff(month(t)).!=0])
    end
end
function eom(t::AbstractArray{Date,1}; cal::Bool=false)
    if cal
        return BitArray{1}(t .== lastdayofmonth(t))
    else
        BitArray{1}([diff(month(t)).!=0;0])
    end
end
function boq(t::AbstractArray{Date,1}; cal::Bool=false)
    if cal
        return BitArray{1}(t .== firstdayofquarter(t))
    else
        return BitArray{1}([0;diff(quarterofyear(t)).!=0])
    end
end
function eoq(t::AbstractArray{Date,1}; cal::Bool=false)
    if cal
        return BitArray{1}(t .== lastdayofquarter(t))
    else
        return BitArray{1}([diff(quarterofyear(t)).!=0;0])
    end
end
function boy(t::AbstractArray{Date,1}; cal::Bool=false)
    if cal
        return BitArray{1}(t .== firstdayofyear(t))
    else
        return BitArray{1}([0;diff(year(t)).!=0])
    end
end
function eoy(t::AbstractArray{Date,1}; cal::Bool=false)
    if cal
        return BitArray{1}(t .== lastdayofyear(t))
    else
        return BitArray{1}([diff(year(t)).!=0;0])
    end
end
function bow(t::AbstractArray{DateTime,1}; cal::Bool=false)
    if cal
        return BitArray{1}(t .== firstdayofweek(t))
    else
        return BitArray{1}([0;diff(week(t)).!=0])
    end
end
function eow(t::AbstractArray{DateTime,1}; cal::Bool=false)
    if cal
        return BitArray{1}(t .== lastdayofweek(t))
    else
        return BitArray{1}([diff(week(t)).!=0;0])
    end
end
function bom(t::AbstractArray{DateTime,1}; cal::Bool=false)
    if cal
        return BitArray{1}(t .== firstdayofmonth(t))
    else
        return BitArray{1}([0;diff(month(t)).!=0])
    end
end
function eom(t::AbstractArray{DateTime,1}; cal::Bool=false)
    if cal
        return BitArray{1}(t .== lastdayofmonth(t))
    else
        BitArray{1}([diff(month(t)).!=0;0])
    end
end
function boq(t::AbstractArray{DateTime,1}; cal::Bool=false)
    if cal
        return BitArray{1}(t .== firstdayofquarter(t))
    else
        return BitArray{1}([0;diff(quarterofyear(t)).!=0])
    end
end
function eoq(t::AbstractArray{DateTime,1}; cal::Bool=false)
    if cal
        return BitArray{1}(t .== lastdayofquarter(t))
    else
        return BitArray{1}([diff(quarterofyear(t)).!=0;0])
    end
end
function boy(t::AbstractArray{DateTime,1}; cal::Bool=false)
    if cal
        return BitArray{1}(t .== firstdayofyear(t))
    else
        return BitArray{1}([0;diff(year(t)).!=0])
    end
end
function eoy(t::AbstractArray{DateTime,1}; cal::Bool=false)
    if cal
        return BitArray{1}(t .== lastdayofyear(t))
    else
        return BitArray{1}([diff(year(t)).!=0;0])
    end
end

bow(x::TS; cal::Bool=false) = x[bow(x.index, cal=cal)]
eow(x::TS; cal::Bool=false) = x[eow(x.index, cal=cal)]
bom(x::TS; cal::Bool=false) = x[bom(x.index, cal=cal)]
eom(x::TS; cal::Bool=false) = x[eom(x.index, cal=cal)]
boq(x::TS; cal::Bool=false) = x[boq(x.index, cal=cal)]
eoq(x::TS; cal::Bool=false) = x[eoq(x.index, cal=cal)]
boy(x::TS; cal::Bool=false) = x[boy(x.index, cal=cal)]
eoy(x::TS; cal::Bool=false) = x[eoy(x.index, cal=cal)]

function collapse{V,T}(x::TS{V,T}, at::AbstractArray{Bool,1}; fun::Function=last, args...)
    @assert size(at,1) == size(x,1) "Arguments `x` and `at` must have same number of rows."
    idx = find(at)
    if fun == last
        return x[idx]
    end
    out = zeros(V, (length(idx)-1, size(x,2)))
    a, b = (idx[1:end-1], idx[2:end])
    @inbounds for j=1:size(x,2), i=1:length(idx)-1
        out[i,j] = fun(x.values[a[i]:b[i], j]; args...)
    end
    at[idx[1]] = false
    return ts(out, x.index[at], x.fields)
end

# Here `at` is a function of the index of `x` which generates a BitArray to be used as above
function collapse{V,T}(x::TS{V,T}, at::Function; fun::Function=last, args...)
    return collapse(x, at(x.index), fun=fun; args...)
end

function apply{V}(x::TS{V}, dim::Int=1; fun=sum)
    if dim == 1
        n = size(x,1)
        out = zeros(V, n)
        @inbounds for i in 1:n
            out[i] = fun(x.values[i,:])
        end
        return ts(out, x.index, Symbol(ucfirst(string(fun))))
    elseif dim == 2
        k = size(x,2)
        out = zeros(V, k)
        @inbounds for j in 1:k
            out[j] = fun(x.values[:,j])
        end
        return out'
    else
        error("Argument `dim` must be either 1 (rows) or 2 (columns).")
    end
end


text:
using CSV
using DelimitedFiles
using DSP
using Random
using LinearAlgebra

function loadtestdata(testCase::String; dataPath::String=(@__DIR__)*"/data_new_testcases")
    #println(pwd()) # to debug github action
    data = readdlm(joinpath(dataPath, "$(testCase)_data.csv"), ',', Float64, '\n')
    data = dropdims(data,dims=1) # convert to vector
    evts = CSV.read(joinpath(dataPath, "$(testCase)_events.csv"), DataFrame)
    return data, evts
end

function gridexpand(conditionA  = [0,1.],continuousA = [-1.,0,1.])

tmp = reshape([ [x,y]  for x=conditionA, y=continuousA ],length(conditionA)*length(continuousA))
evts_grid = DataFrame(hcat(tmp...)')
rename!(evts_grid,["conditionA","continuousA"])
return evts_grid
end

simulate_lmm(args...; kwargs...) = simulate_lmm(Random.GLOBAL_RNG, args...; kwargs...)

function simulate_lmm(rng::AbstractRNG, τ = 1.5,fs = 12;
     β=[0.,-1.],
     σs=[1,1,1],
     σ=1,
     n_sub = 20,
     n_item = 30,
     noise_type="AR-exponential")
    rng_copy = deepcopy(rng)

subj_btwn = item_btwn = both_win = nothing
#subj_btwn = Dict("age" => ["O", "Y"])

# there are no between-item factors in this design so you can omit it or set it to nothing
item_btwn = Dict("stimType" => ["I","II"])

# put within-subject/item factors in a Dict
#both_win = Dict("condition" => ["A", "B"])

# simulate data
evt = DataFrame(simdat_crossed(n_sub, n_item,
                    subj_btwn = subj_btwn,
                    item_btwn = item_btwn,
                    both_win = both_win))

#    f1 = @formula dv ~ 1 + age * condition  + (1+condition|item) + (1+condition|subj);
f1 = @formula dv ~ 1 + stimType  + (1+stimType|subj) + (1|item);
m = MixedModels.fit(MixedModel, f1, evt)

# set the random effects

#gen_han(τ,fs,1)
basis = gen_han(τ,fs,2)

epoch_dat = zeros(Int(τ*fs),size(evt,1))

    #MixedModels doesnt really support σ==0 because all Ranef are scaled residual variance
σ_lmm = 0.0001
σs = σs./σ_lmm

for t = 1:size(epoch_dat,1)
    b = basis[t]

    MixedModelsSim.update!(m,create_re(b .* σs[1],b .*σs[2]),create_re( b.* σs[3]))
    simulate!(deepcopy(rng_copy),m,β = [b .* β[1] ,b .* β[2]], σ = σ_lmm)

    epoch_dat[t,:] = m.y
end

epoch_dat = reshape(epoch_dat,(1,size(epoch_dat)...))

# add some noise
if noise_type == "normal"
    epoch_dat = epoch_dat .+ randn(rng,size(epoch_dat)) .*σ
elseif noise_type == "AR-exponential"
    epoch_dat = epoch_dat .+ gen_noise_exp(rng,size(epoch_dat)).*σ

end

return evt,epoch_dat
end

function gen_han(τ,fs,peak)
    hanLen = Int(τ*fs/3)
    han = hanning(hanLen,zerophase=false)
    sig = zeros(Int(τ*fs))
    sig[1+hanLen*(peak-1):hanLen*peak] .= han
    return sig
end

function circulant(x)
    # Author: Jaromil Frossard
    # returns a symmetric matrix where X was circ-shifted.
    lx = length(x)
    ids = [1:1:(lx-1);]
    a = Array{Float64,2}(undef, lx,lx)
    for i = 1:length(x)
        if i==1
            a[i,:] = x
        else
            a[i,:] = vcat(x[i],a[i-1,ids])
        end
    end
    return Symmetric(a)
end

function exponentialCorrelation(x; nu = 1, length_ratio = 1)
    # Author: Jaromil Frossard
    # generate exponential function
    R = length(x)*length_ratio
    return exp.(-3*(x/R).^nu)
end

function noise_exp(rng,n_t)
    Σ = circulant(exponentialCorrelation([0:1:(n_t-1);],nu=1))
    Ut =  LinearAlgebra.cholesky(Σ).U'
    return (randn(rng,n_t)'*Ut')[1,:]
end

function gen_noise_exp(rng,si)
    n = hcat(map(x->noise_exp(rng,si[2]),1:si[3])...)
    return reshape(n,si)
end


text:
"""
Increments the exponent of a number, without altering it's final value.
The result is a non-normalized OraNumber.

# Example

Given a number "2244" = 22.44 x 100e1 (mantissa 22 44, exponent 1)
If `by_exponent = 2`, the result will be 00.002244 x 100e3 ( 00 00 22 44, exponent 3 )
"""
@inline function inc_exponent(num::OraNumber, by_exponent::Int8) :: OraNumber
    @assert by_exponent >= 0

    if by_exponent == Int8(0)
        return num
    end

    if iszero(num) || is_negative_1e126(num)
        error("inc_exponent not supported for $num.")
    end

    new_len = num.len + by_exponent

    if new_len > UInt8(21)
        error("Overflow on inc_exponent $num by $by_exponent")
    end

    old_exponent_byte = decode_exponent_byte(num)
    new_exponent_byte = old_exponent_byte + by_exponent
    if new_exponent_byte < old_exponent_byte
        error("Overflow on inc_exponent $num by $by_exponent.")
    end

    _isnegative = isnegative(num)
    new_exponent_encoded = encode_exponent_byte(new_exponent_byte, _isnegative)

    @inline function f(i)
        if i <= by_exponent
            return encode_mantissa_byte(0x00, _isnegative)
        else
            return num.mantissa[i-by_exponent]
        end
    end

    new_mantissa = Tuple( f(i) for i in 1:MAX_SIZEOF_MANTISSA )

    return OraNumber(new_len, new_exponent_encoded, new_mantissa)
end

@inline inc_exponent(num::OraNumber, by_exponent::Integer) = inc_exponent(num, Int8(by_exponent))

@inline function isnormalized(num::OraNumber) :: Bool
    if sizeof_mantissa(num) == 0x00
        # 0 of -1e126
        return true
    end

    # the number is considered normalized when the first mantissa byte is not zero
    return decode_mantissa_byte(num, 1) != 0x00
end

@inline function normalize(num::OraNumber) :: OraNumber
    # gets the number of left-hand-size zeros in the mantissa
    local lhs_mantissa_zero_bytes::Int8 = Int8(0)
    _isnegative = isnegative(num)
    _sizeof_mantissa = sizeof_mantissa(num)
    for i in 1:_sizeof_mantissa
        if decode_mantissa_byte(num, i, _isnegative) != 0x00
            break
        end

        lhs_mantissa_zero_bytes += Int8(1)
        @assert lhs_mantissa_zero_bytes != 0x00 "Overflow on normalizing $num."
    end

    if lhs_mantissa_zero_bytes == Int8(0)
        # num is already normalized
        return num
    else
        # will decrease exponent by `lhs_mantissa_zero_bytes`
        new_len = num.len - lhs_mantissa_zero_bytes
        old_exponent_byte = decode_exponent_byte(num)
        new_exponent_byte = old_exponent_byte - lhs_mantissa_zero_bytes

        if new_exponent_byte > old_exponent_byte
            error("Underflow on normalizing $num.")
        end

        new_exponent_encoded = encode_exponent_byte(new_exponent_byte, _isnegative)

        @inline function f(i)
            index = i + lhs_mantissa_zero_bytes
            if index > MAX_SIZEOF_MANTISSA
                return 0x00
            else
                return num.mantissa[i + lhs_mantissa_zero_bytes]
            end
        end

        new_mantissa = Tuple( f(i) for i in 1:MAX_SIZEOF_MANTISSA )
        return OraNumber(new_len, new_exponent_encoded, new_mantissa)
    end
end

function Base.:(-)(num::OraNumber) :: OraNumber
    if iszero(num)
        return num
    elseif is_negative_1e126(num)
        error("Not supported.")
    end

    _isnegative = isnegative(num)
    new_exponent = encode_exponent_byte(decode_exponent_byte(num), !_isnegative)
    _sizeof_mantissa = sizeof_mantissa(num)

    @inline function f(i)
        if i > _sizeof_mantissa
            if !_isnegative && i == _sizeof_mantissa + 1
                # add trailing byte, because we can
                return TRAILING_BYTE_ON_NEGATIVE_NUMBERS
            else
                return 0x00
            end
        else
            return encode_mantissa_byte( decode_mantissa_byte(num, i, _isnegative) , !_isnegative )
        end
    end

    new_mantissa = Tuple(f(i) for i in 1:MAX_SIZEOF_MANTISSA)

    if _isnegative
        # new number is positive
        if num.mantissa[num.len-UInt8(1)] == TRAILING_BYTE_ON_NEGATIVE_NUMBERS
            # old number has a trailing byte. Let's remove it
            return OraNumber(num.len - UInt8(1), new_exponent, new_mantissa)
        end
    else
        # new number is negative
        if num.len < MAX_SIZEOF_MANTISSA + 1
            # there's room for a trailing byte. Let's add it, given that `f` added it to the mantissa.
            return OraNumber(num.len + UInt8(1), new_exponent, new_mantissa)
        end
    end

    return OraNumber(num.len, new_exponent, new_mantissa)
end


text:
# Define `Recyclers.@global`
macro _global(ex::Expr)
    if Meta.isexpr(ex, :(=), 2)
        lhs, rhs = ex.args
        if lhs isa Symbol
            assignment = :(const $(esc(lhs)) = $(esc(rhs)))
            init = :(Recyclers.unsafe_init!($(esc(lhs))))
            InitializerModule = esc(gensym(:InitializerModule))
            expr = quote
                module $InitializerModule
                function $(esc(:__init__)) end
                end  # module $InitializerModule
                $(esc(:_M)) = $InitializerModule
                $(Expr(:block, __source__, assignment))
                $InitializerModule.__init__() = $(Expr(:block, __source__, init))
            end
            return Expr(:toplevel, expr.args...)
        end
    end
    error("require an assignment to a symbol (global name)")
end


text:
using Laguerre
using Test
using LinearAlgebra

include("test_lpoly.jl")
include("test_lfunction.jl")
include("test_poly_fun.jl")
include("test_integration.jl")


text:
function little_2var()
    v = Var(2)
    pos = compile(PlainLogicCircuit, var2lit(v))
    neg = compile(PlainLogicCircuit, -var2lit(v))
    or1 = pos | neg
    or2 = pos | neg

    v = Var(1)
    pos = compile(PlainLogicCircuit, var2lit(v))
    neg = compile(PlainLogicCircuit, -var2lit(v))

    and1 = pos & or1
    and2 = neg & or2
    and1 | and2
end

function little_3var()
    or1 = little_2var()
    v = Var(3)

    pos = compile(PlainLogicCircuit,  var2lit(v))
    neg = compile(PlainLogicCircuit, -var2lit(v))

    or2 = disjoin(children(or1))

    and1 = pos & or1
    and2 = neg & or2
    and1 | and2
end

function little_3var_constants()
    or1 = little_2var()
    v = Var(3)

    t = compile(PlainLogicCircuit, true)
    f = compile(PlainLogicCircuit, false)

    pos = compile(PlainLogicCircuit,  var2lit(v)) & t
    neg = compile(PlainLogicCircuit, -var2lit(v)) & f

    or2 = disjoin(children(or1))

    and1 = pos & or1
    and2 = neg & or2
    and1 | and2
end

function little_4var()
    ors = map(1:4) do v
        v = Var(v)
        pos = compile(PlainLogicCircuit, var2lit(v))
        neg = compile(PlainLogicCircuit, - var2lit(v))
        or = pos | neg
    end
    and1 = ors[1] & ors[2]
    and2 = ors[3] & ors[4]
    or = and1 | and2
end

function little_5var()
    c_4var = little_4var()
    v = Var(5)
    pos = compile(PlainLogicCircuit, var2lit(v))
    neg = compile(PlainLogicCircuit, - var2lit(v))
    or = pos | neg
    and = c_4var & or
    Plain⋁Node([and])
end

function readme_sdd()
    manager = SddMgr(7, :balanced)

    sun, rain, rainbow, cloud, snow, los_angeles, belgium = pos_literals(Sdd, manager, 7)

    sdd = (rainbow & sun & rain) | (-rainbow)
    sdd &= (-los_angeles | -belgium)
    sdd &= (los_angeles ⇒ sun) ∧ (belgium ⇒ cloud)
    sdd &= (¬(rain ∨ snow) ⇐ ¬cloud)
    sdd
end

max_stars_repo_name: Juice-jl/DiscriminativeCircuits

max_issues_repo_name: Juice-jl/DiscriminativeCircuits

max_forks_repo_name: Juice-jl/DiscriminativeCircuits


text:
# Autogenerated wrapper script for Qt5QuickControls_jll for x86_64-apple-darwin
export quickcontrols

using Qt5Declarative_jll
JLLWrappers.@generate_wrapper_header("Qt5QuickControls")
JLLWrappers.@declare_file_product(quickcontrols)
function __init__()
    JLLWrappers.@generate_init_header(Qt5Declarative_jll)
    JLLWrappers.@init_file_product(
        quickcontrols,
        "qml/QtQuick/Controls",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
abstract type Kernel end

struct TopHat <: Kernel end
struct BottomHat <: Kernel  end
struct LocalMean <: Kernel end
struct Erode <: Kernel end
struct Dilate <: Kernel end

max_stars_repo_name: MKAbdElrahman/TopologyOptimizationTools

max_issues_repo_name: MKAbdElrahman/TopologyOptimizationTools

max_forks_repo_name: MKAbdElrahman/TopologyOptimizationTools