text:
module SSmodule
    using LinearAlgebra
    function eigs_SS(A,x)
    end

    function shiftedbicg_Frommer2003(σ,A,b;maxsteps=3000,eps = 1e-15)
        N = length(σ)
        btype = typeof(b)
        q = zero(b)
        x = zero(b)
        r = deepcopy(b)
        rt = deepcopy(b)
        u = zero(b)
        ut = zero(b)
        vec_u = Vector{btype}(undef,N)
        vec_x = Vector{btype}(undef,N)
        for j=1:N
            vec_u[j] = zero(b)
            vec_x[j] = zero(b)
        end
        ρkold = ones(ComplexF64,N)
        ρk = ones(ComplexF64,N)
        πkold = ones(ComplexF64,N)
        πk = ones(ComplexF64,N)
        πknew = ones(ComplexF64,N)
        ρold = 1
        αold = 1
        residual = 0

        for i=1:maxsteps
            ρ = dot(rt,r)
            β = -ρ/ρold
            axpby!(1,r,-β,u)
            axpby!(1,rt,-β',ut)
            mul!(q,A,u)
            qrt = dot(rt,q)
            α = ρ/qrt
            axpy!(α,u,x)

            for k=1:N
                πknew[k] = (1+α*σ[k])*πk[k]+ (α*β/αold)*(πkold[k]- πk[k])
                βk = (πkold[k]/πk[k])^2*β
                αk =  (πk[k]/πknew[k])*α
                axpby!(1/πk[k],r,-βk,vec_u[k])
                axpy!(αk,vec_u[k],vec_x[k])
            end

            axpy!(-α,q,r)
            mul!(q,A',ut)
            axpy!(-α',q,rt)
            αold = α
            ρold  = ρ

            residual = dot(r,r)
            println(dot(r,r))

            if abs(residual) < eps
                println("Converged at $i-th step. eps: $residual")
                println("--------------------------------------")
                return vec_x
            end

            for k=1:N
                πkold[k] = πk[k]
                πk[k] = πknew[k]
            end

        end

        error("bicg is not converged. The residual is $residual")

    end

    function shiftedbicg_Frommer2003_seed(σ,A,b;maxsteps=3000,eps = 1e-15,zseedin = 0im)
        zseed = zseedin
        N = length(σ)
        btype = typeof(b)
        q = zero(b)
        x = zero(b)
        r = deepcopy(b)
        rt = deepcopy(b)
        u = zero(b)
        ut = zero(b)
        vec_u = Vector{btype}(undef,N)
        vec_x = Vector{btype}(undef,N)
        for j=1:N
            vec_u[j] = zero(b)
            vec_x[j] = zero(b)
        end
        ρkold = ones(ComplexF64,N)
        ρk = ones(ComplexF64,N)
        πkold = ones(ComplexF64,N)
        πk = ones(ComplexF64,N)
        πknew = ones(ComplexF64,N)
        ρold = 1
        αold = 1
        residual = 0

        for i=1:maxsteps
            ρ = dot(rt,r)
            β = -ρ/ρold
            axpby!(1,r,-β,u)
            axpby!(1,rt,-β',ut)
            mul!(q,A,u)
            axpby!(zseed,u,1,q)

            qrt = dot(rt,q)
            α = ρ/qrt
            axpy!(α,u,x)

            for k=1:N
                πknew[k] = (1+α*(σ[k]-zseed))*πk[k]+ (α*β/αold)*(πkold[k]- πk[k])
                βk = (πkold[k]/πk[k])^2*β
                αk =  (πk[k]/πknew[k])*α
                axpby!(1/πk[k],r,-βk,vec_u[k])
                axpy!(αk,vec_u[k],vec_x[k])
            end

            axpy!(-α,q,r)
            mul!(q,A',ut)
            axpby!(zseed',ut,1,q)
            axpy!(-α',q,rt)

            residual = dot(r,r)
            println(dot(r,r),"\t",zseed)

            for k=1:N
                πkold[k] = πk[k]
                πk[k] = πknew[k]
            end

            if false

            πkseed = 1
            πkseedold = 1
            kseed = 0
            for k=1:N
                if abs(πk[k]) < abs(πkseed)
                    πkseed =πk[k]
                    πkseedold = πkold[k]
                    zseed = σ[k]
                    kseed = k
                end
            end

            println(kseed)
            if kseed != 0
                axpby!(0,rt,1/πkseed,r)
                axpby!(0,r,1/πkseed,rt)
                axpby!(0,r,1/πkseed,u)
                axpby!(0,r,1/πkseed,ut)
                axpby!(1,vec_x[kseed],0,x)
                ρ = ρ/(πkseedold*πkseedold)
                α = (πkseedold/πkseed)*α

                for k=1:N
                    πkold[k] = πkold[k]/πkseedold
                    πk[k] = πk[k]/πkseed
                end
            end

            end

            αold = α
            ρold  = ρ

            if abs(residual) < eps
                println("Converged at $i-th step. eps: $residual")
                println("--------------------------------------")
                return vec_x
            end

        end

        error("bicg is not converged. The residual is $residual")

    end

    function shiftedbicg_Frommer2003_G_seed(σ,A,b;maxsteps=3000,eps = 1e-15,zseedin = 0im)
        zseed = zseedin
        N = length(σ)
        btype = typeof(b)
        q = zero(b)
        x = zero(b)
        r = deepcopy(b)
        rt = deepcopy(b)
        u = zero(b)
        ut = zero(b)
        vec_u = Vector{btype}(undef,N)
        vec_x = Vector{btype}(undef,N)
        for j=1:N
            vec_u[j] = zero(b)
            vec_x[j] = zero(b)
        end
        ρkold = ones(ComplexF64,N)
        ρk = ones(ComplexF64,N)
        πkold = ones(ComplexF64,N)
        πk = ones(ComplexF64,N)
        πknew = ones(ComplexF64,N)
        ρold = 1
        αold = 1
        residual = 0

        for i=1:maxsteps
            ρ = dot(rt,r)
            β = -ρ/ρold
            axpby!(1,r,-β,u)
            axpby!(1,rt,-β',ut)
            mul!(q,A,u)
            axpby!(zseed,u,-1,q)

            qrt = dot(rt,q)
            α = ρ/qrt
            axpy!(α,u,x)

            for k=1:N
                πknew[k] = (1+α*(σ[k]-zseed))*πk[k]+ (α*β/αold)*(πkold[k]- πk[k])
                βk = (πkold[k]/πk[k])^2*β
                αk =  (πk[k]/πknew[k])*α
                axpby!(1/πk[k],r,-βk,vec_u[k])
                axpy!(αk,vec_u[k],vec_x[k])
            end

            axpy!(-α,q,r)
            mul!(q,A',ut)
            axpby!(zseed',ut,-1,q)
            axpy!(-α',q,rt)

            #ρMAX = maximum(abs.(πknew))^2
            residual = abs(dot(r,r))#*ρMAX)

            #println(residual,"\t",zseed,"\t")#,dot(r,r),"\t",ρMAX)

            for k=1:N
                πkold[k] = πk[k]
                πk[k] = πknew[k]
            end

            if false

            πkseed = 1
            πkseedold = 1
            kseed = 0
            for k=1:N
                if abs(πk[k]) < abs(πkseed)
                    πkseed =πk[k]
                    πkseedold = πkold[k]
                    zseed = σ[k]
                    kseed = k
                end
            end

            println(kseed)
            if kseed != 0
                axpby!(0,rt,1/πkseed,r)
                axpby!(0,r,1/πkseed',rt)
                #axpby!(0,r,1/πkseedold,u)
                #axpby!(0,r,1/πkseedold',ut)
                axpby!(1,vec_x[kseed],0,x)
                ρ = ρ/(πkseedold'*πkseedold)
                α = (πkseedold/πkseed)*α

                for k=1:N
                    πkold[k] = πkold[k]/πkseedold
                    πk[k] = πk[k]/πkseed
                end
            end

            end

            αold = α
            ρold  = ρ

            if abs(residual) < eps
                println("Converged at $i-th step. eps: $residual")
                println("--------------------------------------")
                return vec_x
            end

        end

        error("bicg is not converged. The residual is $residual")

    end

    function shiftedbicgstab_inSS(σjin,A,b;maxsteps=3000,eps = 1e-15,zseed= 0im) #(sigma_i I + A)*x = b
        σj = deepcopy(σjin) .- zseed
        N = length(σj)
        btype = typeof(b)
        r = deepcopy(b)
        temp1 = similar(b)#temps[2]
        mul!(temp1,A,x)
        axpy!(-1,temp1,r) #add!(r,-1,temp1)

        rs = deepcopy(r)#temps[3]
        #substitute_fermion!(rs,r)
        p = deepcopy(r)#temps[4]
        #substitute_fermion!(p,r)
        Ap = similar(b)#temps[5]
        s = similar(b)#temps[6]
        t = similar(b)#temps[7]

        #r = deepcopy(b)
        #p = deepcopy(b)
        #q = similar(b)
        #s = similar(b)
        x = zero(b)

        vec_r = Vector{btype}(undef,N)
        vec_p = Vector{btype}(undef,N)
        vec_x = Vector{btype}(undef,N)
        for j=1:N
            vec_r[j] = deepcopy(b)
            vec_p[j] = deepcopy(b)
            vec_x[j] = zero(b)
        end

        αm = 1.0
        βm = 0.0

        ρm = ones(ComplexF64,N)
        ρ0 = ones(ComplexF64,N)
        ρp = ones(ComplexF64,N)
        τj = deepcopy(ρ0)
        #τjm = deepcopy(ρ0)
        residual = 0

        mul!(p,A',r) #p = (A + zseed*I)'*r
        axpy!(zseed',r,p)
        c1 = dot(p,p)

        for i=1:maxsteps
            c1 = dot(rs,r)
            mul!(Ap,A,p)
            c2 = dot(rs,Ap)
            α = c1/c2
            #s = r - α*A*p
            #axpby!(1,r,β,p) #add!(β,p,1,r)
            # s = r - α Ap
            axpby!(1,r,0,s) #add!(0,s,1,r)
            axpy!(-α,Ap,s) #add!(s,-α,Ap)
            mul!(t,A,s)
            d1 = dot(t,s)
            d2 = dot(t,t)
            ζk = d1/d2

            #r = (1-ω A)s
            axpby!(1,s,0,r)#add!(0,r,1,s)
            axpy!(-ζk,t,r)  #add!(r,-ω,t)

            #x = x + ωs+ αp
            axpy!(ζk,s,x) #add!(x,ω,s)
            axpy!(α,p,x) #add!(x,α,p)

            β = (dot(rs,r)/c1)*(α/ζk)

            #p = r + β*(1-ωA)*p
            axpby!(1,r,β,p) #add!(β,p,1,r)
            axpy!(-ζk*β,Ap,p) #add!(p,-ω*β,Ap)

            αk = alpha
            βk = beta
            println(dot(r,r))

            for j=1:N
                ρkj = ρ0[j]
                if abs(ρkj) < eps
                    continue
                end
                ρkmj =ρm[j]
                ρp[j] = ρkj*ρkmj*αm/(ρkmj*αm*(1.0+αk*σj[j])+αk*βm*(ρkmj-ρkj))
                αkj = (ρp[j]/ρkj)*αk
                ζkj = ζk/(1+ζk*σj[j])

                #! ...  vec_x[j]   = vec_x[j]   + αkj * vec_p[j] +ζkj*τj*ρp[j]*s
                axpy!(αkj,vec_p[j],vec_x[j])
                axpy!(ζkj*τj[j]*ρp[j],s,vec_x[j])

                τjm = τj[j]
                τj[j] = (1/(1+ζk*σj[j]))*τj[j]
                βkj = (ρp[j]/ρkj)^2*βk
                # vec_p[j] = τj[j]*ρp[j]*r + β*(vec_p[j] + (ζkj/αkj)*τjm*(ρp[j]*s-ρkj*r) )
                # vec_p[j] = β*vec_p[j] + β*(τj[j]*ρp[j]-(ζkj/αkj)*τjm*ρkj)*r + β*(ζkj/αkj)*τjm*ρp[j]*s
                axpby!(βkj*(ζkj/αkj)*τjm*ρp[j],s,βkj,vec_p[j])
                axpy!(βkj*(τj[j]*ρp[j]-(ζkj/αkj)*τjm*ρkj),r,vec_p[j])

                #add!(vec_x[j],αkj,vec_p[j])
                #βkj = (ρp[j]/ρkj)^2*βk
                #add!(βkj,vec_p[j],ρp[j],r)
                axpby!(ρp[j],s,βkj,vec_p[j]) #p = beta*p + s
                #add!(βkj,vec_p[j],ρp[j],s)

            end

            ρm .= ρ0
            ρ0 .= ρp
            αm = αk
            βm = βk

            ρMAX = maximum(abs.(ρp))^2
            residual = abs(rr*ρMAX)

            if abs(residual) < eps
                println("Converged at $i-th step. eps: $residual")
                println("--------------------------------------")
                return x,vec_x
            end

        end
        error("bicg is not converged. The residual is $residual")

    end

    function shiftedbicgstab(x,A,b;eps=1e-10,maxsteps = 1000) #Ax=b
        #println_verbose_level3(verbose,"--------------------------------------")
        #println_verbose_level3(verbose,"bicg-stab method")

        #=
        r = deepcopy(b)
        temp1 = similar(x)
        mul!(temp1,A,x)
        add!(r,-1,temp1)

        rs = deepcopy(r)
        p = deepcopy(r)
        Ap = similar(r)
        s = similar(r)
        t = similar(r)
        =#

        #temps = get_temporaryvectors_forCG(A)
        #r = temps[1]
        #substitute_fermion!(r,b)
        r = deepcopy(b)
        temp1 = similar(b)#temps[2]
        mul!(temp1,A,x)
        axpy!(-1,temp1,r) #add!(r,-1,temp1)

        rs = deepcopy(r)#temps[3]
        #substitute_fermion!(rs,r)
        p = deepcopy(r)#temps[4]
        #substitute_fermion!(p,r)
        Ap = similar(b)#temps[5]
        s = similar(b)#temps[6]
        t = similar(b)#temps[7]

        rnorm = real(dot(r,r))

        if rnorm < eps
            return
        end

        for i=1:maxsteps
            c1 = dot(rs,r)
            mul!(Ap,A,p)
            c2 = dot(rs,Ap)
            α = c1/c2
            #s = r - α*A*p
            #axpby!(1,r,β,p) #add!(β,p,1,r)
            # s = r - α Ap
            axpby!(1,r,0,s) #add!(0,s,1,r)
            axpy!(-α,Ap,s) #add!(s,-α,Ap)
            mul!(t,A,s)
            d1 = dot(t,s)
            d2 = dot(t,t)
            ω = d1/d2

            #r = (1-ω A)s
            axpby!(1,s,0,r)#add!(0,r,1,s)
            axpy!(-ω,t,r)  #add!(r,-ω,t)

            #x = x + ωs+ αp
            axpy!(ω,s,x) #add!(x,ω,s)
            axpy!(α,p,x) #add!(x,α,p)

            β = (dot(rs,r)/c1)*(α/ω)

            #p = r + β*(1-ωA)*p
            axpby!(1,r,β,p) #add!(β,p,1,r)
            axpy!(-ω*β,Ap,p) #add!(p,-ω*β,Ap)

            rnorm = real(r ⋅ r)
            println("$i-th eps: $rnorm")

            if rnorm < eps
                println("Converged at $i-th step. eps: $rnorm")
                #println_verbose_level3(verbose,"--------------------------------------")
                return
            end

        end

        error("""
        The BICGstab is not converged! with maxsteps = $(maxsteps)
        residual is $rnorm
        maxsteps should be larger.""")

    end

    function shiftedbicgstab_inSS(x,A,b;eps=1e-10,maxsteps = 1000) #Ax=b

        r = deepcopy(b)
        temp1 = similar(x)
        mul!(temp1,A,x)
        axpy!(-1,temp1,r)
        #add!(r,-1,temp1)

        rnorm = real(dot(r,r))
        println(rnorm)

        if rnorm < eps
            return
        end
        rs = deepcopy(r)
        p = deepcopy(r)
        s = similar(r)

        s = similar(r)
        Ap = similar(x)
        As = similar(x)
        β = 0.0im
        ω = 0.0im

        for i=0:maxsteps
            if i != 0
                #add!(β,p,1,r)
                #add!(p,-ω*β,Ap)
                axpby!(1,r,β,p) #add!(β,p,1,r)
                axpy!(-ω*β,Ap,p) #add!(p,-ω*β,Ap)
            end

            c1 = dot(rs,r)
            mul!(Ap,A,p)
            c2 = dot(rs,Ap)
            α = c1/c2
            #add!(0,s,1,r)
            axpby!(1,r,0,s)
            #add!(s,-α,Ap)
            axpy!(-α,Ap,s)
            mul!(As,A,s)
            s1 = dot(As,s)
            s2 = dot(As,As)

            ω = s1/s2
            axpy!(α,p,x) #add!(x,α,p)

            axpy!(ω,s,x) #add!(x,ω,s)

            axpby!(1,s,0,r) #add!(0,r,1,s)

            axpy!(-ω,As,r) #add!(r,-ω,As)

            rnorm = real(dot(r,r))
            println("$i-th eps: $rnorm")

            if rnorm < eps
                println("Converged at $i-th step. eps: $rnorm")
                #println_verbose_level3(verbose,"--------------------------------------")
                return
            end

            β = (α/ω)*dot(rs,r)/c1

            #=
            axpby!(1,r,β,p) #add!(β,p,1,r)
            axpy!(-ω*β,Ap,p) #add!(p,-ω*β,Ap)
    =#
        end

        error("""
        The BICG is not converged! with maxsteps = $(maxsteps)
        residual is $rnorm
        maxsteps should be larger.""")

    end

    function shiftedbicg_inSS(σjin,A,b;maxsteps=3000,eps = 1e-15,zseedin = 0) #(sigma_i I + A)*x = b
        zseed = zseedin
        σj = deepcopy(σjin) #.- zseed
        N = length(σj)
        btype = typeof(b)
        r = deepcopy(b)
        p = deepcopy(b)
        q = similar(b)
        s = similar(b)
        x = zero(b)

        vec_r = Vector{btype}(undef,N)
        vec_p = Vector{btype}(undef,N)
        vec_x = Vector{btype}(undef,N)
        for j=1:N
            vec_r[j] = deepcopy(b)
            vec_p[j] = deepcopy(b)
            vec_x[j] = zero(b)
        end

        αm = 1.0
        βm = 0.0

        ρm = ones(ComplexF64,N)
        ρ0 = ones(ComplexF64,N)
        ρp = ones(ComplexF64,N)
        residual = 0

        #=
        mul!(p,A',r) #p = (A + zseed*I)'*r
        axpy!(zseed',r,p)
        =#
        c1 = dot(p,p)

        for i=1:maxsteps
            mul!(q,A,p) #q = (A + zseed*I)*p
            axpy!(zseed,p,q)

            c2 = dot(q,q)
            alpha = c1 / c2
            rr = dot(r,r)
            #! ...  x   = x   + alpha * p
            axpy!(alpha,p,x)
            #axpy!(a, X, Y) #Overwrite Y with X*a + Y, where a is a scalar. Return Y.
            #axpby!(a, X, b, Y) #Overwrite Y with X*a + Y*b

            #...  res = res - alpha * q
            axpy!(-alpha,q,r)

            mul!(s,A',r)#s = (A + zseed*I)'*r
            axpy!(zseed',r,s)

            c3 = dot(s,s)

            beta = c3 / c1
            c1 = c3
            #axpby!(a, X, b, Y) #Y = X*a + Y*b
            axpby!(1,s,beta,p) #p = beta*p + s

            αk = alpha
            βk = beta
            println(dot(r,r))
            ρmin = 1
            jmin = 0
            ρkjmin = 1

            for j=1:N
                ρkj = ρ0[j]
                if abs(ρkj) < eps
                    continue
                end
                ρkmj =ρm[j]
                ρp[j] = ρkj*ρkmj*αm/(ρkmj*αm*(1.0+αk*(σj[j]-zseed))+αk*βm*(ρkmj-ρkj))
                if abs(ρmin) > abs(ρp[j])
                    ρmin = ρp[j]
                    jmin = j
                    ρkjmin = ρkj
                end
                #println(ρp[j])
                αkj = (ρp[j]/ρkj)*αk
                #println("$αk alpha = ",αkj)

                #! ...  vec_x[j]   = vec_x[j]   + αkj * vec_p[j]
                axpy!(αkj,vec_p[j],vec_x[j])
                #println("dot ",dot(x,vec_x[j]))
                #add!(vec_x[j],αkj,vec_p[j])
                βkj = (ρp[j]/ρkj)^2*βk
                #add!(βkj,vec_p[j],ρp[j],r)
                #println("ρp[j] ",ρp[j])
                axpby!(ρp[j],s,βkj,vec_p[j]) #p = beta*p + s
                #add!(βkj,vec_p[j],ρp[j],s)
                #println("dot p ",dot(p,vec_p[j]))

            end

            ρm .= ρ0
            ρ0 .= ρp
            αm = αk
            βm = βk

            ρMAX = maximum(abs.(ρp))^2

            if jmin != 0
                zseed = σj[jmin]
                axpby!(0,s,1/ρmin,p)
                axpby!(1,vec_x[jmin],0,x)
                axpby!(0,p,1/ρmin,s)
                c1 /= ρmin
            end

            residual = abs(rr*ρMAX)

            #println(abs.(ρp))

            if abs(residual) < eps
                println("Converged at $i-th step. eps: $residual")
                println("--------------------------------------")
                return x,vec_x
            end

        end
        error("bicg is not converged. The residual is $residual")

    end

end


text:
using TimeSeries
using Dates
using Random
using PowerSystems

DayAhead  = collect(DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S"):Hour(1):DateTime("1/1/2024  23:00:00", "d/m/y  H:M:S"))
#Dispatch_11am =  collect(DateTime("1/1/2024  0:11:00", "d/m/y  H:M:S"):Minute(15):DateTime("1/1/2024  12::00", "d/m/y  H:M:S"))

nodes5() = [Bus(1,"nodeA", "PV", 0, 1.0, (min = 0.9, max=1.05), 230),
             Bus(2,"nodeB", "PQ", 0, 1.0, (min = 0.9, max=1.05), 230),
             Bus(3,"nodeC", "PV", 0, 1.0, (min = 0.9, max=1.05), 230),
             Bus(4,"nodeD", "REF", 0, 1.0, (min = 0.9, max=1.05), 230),
             Bus(5,"nodeE", "PV", 0, 1.0, (min = 0.9, max=1.05), 230),
            ];

branches5_dc(nodes5) = [Line("1", true, 0.0, 0.0, Arc(from=nodes5[1],to=nodes5[2]), 0.00281, 0.0281, (from=0.00356, to=0.00356), 2.0, (min = -0.7, max = 0.7)),
             HVDCLine("DCL2", true, 0.0, Arc(from=nodes5[1],to=nodes5[4]), (min=-3000.0, max=3000.0), (min=-3000, max=3000), (min=-3000.0, max=3000.0), (min=-3000.0, max=3000.0), (l0=0.0, l1=0.01)),
             Line("3", true, 0.0, 0.0, Arc(from=nodes5[1],to=nodes5[5]), 0.00064, 0.0064, (from=0.01563, to=0.01563), 18.8120, (min = -0.7, max = 0.7)),
             Line("4", true, 0.0, 0.0, Arc(from=nodes5[2],to=nodes5[3]), 0.00108, 0.0108, (from=0.00926, to=0.00926), 11.1480, (min = -0.7, max = 0.7)),
             Line("5", true, 0.0, 0.0, Arc(from=nodes5[3],to=nodes5[4]), 0.00297, 0.0297, (from=0.00337, to=0.00337), 40.530, (min = -0.7, max = 0.7)),
             Line("6", true, 0.0, 0.0, Arc(from=nodes5[4],to=nodes5[5]), 0.00297, 0.0297, (from=0.00337, to=00.00337), 2.00, (min = -0.7, max = 0.7))
];

branches5(nodes5) = [Line("1", true, 0.0, 0.0, Arc(from=nodes5[1],to=nodes5[2]), 0.00281, 0.0281, (from=0.00356, to=0.00356), 2.0, (min = -0.7, max = 0.7)),
             Line("2", true, 0.0, 0.0, Arc(from=nodes5[1],to=nodes5[4]), 0.00304, 0.0304, (from=0.00329, to=0.00329), 2.0, (min = -0.7, max = 0.7)),
             Line("3", true, 0.0, 0.0, Arc(from=nodes5[1],to=nodes5[5]), 0.00064, 0.0064, (from=0.01563, to=0.01563), 18.8120, (min = -0.7, max = 0.7)),
             Line("4", true, 0.0, 0.0, Arc(from=nodes5[2],to=nodes5[3]), 0.00108, 0.0108, (from=0.00926, to=0.00926), 11.1480, (min = -0.7, max = 0.7)),
             Line("5", true, 0.0, 0.0, Arc(from=nodes5[3],to=nodes5[4]), 0.00297, 0.0297, (from=0.00337, to=0.00337), 40.530, (min = -0.7, max = 0.7)),
             Line("6", true, 0.0, 0.0, Arc(from=nodes5[4],to=nodes5[5]), 0.00297, 0.0297, (from=0.00337, to=00.00337), 2.00, (min = -0.7, max = 0.7))
];

branches5_ml(nodes5) = [MonitoredLine("1", true, 0.0, 0.0, Arc(from=nodes5[1],to=nodes5[2]), 0.00281, 0.0281, (from=0.00356, to=0.00356), (from_to=1.0, to_from=1.0), 2.0, (min = -0.7, max = 0.7)),
                Line("2", true, 0.0, 0.0, Arc(from=nodes5[1],to=nodes5[4]), 0.00304, 0.0304, (from=0.00329, to=0.00329), 2.0, (min = -0.7, max = 0.7)),
                Line("3", true, 0.0, 0.0, Arc(from=nodes5[1],to=nodes5[5]), 0.00064, 0.0064, (from=0.01563, to=0.01563), 18.8120, (min = -0.7, max = 0.7)),
                Line("4", true, 0.0, 0.0, Arc(from=nodes5[2],to=nodes5[3]), 0.00108, 0.0108, (from=0.00926, to=0.00926), 11.1480, (min = -0.7, max = 0.7)),
                Line("5", true, 0.0, 0.0, Arc(from=nodes5[3],to=nodes5[4]), 0.00297, 0.0297, (from=0.00337, to=0.00337), 40.530, (min = -0.7, max = 0.7)),
                Line("6", true, 0.0, 0.0, Arc(from=nodes5[4],to=nodes5[5]), 0.00297, 0.0297, (from=0.00337, to=00.00337), 2.00, (min = -0.7, max = 0.7))
                ];

solar_ts_DA = [0
               0
               0
               0
               0
               0
               0
               0
               0
               0.351105684
               0.632536266
               0.99463925
               1
               0.944237283
               0.396681234
               0.366511428
               0.155125829
               0.040872694
               0
               0
               0
               0
               0
               0]

wind_ts_DA = [0.985205412
           0.991791369
           0.997654144
           1
           0.998663733
           0.995497149
           0.992414567
           0.98252418
           0.957203427
           0.927650911
           0.907181989
           0.889095913
           0.848186718
           0.766813846
           0.654052531
           0.525336131
           0.396098004
           0.281771509
           0.197790004
           0.153241012
           0.131355854
           0.113688144
           0.099302656
           0.069569628]

thermal_generators5(nodes5) = [ThermalStandard("Alta", true, nodes5[1], 0.40, 0.010,
           TechThermal(0.5, PowerSystems.ST, PowerSystems.COAL, (min=0.0, max=0.40),  (min = -0.30, max = 0.30), nothing, nothing),
           ThreePartCost((0.0, 1400.0), 0.0, 4.0, 2.0)
           ),
           ThermalStandard("Park City", true, nodes5[1], 1.70, 0.20,
               TechThermal(2.2125, PowerSystems.ST, PowerSystems.COAL, (min=0.0, max=1.70), (min =-1.275, max=1.275), (up=0.02, down=0.02), (up=2.0, down=1.0)),
               ThreePartCost((0.0, 1500.0), 0.0, 1.5, 0.75)
           ),
           ThermalStandard("Solitude", true, nodes5[3], 5.2, 1.00,
               TechThermal(5.20, PowerSystems.ST, PowerSystems.COAL, (min=0.0, max=5.20), (min =-3.90, max=3.90), (up=0.012, down=0.012), (up=3.0, down=2.0)),
               ThreePartCost((0.0, 3000.0), 0.0, 3.0, 1.5)
           ),
           ThermalStandard("Sundance", true, nodes5[4], 2.0, 0.40,
               TechThermal(2.5, PowerSystems.ST, PowerSystems.COAL, (min=0.0, max=2.0), (min =-1.5, max=1.5), (up=0.015, down=0.015), (up=2.0, down=1.0)),
               ThreePartCost((0.0, 4000.0), 0.0, 4.0, 2.0)
           ),
           ThermalStandard("Brighton", true, nodes5[5], 6.0, 1.50,
               TechThermal(7.5, PowerSystems.ST, PowerSystems.COAL, (min=0.0, max=6.0), (min =-4.50, max=4.50), (up=0.015, down=0.015), (up=5.0, down=3.0)),
               ThreePartCost((0.0, 1000.0), 0.0, 1.5, 0.75)
           )];

renewable_generators5(nodes5) = [RenewableDispatch("WindBusA", true, nodes5[5], 0.0, 0.0, PowerSystems.WT, 1.200, TwoPartCost(22.0, 0.0)),
                         RenewableDispatch("WindBusB", true, nodes5[4], 0.0, 0.0, PowerSystems.WT, 1.200, TwoPartCost(22.0, 0.0)),
                         RenewableDispatch("WindBusC", true, nodes5[3], 0.0, 0.0, TechRenewable(1.20, PowerSystems.WT, (min = -0.800, max = 0.800), 1.0), TwoPartCost(22.0, 0.0))];

hydro_generators5(nodes5) = [
                    HydroFix("HydroFix", true, nodes5[2], 0.0, 0.0,
                        TechHydro(0.600, PowerSystems.HY, (min = 0.0, max = 60.0), (min = 0.0, max = 60.0), nothing, nothing)
                    ),
                    HydroDispatch("HydroDispatch", true, nodes5[3], 0.0, 0.0,
                        TechHydro(0.600, PowerSystems.HY, (min = 0.0, max = 60.0), (min = 0.0, max = 60.0), (up = 10.0, down = 10.0), nothing),
                        TwoPartCost(15.0, 0.0))
                    ];

battery5(nodes5) = [GenericBattery(name = "Bat",
                            primemover = PowerSystems.BA,
                            available = true,
                            bus = nodes5[1],
                            energy = 5.0,
                            capacity = (min = 5.0, max = 100.0),
                            rating = 70,
                            activepower = 10.0,
                            inputactivepowerlimits = (min = 0.0, max = 50.0),
                            outputactivepowerlimits = (min = 0.0, max = 50.0),
                            reactivepower = 0.0,
                            reactivepowerlimits = (min = -50.0, max = 50.0),
                            efficiency = (in = 0.80, out = 0.90),
                            )];

loadbus2_ts_DA = [ 0.792729978
                    0.723201574
                    0.710952098
                    0.677672816
                    0.668249175
                    0.67166919
                    0.687608809
                    0.711821241
                    0.756320618
                    0.7984057
                    0.827836527
                    0.840362459
                    0.84511032
                    0.834592803
                    0.822949221
                    0.816941743
                    0.824079963
                    0.905735139
                    0.989967048
                    1
                    0.991227765
                    0.960842114
                    0.921465115
                    0.837001437 ]

loadbus3_ts_DA = [ 0.831093782
                0.689863228
                0.666058513
                0.627033103
                0.624901388
                0.62858924
                0.650734211
                0.683424321
                0.750876413
                0.828347191
                0.884248576
                0.888523615
                0.87752169
                0.847534405
                0.8227661
                0.803809323
                0.813282799
                0.907575962
                0.98679848
                1
                0.990489904
                0.952520972
                0.906611479
                0.824307054]

loadbus4_ts_DA = [ 0.871297342
                0.670489749
                0.642812243
                0.630092987
                0.652991383
                0.671971681
                0.716278493
                0.770885833
                0.810075243
                0.85562361
                0.892440566
                0.910660449
                0.922135467
                0.898416969
                0.879816542
                0.896390855
                0.978598576
                0.96523761
                1
                0.969626503
                0.901212601
                0.81894251
                0.771004923
                0.717847996]

loads5(nodes5) = [ PowerLoad("Bus2", true, nodes5[2], PowerSystems.ConstantPower, 3.0, 0.9861, 3.0, 0.9861),
           PowerLoad("Bus3", true, nodes5[3], PowerSystems.ConstantPower, 3.0, 0.9861, 3.0, 0.9861),
           PowerLoad("Bus4", true, nodes5[4], PowerSystems.ConstantPower, 4.0, 1.3147, 4.0, 1.3147),
        ];

interruptible(nodes5) = [InterruptibleLoad("IloadBus4", true, nodes5[4], PowerSystems.ConstantPower, 0.10, 0.0,  0.10, 0.0, TwoPartCost(150.0, 2400.0))]

#reserve5(thermal_generators5) = [StaticReserve("Reserve1", thermal_generators5, 0.6, maximum([gen.tech.activepowerlimits[:max] for gen in thermal_generators5])),
#StaticReserve("Reserve2", [collect(thermal_generators5)[end]], 0.6, maximum([gen.tech.activepowerlimits[:max] for gen in thermal_generators5]))]

Reserve_ts = [[TimeArray(DayAhead, rand(24))],
              [TimeArray(DayAhead+Day(1), rand(24))]]

hydro_timeseries_DA = [[TimeSeries.TimeArray(DayAhead,wind_ts_DA)],
                     [TimeSeries.TimeArray(DayAhead + Day(1), rand(24)*0.1 + wind_ts_DA)]]

RealTime = collect(DateTime("1/1/2024 0:00:00", "d/m/y H:M:S"):Minute(5):DateTime("1/1/2024 23:55:00", "d/m/y H:M:S"))

load_timeseries_RT = [[TimeArray(RealTime, repeat(loadbus2_ts_DA,inner=12)),
                     TimeArray(RealTime, repeat(loadbus3_ts_DA,inner=12)),
                     TimeArray(RealTime, repeat(loadbus4_ts_DA,inner=12))],
                    [TimeArray(RealTime+Day(1), rand(288)*0.1 + repeat(loadbus2_ts_DA,inner=12)),
                     TimeArray(RealTime+Day(1), rand(288)*0.1 + repeat(loadbus3_ts_DA,inner=12)),
                     TimeArray(RealTime+Day(1), rand(288)*0.1 + repeat(loadbus4_ts_DA,inner=12))]]

load_timeseries_DA = [[TimeSeries.TimeArray(DayAhead,loadbus2_ts_DA),
                    TimeSeries.TimeArray(DayAhead,loadbus3_ts_DA),
                    TimeSeries.TimeArray(DayAhead,loadbus4_ts_DA)],
                   [TimeSeries.TimeArray(DayAhead + Day(1), rand(24)*0.1 + loadbus2_ts_DA),
                    TimeSeries.TimeArray(DayAhead + Day(1), rand(24)*0.1 + loadbus3_ts_DA),
                    TimeSeries.TimeArray(DayAhead + Day(1), rand(24)*0.1 + loadbus4_ts_DA)
                  ]];

ren_timeseries_DA = [[TimeSeries.TimeArray(DayAhead,solar_ts_DA),
                    TimeSeries.TimeArray(DayAhead,wind_ts_DA),
                    TimeSeries.TimeArray(DayAhead,wind_ts_DA)],
                   [TimeSeries.TimeArray(DayAhead + Day(1), rand(24)*0.1 + solar_ts_DA),
                    TimeSeries.TimeArray(DayAhead + Day(1), rand(24)*0.1 + wind_ts_DA),
                    TimeSeries.TimeArray(DayAhead + Day(1), rand(24)*0.1 + wind_ts_DA)
                  ]];

Iload_timeseries_DA = [[TimeArray(DayAhead, loadbus4_ts_DA)],
                      [TimeArray(DayAhead+Day(1), loadbus4_ts_DA + 0.1*rand(24))]]

ext: jl


text:
using Unitful

include("http_parser_benchmark.jl")

using HTTP
using HTTP.IOExtras
using HTTP.ConnectionPool
using HTTP.Streams
using HTTP.Messages

responses = [
"github" => """
HTTP/1.1 200 OK\r
Date: Mon, 15 Jan 2018 00:57:33 GMT\r
Content-Type: text/html; charset=utf-8\r
Transfer-Encoding: chunked\r
Server: GitHub.com\r
Status: 200 OK\r
Cache-Control: no-cache\r
Vary: X-PJAX\r
X-UA-Compatible: IE=Edge,chrome=1\r
Set-Cookie: logged_in=no; domain=.github.com; path=/; expires=Fri, 15 Jan 2038 00:57:33 -0000; secure; HttpOnly\r
Set-Cookie: _gh_sess=XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX; path=/; secure; HttpOnly\r
X-Request-Id: fa15328ce7eb90bfe0274b6c97ed4c59\r
X-Runtime: 0.323267\r
Expect-CT: max-age=2592000, report-uri="https://api.github.com/_private/browser/errors"\r
Content-Security-Policy: default-src 'none'; base-uri 'self'; block-all-mixed-content; child-src render.githubusercontent.com; connect-src 'self' uploads.github.com status.github.com collector.githubapp.com api.github.com www.google-analytics.com github-cloud.s3.amazonaws.com github-production-repository-file-5c1aeb.s3.amazonaws.com github-production-upload-manifest-file-7fdce7.s3.amazonaws.com github-production-user-asset-6210df.s3.amazonaws.com wss://live.github.com; font-src assets-cdn.github.com; form-action 'self' github.com gist.github.com; frame-ancestors 'none'; img-src 'self' data: assets-cdn.github.com identicons.github.com collector.githubapp.com github-cloud.s3.amazonaws.com *.githubusercontent.com; media-src 'none'; script-src assets-cdn.github.com; style-src 'unsafe-inline' assets-cdn.github.com\r
Strict-Transport-Security: max-age=31536000; includeSubdomains; preload\r
X-Content-Type-Options: nosniff\r
X-Frame-Options: deny\r
X-XSS-Protection: 1; mode=block\r
X-Runtime-rack: 0.331298\r
Vary: Accept-Encoding\r
X-GitHub-Request-Id: EB38:2597E:1126615:197B74F:5A5BFC7B\r
\r
"""
,
"wikipedia" => """
HTTP/1.1 200 OK\r
Date: Mon, 15 Jan 2018 01:05:05 GMT\r
Content-Type: text/html\r
Content-Length: 74711\r
Connection: keep-alive\r
Server: mw1327.eqiad.wmnet\r
Cache-Control: s-maxage=86400, must-revalidate, max-age=3600\r
ETag: W/"123d7-56241bf703864"\r
Last-Modified: Mon, 08 Jan 2018 11:03:27 GMT\r
Backend-Timing: D=257 t=1515662958810464\r
Vary: Accept-Encoding\r
X-Varnish: 581763780 554073291, 1859700 2056157, 154619657 149630877, 484912801 88058768\r
Via: 1.1 varnish-v4, 1.1 varnish-v4, 1.1 varnish-v4, 1.1 varnish-v4\r
Age: 56146\r
X-Cache: cp1054 hit/10, cp2016 hit/1, cp4029 hit/2, cp4031 hit/345299\r
X-Cache-Status: hit-front\r
Strict-Transport-Security: max-age=106384710; includeSubDomains; preload\r
Set-Cookie: WMF-Last-Access=15-Jan-2018;Path=/;HttpOnly;secure;Expires=Fri, 16 Feb 2018 00:00:00 GMT\r
Set-Cookie: WMF-Last-Access-Global=15-Jan-2018;Path=/;Domain=.wikipedia.org;HttpOnly;secure;Expires=Fri, 16 Feb 2018 00:00:00 GMT\r
X-Analytics: https=1;nocookies=1\r
X-Client-IP: 14.201.217.194\r
Set-Cookie: GeoIP=AU:VIC:Frankston:-38.14:145.12:v4; Path=/; secure; Domain=.wikipedia.org\r
Accept-Ranges: bytes\r
\r
"""
,
"netflix" => """
HTTP/1.1 200 OK\r
Cache-Control: no-cache, no-store, must-revalidate\r
Content-Type: text/html; charset=utf-8\r
Date: Mon, 15 Jan 2018 01:06:12 GMT\r
Expires: 0\r
Pragma: no-cache\r
req_id: 5ee36317-8bbc-448d-88e8-d19a4fb20331\r
Server: shakti-prod i-0a9188462a5e21164\r
Set-Cookie: SecureNetflixId=v%3D2%26mac%3DAQEAEQABABTLdDaklPBpPR3Gwmj9zmjDfOlSLqOftDk.%26dt%3D1515978372532; Domain=.netflix.com; Path=/; Expires=Fri, 01 Jan 2038 00:00:00 GMT; HttpOnly; Secure\r
Set-Cookie: NetflixId=v%3D2%26ct%3DBQAOAAEBEMuUr7yphqTQpaSgxbmFKYeBoN-TtAcb2E5hDcV73R-Qu8VVqjODTPk7URvqrvzZ1bgpPpiS1DlPOn0ePBInOys9GExY-DUnxkMb99H9ogWoz_Ibi3zWJm7tBsizk3NRHdl4dLZJzrrUGy2HC1JJf3k52FBRx-gFf7UyhseihLocdFS579_IQA1xGugB0pAEadI14eeQkjiZadQ0DwHiZJqjxK8QsUX3_MqMFL9O0q2r5oVMO3sq9VLuxSZz3c_XCFPKTvYSCR_sZPCSPp4kQyIkPUzpQLDTC_dEFapTWGKRNnJBtcH6MJXICcB19SXi83lV26gxFRVyd7MVwZUtGvX_jRGkWgZeRXQuS1YHx0GefQF64y7uEwEgGL49fyKLafF9cHH0tGewQm0iPU3NJktKegMzOx1o4j0-HnWRQzDwgiqQqCLy5sqDCGJyxjagvGdQfc0TiHIVvAeuztEP9XPT1IMvydE8F9C7qzVgpwIaVjkrEzSEG4sazqy7xj5y_dTfOeoHPsQEO2IazjiM1SSMPlDK9SEpiQl18B0sR-tNbYGZgRw6KvSJNBSd9KcsfXsf%26bt%3Ddbl%26ch%3DAQEAEAABABTN0reTusSIN2zATnuy-b6hYllVzhgAaq0.%26mac%3DAQEAEAABABRKcEvZN9gx2qYCWw6GyrTFp5r9efkDmFg.; Domain=.netflix.com; Path=/; Expires=Tue, 15 Jan 2019 06:54:58 GMT; HttpOnly\r
Set-Cookie: clSharedContext=f6f2fc4b-427d-4ac1-b21a-45994d6aa3be; Domain=.netflix.com; Path=/\r
Set-Cookie: memclid=73a77399-8bb5-4a12-bdd8-e07f194d6e3d; Max-Age=31536000; Expires=Tue, 15 Jan 2019 01:06:12 GMT; Path=/; Domain=.netflix.com\r
Set-Cookie: nfvdid=BQFmAAEBEFBHN%2BNNVgwjmGFtH09kBA9AsMDpgxuORCUvQF8ZamZmnBCC1dU4zuJausLoGBA%2FqhGfeoY%2FN0%2F6KN9VlAaHaTU%2BLjtMOr8WsF5dSfyQVjJbNg%3D%3D; Max-Age=26265600; Expires=Thu, 15 Nov 2018 01:06:12 GMT; Path=/; Domain=.netflix.com\r
Strict-Transport-Security: max-age=31536000\r
Vary: Accept-Encoding\r
Via: 1.1 i-0842dec9615b5afe9 (us-west-2)\r
X-Content-Type-Options: nosniff\r
X-Frame-Options: DENY\r
X-Netflix-From-Zuul: true\r
X-Netflix.nfstatus: 1_1\r
X-Netflix.proxy.execution-time: 240\r
X-Originating-URL: https://www.netflix.com/au/\r
X-Xss-Protection: 1; mode=block; report=https://ichnaea.netflix.com/log/freeform/xssreport\r
transfer-encoding: chunked\r
Connection: keep-alive\r
\r
"""
,
"twitter" => """
HTTP/1.1 200 OK\r
cache-control: no-cache, no-store, must-revalidate, pre-check=0, post-check=0\r
content-length: 324899\r
content-type: text/html;charset=utf-8\r
date: Mon, 15 Jan 2018 01:08:31 GMT\r
expires: Tue, 31 Mar 1981 05:00:00 GMT\r
last-modified: Mon, 15 Jan 2018 01:08:31 GMT\r
pragma: no-cache\r
server: tsa_l\r
set-cookie: fm=0; Expires=Mon, 15 Jan 2018 01:08:22 UTC; Path=/; Domain=.twitter.com; Secure; HTTPOnly, _twitter_sess=BAh7CSIKZmxhc2hJQzonQWN0aW9uQ29udHJvbGxlcjo6Rmxhc2g6OkZsYXNo%250ASGFzaHsABjoKQHVzZWR7ADoPY3JlYXRlZF9hdGwrCDxUXPdgAToMY3NyZl9p%250AZCIlOGNhY2U5NjdjOWUyYzkwZWIxMGRiZTgyYjI5NDRkNjY6B2lkIiU4NDdl%250AOTIwNGZhN2JhZjg2NTE5ZjY4ZWJhZTEyOGNjNQ%253D%253D--09d63c3cb46541e33a903d83241480ceab65b33e; Path=/; Domain=.twitter.com; Secure; HTTPOnly, personalization_id="v1_7Q083SuV/G6hK5GjmjV/JQ=="; Expires=Wed, 15 Jan 2020 01:08:31 UTC; Path=/; Domain=.twitter.com, guest_id=v1%3A151597851141844671; Expires=Wed, 15 Jan 2020 01:08:31 UTC; Path=/; Domain=.twitter.com, ct0=d81bcd3211730706c80f4b6cb7b793e8; Expires=Mon, 15 Jan 2018 07:08:31 UTC; Path=/; Domain=.twitter.com; Secure\r
status: 200 OK\r
strict-transport-security: max-age=631138519\r
x-connection-hash: 9cc6e603545c8ee9e1895a9c2acddb89\r
x-content-type-options: nosniff\r
x-frame-options: SAMEORIGIN\r
x-response-time: 531\r
x-transaction: 0019e08a0043f702\r
x-twitter-response-tags: BouncerCompliant\r
x-ua-compatible: IE=edge,chrome=1\r
x-xss-protection: 1; mode=block; report=https://twitter.com/i/xss_report\r
\r
"""
,
"akamai" => """
HTTP/1.1 403 Forbidden\r
Server: AkamaiGHost\r
Mime-Version: 1.0\r
Content-Type: text/html\r
Content-Length: 270\r
Expires: Mon, 15 Jan 2018 01:11:29 GMT\r
Date: Mon, 15 Jan 2018 01:11:29 GMT\r
Connection: close\r
Vary: User-Agent\r
\r
"""
,
"nytimes" => """
HTTP/1.1 200 OK\r
Server: Apache\r
Cache-Control: no-cache\r
X-ESI: 1\r
X-App-Response-Time: 0.62\r
Content-Type: text/html; charset=utf-8\r
X-PageType: homepage\r
X-Age: 69\r
X-Origin-Time: 2018-01-14 20:47:16 EDT\r
Content-Length: 213038\r
Accept-Ranges: bytes\r
Date: Mon, 15 Jan 2018 01:13:41 GMT\r
Age: 1585\r
X-Frame-Options: DENY\r
Set-Cookie: vi_www_hp=z00; path=/; domain=.nytimes.com; expires=Wed, 01 Jan 2020 00:00:00 GMT\r
Set-Cookie: vistory=z00; path=/; domain=.nytimes.com; expires=Wed, 01 Jan 2020 00:00:00 GMT\r
Set-Cookie: nyt-a=4HpgpyMVIOWB1wkVLcX98Q; Expires=Tue, 15 Jan 2019 01:13:41 GMT; Path=/; Domain=.nytimes.com\r
Connection: close\r
X-API-Version: F-5-5\r
Content-Security-Policy: default-src data: 'unsafe-inline' 'unsafe-eval' https:; script-src data: 'unsafe-inline' 'unsafe-eval' https: blob:; style-src data: 'unsafe-inline' https:; img-src data: https: blob:; font-src data: https:; connect-src https: wss:; media-src https: blob:; object-src https:; child-src https: data: blob:; form-action https:; block-all-mixed-content;\r
X-Served-By: cache-mel6524-MEL\r
X-Cache: HIT\r
X-Cache-Hits: 2\r
X-Timer: S1515978821.355774,VS0,VE0\r
Vary: Accept-Encoding, Fastly-SSL\r
\r
"""
]

pipeline_limit = 100

randbody(l) = Vector{UInt8}(rand('A':'Z', l))

const tunit = u"μs"
const tmul = Int(1u"s"/1tunit)
delta_t(a, b) = round((b - a) * tmul, 0)tunit

fields = [:setup, :head, :body, :close, :joyent_head]

function go(count::Int)

    times = Dict()
    for (name, bytes) in responses
        times[name] = Dict()
        for f in fields
            times[name][f] = []
        end
    end

                                                           t_init_start = time()
    io = Base.BufferStream()
    c = Connection("", "", pipeline_limit, 0, true, io)

                                                            t_init_done = time()
    for rep in 1:count
    gc()
    for (name, bytes) in shuffle(responses)

        write(io, bytes)

        r = Request().response
        readheaders(IOBuffer(bytes), r)
        l = bodylength(r)
        if l == unknown_length
            for i = 1:100
                l = 10000
                chunk = randbody(l)
                write(io, hex(l), "\r\n", chunk, "\r\n")
            end
            write(io, "0\r\n\r\n")
        else
            write(io, randbody(l))
        end

                                                                t_start = time()
        t = Transaction(c)
        r = Request()
        s = Stream(r.response, t)
        #startread(s)
        #function IOExtras.startread(http::Stream)

            http = s
            startread(http.stream)

            # Ensure that the header and body are buffered in the Connection
            # object. Otherwise, the time spent in readheaders below is
            # dominated by readavailable() copying huge body data from the
            # Base.BufferStream. We want to measure the parsing performance.
            unread!(http.stream, readavailable(http.stream))

                                                                t_setup = time()
            HTTP.Streams.readheaders(http.stream, http.message)
                                                         t_headers_done = time()
            HTTP.Streams.handle_continue(http)

            http.readchunked = HTTP.Messages.ischunked(http.message)
            http.ntoread = HTTP.Messages.bodylength(http.message)

        #    return http.message
        #end

        while !eof(s)
            readavailable(s)
        end
                                                            t_body_done = time()
        closeread(s)
                                                                 t_done = time()
        push!(times[name][:setup], delta_t(t_start, t_setup))
        push!(times[name][:head],  delta_t(t_setup, t_headers_done))
        push!(times[name][:body],  delta_t(t_headers_done, t_body_done))
        push!(times[name][:close], delta_t(t_body_done, t_done))

                                                                t_start = time()
        r = HttpParserTest.parse(bytes)
                                                                 t_done = time()
        push!(times[name][:joyent_head], delta_t(t_start, t_done))
    end
    end

    if count <= 10
        return
    end

    w = 12
    print("| ")
    print(lpad("Response",w))
    for f in fields
        print(" | ")
        print(lpad(f,w))
    end
    print(" | ")
    print(lpad("jl faster",w))
    println(" |")

    print("| ")
    print(repeat("-",w))
    for f in [fields..., ""]
        print(":| ")
        print(repeat("-",w))
    end
    println(":|")

    for (name, bytes) in responses
        print("| ")
        print(lpad("$name: ",w))
        for f in fields
            print(" | ")
            print(lpad(mean(times[name][f]), w))
        end
        faster = mean(times[name][:joyent_head]) / mean(times[name][:head])
        print(" | ")
        print(lpad("x $(round(faster, 1))", w))
        println(" |")
    end
end

for r in [10, 100]
    go(r)
end


text:
const sf = ScalarReMat(ds[:Batch], ones(size(ds,1)), :Batch, Compat.ASCIIString["(Intercept)"])

const Yield = Array(ds[:Yield])

# check size methods
@test size(sf) == (30, 6)
@test size(sf,1) == 30
@test size(sf,2) == 6
@test size(sf,3) == 1

dd = fill(5., 6)
@test sf'ones(30) == dd

@test ones(30)'sf == dd'
const crp = sf'sf
@test isa(crp, Diagonal{Float64})
const crp1 = copy(crp)
@test crp1 == crp
@test crp[2,6] == 0
@test crp[6,6] == 5
@test size(crp) == (6,6)
@test crp.diag == fill(5.,6)
const rhs = sf'Yield
@test rhs == [7525.0,7640.0,7820.0,7490.0,8000.0,7350.0]
@test A_ldiv_B!(crp,copy(rhs)) == [1505.,1528.,1564.,1498.,1600.,1470.]

const D = Diagonal(ones(30))
const csf = D * sf
@test sf == csf
@test sf == LinAlg.A_mul_B!(csf, D, sf)

@test sf == copy!(csf, sf)

const L = MixedModels.LT(sf)
setθ!(L, [0.5])

@test isa(MixedModels.tscale!(L, crp), Diagonal)
@test crp.diag == fill(2.5, 6)
@test copy!(crp1, crp) == crp

const sf1 = ScalarReMat(psts[:Sample], ones(size(psts,1)), :Sample, Compat.ASCIIString["(Intercept)"])
const sf2 = ScalarReMat(psts[:Batch], ones(size(psts, 1)), :Batch, Compat.ASCIIString["(Intercept)"])
@test size(sf1) == (60, 30)
@test size(sf2) == (60, 10)

const crp11 = sf1'sf1
const pr21 = sf2'sf1
const crp22 = sf2'sf2

@test isa(crp11,Diagonal{Float64})
@test isa(crp22,Diagonal{Float64})
@test isa(pr21,SparseMatrixCSC{Float64})


text:
# Autogenerated wrapper script for xrootdgo_jll for x86_64-apple-darwin
export xrootdgo

JLLWrappers.@generate_wrapper_header("xrootdgo")
JLLWrappers.@declare_library_product(xrootdgo, "@rpath/xrootdgo.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        xrootdgo,
        "lib/xrootdgo.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module SPTestAllGood
using SpatialBoundaries
using Test

@test true

end


text:
__precompile__()

module FredData

using DataFrames
using Dates
using Printf
using TimeZones
import HTTP
import JSON

export
       # Fred object
       Fred, get_api_url, set_api_url!, get_api_key,

       # FredSeries object
       FredSeries,

       # Download data
       get_data

const MAX_ATTEMPTS       = 3
const FIRST_REALTIME     = Date(1776,07,04)
const LAST_REALTIME      = Date(9999,12,31)
const EARLY_VINTAGE_DATE = "1991-01-01"
const FRED_DATE_FORMAT   = DateFormat("yyyy-mm-dd HH:MM:SSzz")
const OUTPUT_TZ_TYPE     = UTC
const DEFAULT_API_URL    = "https://api.stlouisfed.org/fred/"
const API_KEY_LENGTH     = 32
const KEY_ENV_NAME       = "FRED_API_KEY"
const KEY_FILE_NAME      = ".freddatarc"

# Fred connection type
"""
A connection to the Fred API.

Constructors
------------
- `Fred()`: Key detected automatically. First, looks for the environment variable
    `FRED_API_KEY`, then looks for the file `~/.freddatarc`.
- `Fred(key::AbstractString)`: User specifies key directly

Arguments
---------
- `key`: Registration key provided by FRED.

Notes
-----
- Set the API url with `set_api_url!(f::Fred, url::AbstractString)`
"""
mutable struct Fred
    key::AbstractString
    url::AbstractString
    function Fred(key, url)
        # Key validation
        if length(key) > API_KEY_LENGTH
            key = key[1:API_KEY_LENGTH]
            @warn("FRED API key too long. First $(API_KEY_LENGTH) chars used.")
        elseif length(key) < API_KEY_LENGTH
            error("Invalid FRED API key -- key too short: $(key)")
        end
        if !all(isxdigit, key)
            error("Invalid FRED API key -- invalid characters: $(key)")
        end
        return new(key, url)
    end
end

Fred(key::AbstractString) = Fred(key, DEFAULT_API_URL)

key_file() = joinpath(homedir(), KEY_FILE_NAME)

function load_fred_key()
    if KEY_ENV_NAME in keys(ENV)
        ENV[KEY_ENV_NAME]
    elseif isfile(key_file())
        open(key_file(), "r") do file
            rstrip(read(file, String))
        end
    else
        error("FRED API Key not detected.")
    end
end

has_fred_key() = KEY_ENV_NAME in keys(ENV) || isfile(key_file())

function Fred()
    key = load_fred_key()
    println("API key loaded.")
    return Fred(key)
end

"""Get the FRED API key that is used for this connection"""
get_api_key(f::Fred) = f.key

"""Get the base URL used to connect to the FRED server"""
get_api_url(f::Fred) = f.url

"""Set the base URL used to connect to the FRED server"""
set_api_url!(f::Fred, url::AbstractString) = setfield!(f, :url, url)

function Base.show(io::IO, f::Fred)
    @printf io "FRED API Connection\n"
    @printf io "\turl: %s\n" get_api_url(f)
    @printf io "\tkey: %s\n" get_api_key(f)
end

"""
```
FredSeries(...)
```

Represent a single data series, and all associated metadata, as queried from FRED.

The following fields are available:
- `id`: Series ID
- `title`: Series title
- `units_short`: Units (abbr.)
- `units`: Units
- `seas_adj_short`: Seasonal adjustment (abbr.)
- `seas_adj`:Seasonal adjustment
- `freq_short`:*Native* frequency (abbr.)
- `freq`:*Native* frequency
- `realtime_start`:Date realtime period starts
- `realtime_end`:Date realtime period ends
- `last_updated`:Date series last updated
- `notes`:Series notes
- `trans_short`:Transformation of queried data (abbr.)
- `data`:The actual data; DataFrame with columns `:realtime_start`,
  `:realtime_end`, `:date`, `:value`

"""
struct FredSeries
    # From series query
    id::AbstractString
    title::AbstractString
    units_short::AbstractString
    units::AbstractString
    seas_adj_short::AbstractString
    seas_adj::AbstractString
    freq_short::AbstractString
    freq::AbstractString
    realtime_start::AbstractString
    realtime_end::AbstractString
    last_updated::DateTime
    notes::AbstractString

    # From series/observations query
    trans_short::AbstractString # "units"
    data::DataFrames.DataFrame

    # deprecated
    df::DataFrames.DataFrame
end

function Base.show(io::IO, s::FredSeries)
    @printf io "FredSeries\n"
    @printf io "\tid: %s\n"                s.id
    @printf io "\ttitle: %s\n"             s.title
    @printf io "\tunits: %s\n"             s.units
    @printf io "\tseas_adj (native): %s\n" s.seas_adj
    @printf io "\tfreq (native): %s\n"     s.freq
    @printf io "\trealtime_start: %s\n"    s.realtime_start
    @printf io "\trealtime_end: %s\n"      s.realtime_end
    @printf io "\tlast_updated: %s\n"      s.last_updated
    @printf io "\tnotes: %s\n"             s.notes
    @printf io "\ttrans_short: %s\n"       s.trans_short
    @printf io "\tdata: %dx%d DataFrame with columns %s\n" size(s.data)... names(s.data)
end

# old, deprecated accessors
export
    id, title, units_short, units, seas_adj_short, seas_adj, freq_short,
    freq, realtime_start, realtime_end, last_updated, notes, trans_short,
    df
@deprecate id(f::FredSeries) getfield(f, :id)
@deprecate title(f::FredSeries) getfield(f, :title)
@deprecate units_short(f::FredSeries) getfield(f, :units_short)
@deprecate units(f::FredSeries) getfield(f, :units)
@deprecate seas_adj_short(f::FredSeries) getfield(f, :seas_adj_short)
@deprecate seas_adj(f::FredSeries) getfield(f, :seas_adj)
@deprecate freq_short(f::FredSeries) getfield(f, :freq_short)
@deprecate freq(f::FredSeries) getfield(f, :freq)
@deprecate realtime_start(f::FredSeries) getfield(f, :realtime_start)
@deprecate realtime_end(f::FredSeries) getfield(f, :realtime_end)
@deprecate last_updated(f::FredSeries) getfield(f, :last_updated)
@deprecate notes(f::FredSeries) getfield(f, :notes)
@deprecate trans_short(f::FredSeries) getfield(f, :trans_short)
@deprecate df(f::FredSeries) getfield(f, :data)

include("get_data.jl")

end # module


text:
function apply_dissipation_3D!(var_out, var, sys::System)
    Nu, Nx, Ny = size(sys)

    DKOx = sys.DKOx
    DKOy = sys.DKOy

    copyto!(var_out, var)

    @fastmath @inbounds @threads for j in 1:Ny
        @inbounds for i in 1:Nx
            @inbounds for a in 1:Nu
                var_out[a,i,j] += DKOx(var, a,i,j) + DKOy(var, a,i,j)
            end
        end
    end

    nothing
end

function apply_dissipation_2D!(var_out, var, sys::System)
    _, Nx, Ny = size(sys)

    DKOx = sys.DKOx
    DKOy = sys.DKOy

    copyto!(var_out, var)

    @fastmath @inbounds @threads for j in 1:Ny
        @inbounds for i in 1:Nx
            var_out[1,i,j] += DKOx(var, 1,i,j) + DKOy(var, 1,i,j)
        end
    end

    nothing
end

function apply_dissipation!(boundary::Boundary, cache::Boundary, sys::System)
    a4  = geta4(boundary)
    fx2 = getfx2(boundary)
    fy2 = getfy2(boundary)

    a4_cache  = geta4(cache)
    fx2_cache = getfx2(cache)
    fy2_cache = getfy2(cache)

    # the loops in these functions are threaded, so it's probably not worth it
    # to @spawn here
    apply_dissipation_2D!( a4_cache,  a4, sys)
    apply_dissipation_2D!(fx2_cache, fx2, sys)
    apply_dissipation_2D!(fy2_cache, fy2, sys)

    copyto!( a4,  a4_cache)
    copyto!(fx2, fx2_cache)
    copyto!(fy2, fy2_cache)

    nothing
end

function apply_dissipation!(gauge::Gauge, cache::Gauge, sys::System)
    xi = getxi(gauge)
    xi_cache = getxi(cache)

    apply_dissipation_2D!(xi_cache, xi, sys)
    copyto!(xi, xi_cache)

    nothing
end

function apply_dissipation!(bulkevol::BulkEvolved, cache::BulkEvolved,
                            sys::System)
    B1  = getB1(bulkevol)
    B2  = getB2(bulkevol)
    G   = getG(bulkevol)
    phi = getphi(bulkevol)

    B1_cache  = getB1(cache)
    B2_cache  = getB2(cache)
    G_cache   = getG(cache)
    phi_cache = getphi(cache)

    # the loops in these functions are threaded, so it's probably not worth it
    # to @spawn here
    apply_dissipation_3D!( B1_cache,  B1, sys)
    apply_dissipation_3D!( B2_cache,  B2, sys)
    apply_dissipation_3D!(  G_cache,   G, sys)
    apply_dissipation_3D!(phi_cache, phi, sys)

    copyto!( B1,  B1_cache)
    copyto!( B2,  B2_cache)
    copyto!(  G,   G_cache)
    copyto!(phi, phi_cache)

    nothing
end

# exponential filtering
function (filters::Filters)(bulkevol::BulkEvolved)
    @sync begin
        @spawn filters.exp_filter(bulkevol.B1)
        @spawn filters.exp_filter(bulkevol.B2)
        @spawn filters.exp_filter(bulkevol.G)
        @spawn filters.exp_filter(bulkevol.phi)
    end
    nothing
end


text:
################################################################################
#                      TEST: Compose Rotations using DCMs
################################################################################

for i = 1:samples
    # Sample 8 DCMs.
    D1 = angle2dcm(EulerAngles(-pi + 2*pi*rand(),
                               -pi + 2*pi*rand(),
                               -pi + 2*pi*rand(),
                               rand(rot_seq_array)))

    D2 = angle2dcm(EulerAngles(-pi + 2*pi*rand(),
                               -pi + 2*pi*rand(),
                               -pi + 2*pi*rand(),
                               rand(rot_seq_array)))

    D3 = angle2dcm(EulerAngles(-pi + 2*pi*rand(),
                               -pi + 2*pi*rand(),
                               -pi + 2*pi*rand(),
                               rand(rot_seq_array)))

    D4 = angle2dcm(EulerAngles(-pi + 2*pi*rand(),
                               -pi + 2*pi*rand(),
                               -pi + 2*pi*rand(),
                               rand(rot_seq_array)))

    D5 = angle2dcm(EulerAngles(-pi + 2*pi*rand(),
                               -pi + 2*pi*rand(),
                               -pi + 2*pi*rand(),
                               rand(rot_seq_array)))

    D6 = angle2dcm(EulerAngles(-pi + 2*pi*rand(),
                               -pi + 2*pi*rand(),
                               -pi + 2*pi*rand(),
                               rand(rot_seq_array)))

    D7 = angle2dcm(EulerAngles(-pi + 2*pi*rand(),
                               -pi + 2*pi*rand(),
                               -pi + 2*pi*rand(),
                               rand(rot_seq_array)))

    D8 = angle2dcm(EulerAngles(-pi + 2*pi*rand(),
                               -pi + 2*pi*rand(),
                               -pi + 2*pi*rand(),
                               rand(rot_seq_array)))

    # Test the function `compose_rotation`.
    @test                   D2*D1 ≈ compose_rotation(D1,D2)                   atol=1e-14
    @test                D3*D2*D1 ≈ compose_rotation(D1,D2,D3)                atol=1e-14
    @test             D4*D3*D2*D1 ≈ compose_rotation(D1,D2,D3,D4)             atol=1e-14
    @test          D5*D4*D3*D2*D1 ≈ compose_rotation(D1,D2,D3,D4,D5)          atol=1e-14
    @test       D6*D5*D4*D3*D2*D1 ≈ compose_rotation(D1,D2,D3,D4,D5,D6)       atol=1e-14
    @test    D7*D6*D5*D4*D3*D2*D1 ≈ compose_rotation(D1,D2,D3,D4,D5,D6,D7)    atol=1e-14
    @test D8*D7*D6*D5*D4*D3*D2*D1 ≈ compose_rotation(D1,D2,D3,D4,D5,D6,D7,D8) atol=1e-14

end


text:
@testset "Rician" begin

    d1 = Rician(0.0, 10.0)
    @test d1 isa Rician{Float64}
    @test params(d1) == (0.0, 10.0)
    @test shape(d1) == 0.0
    @test scale(d1) == 200.0
    @test partype(d1) === Float64
    @test eltype(d1) === Float64
    @test rand(d1) isa Float64

    d2 = Rayleigh(10.0)
    @test mean(d1) ≈ mean(d2)
    @test var(d1) ≈ var(d2)
    @test mode(d1) ≈ mode(d2)
    @test median(d1) ≈ median(d2)
    @test quantile.(d1, [0.25, 0.45, 0.60, 0.80, 0.90]) ≈ quantile.(d2, [0.25, 0.45, 0.60, 0.80, 0.90])
    @test pdf.(d1, 0.0:0.1:1.0) ≈ pdf.(d2, 0.0:0.1:1.0)
    @test cdf.(d1, 0.0:0.1:1.0) ≈ cdf.(d2, 0.0:0.1:1.0)

    d1 = Rician(10.0, 10.0)
    @test median(d1) == quantile(d1, 0.5)
    x = quantile.(d1, [0.25, 0.45, 0.60, 0.80, 0.90])
    @test all(cdf.(d1, x) .≈ [0.25, 0.45, 0.60, 0.80, 0.90])

    x = rand(Rician(5.0, 5.0), 100000)
    d1 = fit(Rician, x)
    @test d1 isa Rician{Float64}
    @test params(d1)[1] ≈ 5.0 atol=0.2
    @test params(d1)[2] ≈ 5.0 atol=0.2

    d1 = Rician(10.0f0, 10.0f0)
    @test d1 isa Rician{Float32}
    @test params(d1) == (10.0f0, 10.0f0)
    @test shape(d1) == 0.5f0
    @test scale(d1) == 300.0f0
    @test partype(d1) === Float32
    @test eltype(d1) === Float64
    @test rand(d1) isa Float64

    d1 = Rician()
    @test d1 isa Rician{Float64}
    @test params(d1) == (0.0, 1.0)

    @test pdf(d1, -Inf) == 0.0
    @test pdf(d1, -1) == 0.0
    @test pdf(d1, Inf) == 0.0
    @test isnan(pdf(d1, NaN))

    @test logpdf(d1, -Inf) == -Inf
    @test logpdf(d1, -1) == -Inf
    @test logpdf(d1, Inf) == -Inf
    @test isnan(logpdf(d1, NaN))

    @test cdf(d1, -Inf) == 0.0
    @test cdf(d1, -1) == 0.0
    @test cdf(d1, Inf) == 1.0
    @test isnan(cdf(d1, NaN))

    @test logcdf(d1, -Inf) == -Inf
    @test logcdf(d1, -1) == -Inf
    @test logcdf(d1, Inf) == 0.0
    @test isnan(logcdf(d1, NaN))

    @inferred pdf(d1, -Inf32)
    @inferred pdf(d1, -1)
    @inferred pdf(d1, 1.0)
    @inferred pdf(d1, 1.0f0)
    @inferred pdf(d1, 1)
    @inferred pdf(d1, 1//2)
    @inferred pdf(d1, Inf)

    @inferred logpdf(d1, -Inf32)
    @inferred logpdf(d1, -1)
    @inferred logpdf(d1, 1.0)
    @inferred logpdf(d1, 1.0f0)
    @inferred logpdf(d1, 1)
    @inferred logpdf(d1, 1//2)
    @inferred logpdf(d1, Inf)

    @inferred cdf(d1, -Inf32)
    @inferred cdf(d1, -1)
    @inferred cdf(d1, 1.0)
    @inferred cdf(d1, 1.0f0)
    @inferred cdf(d1, 1)
    @inferred cdf(d1, 1//2)
    @inferred cdf(d1, Inf)

    @inferred logcdf(d1, -Inf32)
    @inferred logcdf(d1, -1)
    @inferred logcdf(d1, 1.0)
    @inferred logcdf(d1, 1.0f0)
    @inferred logcdf(d1, 1)
    @inferred logcdf(d1, 1//2)
    @inferred logcdf(d1, Inf)

end


text:
#!/usr/bin/env julia
################################################################################
#
#  Copyright (c) 2015 Wojciech Migda
#  All rights reserved
#  Distributed under the terms of the Apache 2.0 license
#
################################################################################
#
#  Filename: MFCC.jl
#
#  Decription:
#       Mel Frequency Cepstral Coefficients
#
#  Authors:
#       Wojciech Migda
#
################################################################################
#
#  History:
#  --------
#  Date         Who  Ticket     Description
#  ----------   ---  ---------  ------------------------------------------------
#  2015-09-20   wm              Initial version
#
################################################################################

module MFCC

export mfcc

"""
MFCC code ported from https://github.com/jameslyons/python_speech_features
The MIT License (MIT)
"""

function hz2mel(hz::Float64)
    """Convert a value in Hertz to Mels

    :param hz: a value in Hz. This can also be a numpy array, conversion proceeds element-wise.
    :returns: a value in Mels. If an array was passed in, an identical sized array is returned.
    """
    return 2595. * log10(1. + hz / 700.0)
end

function mel2hz(mel::Float64)
    """Convert a value in Mels to Hertz

    :param mel: a value in Mels. This can also be a numpy array, conversion proceeds element-wise.
    :returns: a value in Hertz. If an array was passed in, an identical sized array is returned.
    """
    return 700. * (10. ^ (mel / 2595.0) - 1.)
end

function preemphasis(signal::Vector{Float64}, coeff=0.95)
    """perform preemphasis on the input signal.

    :param signal: The signal to filter.
    :param coeff: The preemphasis coefficient. 0 is no filter, default is 0.95.
    :returns: the filtered signal.
    """
    return vcat(signal[1], signal[2:end] - signal[1:end - 1] * coeff)
end

function preemphasis!(signal::Vector{Float64}, coeff=0.95)
    """perform preemphasis on the input signal.

    :param signal: The signal to filter.
    :param coeff: The preemphasis coefficient. 0 is no filter, default is 0.95.
    """
    #@inbounds for i = -length(signal):-2 signal[-i] -= signal[-i - 1] * coeff end
    #"""
    @inbounds begin
        m = signal[1]
        @simd for i = 2:length(signal)
            n = signal[i] - m * coeff
            m = signal[i]
            signal[i] = n
        end
    end
    #"""
end

function lifter(cepstra::Array{Float64, 2}, L=22.)
    """Apply a cepstral lifter the the matrix of cepstra. This has the effect of increasing the
    magnitude of the high frequency DCT coeffs.

    :param cepstra: the matrix of mel-cepstra, will be numframes * numcep in size.
    :param L: the liftering coefficient to use. Default is 22. L <= 0 disables lifter.
    """
    if L > 0
        ncoeff = size(cepstra, 2)
        lift = 1. + (L / 2.) * sin(pi * [0:ncoeff - 1] / L)
        return cepstra .* lift'
    else
        # values of L <= 0, do nothing
        return cepstra
    end
end

function lifter!(cepstra::Array{Float64, 2}, L=22.)
    """Apply a cepstral lifter the the matrix of cepstra. This has the effect of increasing the
    magnitude of the high frequency DCT coeffs.

    :param cepstra: the matrix of mel-cepstra, will be numframes * numcep in size.
    :param L: the liftering coefficient to use. Default is 22. L <= 0 disables lifter.
    """
    if L > 0
        ncoeff = size(cepstra, 2)
        for j = 1:size(cepstra, 2)
            lift = 1. + (L / 2.) * sin(pi * (j - 1) / L)
            for i = 1:size(cepstra, 1)
                @inbounds cepstra[i, j] *= lift
            end
        end
    end
end

function magspec(frames::Array{Float64, 2}, NFFT)
    """Compute the magnitude spectrum of each frame in frames. If frames is an NxD matrix, output will be NxNFFT.

    :param frames: the array of frames. Each row is a frame.
    :param NFFT: the FFT length to use. If NFFT > frame_len, the frames are zero-padded.
    :returns: If frames is an NxD matrix, output will be NxNFFT. Each row will be the magnitude spectrum of the corresponding frame.
    """
    padded_frames = hcat(frames, repeat(zeros(size(frames, 1)), outer=[1, NFFT - size(frames, 2) % NFFT]))
    complex_spec = rfft(padded_frames, 2)
    #TODO magspec2 -> abs2
    return abs(complex_spec)
end

function magspec2(frames::Array{Float64, 2}, NFFT)
    """Compute the squared magnitude spectrum of each frame in frames. If frames is an NxD matrix, output will be NxNFFT.

    :param frames: the array of frames. Each row is a frame.
    :param NFFT: the FFT length to use. If NFFT > frame_len, the frames are zero-padded.
    :returns: If frames is an NxD matrix, output will be NxNFFT. Each row will be the magnitude spectrum of the corresponding frame.
    """
    #padded_frames = hcat(frames, repeat(zeros(size(frames, 1)), outer=[1, NFFT - size(frames, 2) % NFFT]))
    @time padded_frames = reshape(resize!(vec(frames), size(frames, 1) * NFFT), size(frames, 1), NFFT)

    return abs2(rfft(padded_frames, 2))
end

function powmagspec(frames::Array{Float64, 2}, NFFT)
    """Compute the squared magnitude spectrum of each frame in frames. If frames is an NxD matrix, output will be NxNFFT.

    :param frames: the array of frames. Each row is a frame.
    :param NFFT: the FFT length to use. If NFFT > frame_len, the frames are zero-padded.
    :returns: If frames is an NxD matrix, output will be NxNFFT. Each row will be the magnitude spectrum of the corresponding frame.
    """
    #padded_frames = hcat(frames, repeat(zeros(size(frames, 1)), outer=[1, NFFT - size(frames, 2) % NFFT]))
    padded_vec = resize!(vec(frames), size(frames, 1) * NFFT)
    padded_vec[length(frames) + 1:end] = 0.
    padded_frames = reshape(padded_vec, size(frames, 1), NFFT)

    return abs2(rfft(padded_frames, 2)) / NFFT
end

function magspec2!(spec::Array{Float64, 2}, frames::Array{Float64, 2}, NFFT)
    """Compute the squared magnitude spectrum of each frame in frames. If frames is an NxD matrix, output will be NxNFFT.

    :param frames: the array of frames. Each row is a frame.
    :param NFFT: the FFT length to use. If NFFT > frame_len, the frames are zero-padded.
    :returns: If frames is an NxD matrix, output will be NxNFFT. Each row will be the magnitude spectrum of the corresponding frame.
    """
    padded_frames = hcat(frames, repeat(zeros(size(frames, 1)), outer=[1, NFFT - size(frames, 2) % NFFT]))
    copy!(spec, abs2(rfft(padded_frames, 2)))
end

function powspec(frames::Array{Float64, 2}, NFFT)
    """Compute the power spectrum of each frame in frames. If frames is an NxD matrix, output will be NxNFFT.

    :param frames: the array of frames. Each row is a frame.
    :param NFFT: the FFT length to use. If NFFT > frame_len, the frames are zero-padded.
    :returns: If frames is an NxD matrix, output will be NxNFFT. Each row will be the power spectrum of the corresponding frame.
    """
    return 1.0 / NFFT * magspec(frames, NFFT) .^ 2
end

function powspec2(frames::Array{Float64, 2}, NFFT)
    """Compute the power spectrum of each frame in frames. If frames is an NxD matrix, output will be NxNFFT.

    :param frames: the array of frames. Each row is a frame.
    :param NFFT: the FFT length to use. If NFFT > frame_len, the frames are zero-padded.
    :returns: If frames is an NxD matrix, output will be NxNFFT. Each row will be the power spectrum of the corresponding frame.
    """
    return 1.0 / NFFT * magspec2(frames, NFFT)
end

function powspec2!(pspec::Array{Float64, 2}, frames, NFFT)
    """Compute the power spectrum of each frame in frames. If frames is an NxD matrix, output will be NxNFFT.

    :param frames: the array of frames. Each row is a frame.
    :param NFFT: the FFT length to use. If NFFT > frame_len, the frames are zero-padded.
    :returns: If frames is an NxD matrix, output will be NxNFFT. Each row will be the power spectrum of the corresponding frame.
    """
    magspec2!(pspec, frames, NFFT)
    pspec /= float(NFFT)
end

function framesig(sig::Vector{Float64}, frame_len, frame_step, winfunc=ones)
    """Frame a signal into overlapping frames.

    :param sig: the audio signal to frame.
    :param frame_len: length of each frame measured in samples.
    :param frame_step: number of samples after the start of the previous frame that the next frame should begin.
    :param winfunc: the analysis window to apply to each frame. By default no window is applied.
    :returns: an array of frames. Size is NUMFRAMES by frame_len.
    """
    const slen = size(sig, 1)
    frame_len = int(round(frame_len))
    frame_step = int(round(frame_step))
    if slen <= frame_len
        numframes = 1
    else
        numframes = 1 + int(ceil((1.0 * slen - frame_len) / frame_step))
    end

    const padlen = int((numframes - 1) * frame_step + frame_len)

    # TODO, check which is faster
    #padsignal = vcat(sig, zeros(padlen - slen))
    resize!(sig, padlen)
    sig[padlen + 1:end] = 0

    frames = Array(Float64, numframes, frame_len)
    @simd for q = 1:frame_len
        i = q
        @simd for p = 0:numframes - 1
            @inbounds frames[p + 1, q] = sig[i]
            i += frame_step
        end
    end

    if winfunc != ones
        win = repeat(winfunc(frame_len)', outer=[numframes, 1])
        return frames .* win
    else
        return frames
    end
end

function get_filterbanks(nfilt=20, nfft=512, samplerate=16000, lowfreq=0, highfreq=samplerate / 2)
    """Compute a Mel-filterbank. The filters are stored in the rows, the columns correspond
    to fft bins. The filters are returned as an array of size nfilt * (nfft/2 + 1)

    :param nfilt: the number of filters in the filterbank, default 20.
    :param nfft: the FFT size. Default is 512.
    :param samplerate: the samplerate of the signal we are working with. Affects mel spacing.
    :param lowfreq: lowest band edge of mel filters, default 0 Hz
    :param highfreq: highest band edge of mel filters, default samplerate/2
    :returns: A numpy array of size nfilt * (nfft/2 + 1) containing filterbank. Each row holds 1 filter.
    """
    @assert highfreq <= (samplerate / 2) "highfreq is greater than samplerate/2"

    # compute points evenly spaced in mels
    lowmel = hz2mel(float(lowfreq))
    highmel = hz2mel(float(highfreq))

    melpoints = linspace(lowmel, highmel, nfilt + 2)
    # our points are in Hz, but we use fft bins, so we have to convert
    #  from Hz to fft bin number
    bin = int(floor((nfft + 1) * map(mel2hz, melpoints) / samplerate)) + 1

    fbank = zeros(div(nfft, 2) + 1, nfilt)
    @inbounds (
        @simd for j in 1:nfilt
            @simd for i in int(bin[j]):int(bin[j + 1])
                fbank[i, j] = (i - bin[j]) / (bin[j + 1] - bin[j])
            end
        end
    )
    @inbounds (
        @simd for j in 1:nfilt
            @simd for i in int(bin[j + 1]):int(bin[j + 2])
                fbank[i, j] = (bin[j + 2] - i) / (bin[j + 2] - bin[j + 1])
            end
        end
    )
    return fbank
end

function fbank(signal::Vector{Float64}, samplerate=16000, winlen=0.025, winstep=0.01,
    nfilt=26, nfft=512, lowfreq=0, highfreq=samplerate / 2, preemph=0.97)
    """Compute Mel-filterbank energy features from an audio signal.

    :param signal: the audio signal from which to compute features. Should be an N*1 array
    :param samplerate: the samplerate of the signal we are working with.
    :param winlen: the length of the analysis window in seconds. Default is 0.025s (25 milliseconds)
    :param winstep: the step between successive windows in seconds. Default is 0.01s (10 milliseconds)
    :param nfilt: the number of filters in the filterbank, default 26.
    :param nfft: the FFT size. Default is 512.
    :param lowfreq: lowest band edge of mel filters. In Hz, default is 0.
    :param highfreq: highest band edge of mel filters. In Hz, default is samplerate/2
    :param preemph: apply preemphasis filter with preemph as coefficient. 0 is no filter. Default is 0.97.
    :returns: 2 values. The first is a numpy array of size (NUMFRAMES by nfilt) containing features. Each row holds 1 feature vector. The
        second return value is the energy in each frame (total energy, unwindowed)
    """
    preemphasis!(signal, preemph)

    frames = framesig(signal, winlen * samplerate, winstep * samplerate)
    pspec = powmagspec(frames, nfft)

    energy = sum(pspec, 2) # this stores the total energy in each frame
    @inbounds energy[find(x -> x == 0, energy)] = eps(eltype(energy)) # if energy is zero, we get problems with log

    fb = get_filterbanks(nfilt, nfft, samplerate, lowfreq, highfreq)
    feat = pspec * fb # compute the filterbank energies
    @inbounds feat[find(x -> x == 0, feat)] = eps(eltype(feat)) # if feat is zero, we get problems with log

    return feat, energy
end

################################################################################

function mfcc(signal::Vector{Float64};
              samplerate=16000., winlen=0.025, winstep=0.01, numcep=13,
              nfilt=26, nfft=512, lowfreq=0, highfreq=samplerate/2, preemph=0.97, ceplifter=22, appendEnergy=true)
    """Compute MFCC features from an audio signal.

    :param signal: the audio signal from which to compute features. Should be an N*1 array
    :param samplerate: the samplerate of the signal we are working with.
    :param winlen: the length of the analysis window in seconds. Default is 0.025s (25 milliseconds)
    :param winstep: the step between successive windows in seconds. Default is 0.01s (10 milliseconds)
    :param numcep: the number of cepstrum to return, default 13
    :param nfilt: the number of filters in the filterbank, default 26.
    :param nfft: the FFT size. Default is 512.
    :param lowfreq: lowest band edge of mel filters. In Hz, default is 0.
    :param highfreq: highest band edge of mel filters. In Hz, default is samplerate/2
    :param preemph: apply preemphasis filter with preemph as coefficient. 0 is no filter. Default is 0.97.
    :param ceplifter: apply a lifter to final cepstral coefficients. 0 is no lifter. Default is 22.
    :param appendEnergy: if this is true, the zeroth cepstral coefficient is replaced with the log of the total frame energy.
    :returns: A numpy array of size (NUMFRAMES by numcep) containing features. Each row holds 1 feature vector.
    """
    feat, energy = fbank(signal, samplerate, winlen, winstep, nfilt, nfft, lowfreq, highfreq, preemph)

    feat = log(feat)
    dct!(feat, 2)
    feat = feat[:, 1:numcep]

    lifter!(feat, ceplifter)
    if appendEnergy
        feat[:,1] = log(energy) # replace first cepstral coefficient with log of frame energy
    end

    return feat
end

end # module MFCC


text:
using JLD

jldopen("85.jld") do file
    global X, vocab, vocab2idx, idx2vocab
    X = read(file, "X")
    vocab = read(file, "vocab")
    vocab2idx = read(file, "vocab2idx")
    idx2vocab = read(file, "idx2vocab")
end

println(X[vocab2idx["united_states"],:])


text:
using Documenter
using MLJBase

makedocs(;
    modules=[MLJBase],
    format=Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages=[
        "Home" => "index.md",
        "Measures" => "measures.md",
        "Resampling" => "resampling.md",
        "Composition" => "composition.md",
        "Datasets" => "datasets.md",
        "Distributions" => "distributions.md",
        "OpenML" => "openml.md",
        "Utilities" => "utilities.md"
    ],
    repo="https://github.com/alan-turing-institute/MLJBase.jl/blob/{commit}{path}#L{line}",
    sitename="MLJBase.jl"
)

# By default Documenter does not deploy docs just for PR
# this causes issues with how we're doing things and ends
# up choking the deployment of the docs, so  here we
# force the environment to ignore this so that Documenter
# does indeed deploy the docs
ENV["TRAVIS_PULL_REQUEST"] = "false"

deploydocs(;
    repo="github.com/alan-turing-institute/MLJBase.jl.git",
)


text:
#   This file is part of Maxima.jl. It is licensed under the MIT license
#   Copyright (c) 2016 Nathan Smith

export subst,
       logexpand

import Base: diff,
             float

simfun = [
    :ratsimp,:radcan,
    :factor,:gfactor,:expand,
    :logcontract,
    :makefact,:makegamma,
    :trigsimp,:trigreduce,:trigexpand,:trigrat,
    :rectform,:polarform,
    :realpart,:imagpart,
    :demoivre,
    :exponentialize,
    :float
]

:(export $(simfun...)) |> eval

ty = [:(String),:Expr]

for fun in simfun
    for T in ty
        quote
            function $fun(expr::$T)
                convert($T, $fun(MExpr(expr)))
            end
        end |> eval
    end
    quote
        function $fun(m::MExpr)
            nsr = String[]
            sexpr = split(m).str
            for h in 1:length(sexpr)
                if occursin( ":=", sexpr[h])
                    sp = split(sexpr[h], ":=")
                    push!(nsr,String(sp[1]) * ":=" * string(sp[2] |> String |> MExpr |> $fun))
                elseif occursin( "block([],", sexpr[h])
                    rp = replace(sexpr[h], "block([],", "") |> chop
                    sp = split(rp, ",")
                    ns = "block([],"
                    for u in 1:length(sp)
                        ns = ns * string(sp[u] |> String |> MExpr |> $fun)
                    end
                    ns = ns * ")"
                    push!(nsr, ns)
                elseif occursin( ":", sexpr[h])
                    sp = split(sexpr[h], ":")
                    push!(nsr, String(sp[1]) * ":" * string(sp[2] |> String |> MExpr |> $fun))
                else
                    push!(nsr, $(string(fun)) * "($(sexpr[h]))" |> MExpr |> mcall)
                end
            end
            return MExpr(nsr)
        end
    end |> eval
end

@doc """
    ratsimp(expr::T) where T

Simplify expression.

# Examples

```julia
julia> ratsimp(\"a + b/c\")
\"(a*c+b)/c\"

julia> ratsimp(:(sin(asin(a + b/c))))
:((a * c + b) / c)

julia> ratsimp(m\"%e^log(x)\")

                                       x
```
""" ratsimp

@doc """
    radcan(expr::T) where T

Simplify radicals in expression.

## Examples
```julia
julia> radcan(:(sqrt(x/y)))
:(sqrt(x)/sqrt(y))

julia> radcan(m\"sqrt(x/y)\")

                                    sqrt(x)
                                    -------
                                    sqrt(y)

```
""" radcan

@doc """
    factor(expr::T) where T

Factorize polynomial expression

## Examples

```julia
julia> factor(:(x^2 + 2x + 1))
:((x + 1) ^ 2)

julia> factor(MExpr(\"a^2 - b^2\"))

                                (a - b) (b + a)

```
""" factor

@doc """
    gfactor(expr::T) where T

Factorize complex polynomial expression

## Examples

```julia
julia> gfactor(:(z^2 + 2*im*z - 1))
:((z + im) ^ 2)

```
""" gfactor

"""
    expand(expr)

Expand expression.

## Examples
```julia
julia> expand(m\"(a + b)^2\")

                                 2            2
                                b  + 2 a b + a

```
""" expand

@doc """
    logcontract(expr::T) where T

Contract logarithms in expression
""" logcontract

"""
    logexpand(expr::T) where T

Expand logarithm terms in an expression
""" logexpand
for t in ty
    quote
        function logexpand(expr::$t)
            convert($t, logexpand(MExpr(expr)))
        end
    end |> eval
end

function logexpand(m::MExpr)
    nsr = String[]
    sexpr = split(m).str
    for h in 1:length(sexpr)
        if occursin( ":=", sexpr[h])
            sp = split(sexpr[h], ":=")
            push!(nsr, String(sp[1])*":="*string(sp[2]|>String|>MExpr|>logexpand))
        elseif occursin( "block([],", sexpr[h])
            rp = replace(sexpr[h], "block([],","") |> chop
            sp = split(rp, ",")
            ns = "block([],"
            for u in 1:length(sp)
                ns = ns*string(sp[u] |> String |> MExpr |> logexpand)
            end
            ns = ns * ")"
            push!(nsr,ns)
        elseif occursin(":", sexpr[h])
            sp = split(sexpr[h],":")
            push!(nsr, String(sp[1])*":"*string(sp[2]|>String|>MExpr|>logexpand))
        else
            push!(nsr, "$(sexpr[h]), logexpand=super" |> MExpr |> mcall)
        end
    end
    return MExpr(nsr)
end

@doc """
    makefact(expr::T) where T

Convert expression into factorial form.
""" makefact

@doc """
    makegamma(expr::T) where T

Convert factorial to gamma functions in expression
""" makegamma

@doc """
    trigsimp(expr::T) where T

Simplify trigonometric expression

## Examples
```julia
julia> trigsimp(m\"sin(x)^2 + cos(x)^2\")

                                       1
```
""" trigsimp

@doc """
    trigreduce(expr)

Contract trigonometric functions
""" trigreduce

@doc """
    trigexpand(expr)

Expand out trig functions in expression
""" trigexpand

@doc """
    trigrat(expr)

Convert expression into a canonical trigonometric form

## Examples
```julia
julia> trigrat(:(exp(im*x) + exp(-im*x)))
:(2 * cos(x))

```
""" trigrat

@doc """
    rectform(expr)

Put complex expression in rectangular form

## Examples
```julia
julia> rectform(:(R*e^(im*θ)))
:(R * im * sin(θ) + R * cos(θ))

```
""" rectform

@doc """
    polarform(expr)

Put a complex expression into polarform

## Examples
```julia
julia> polarform(m\"a + %i*b\")

                              2    2    %i atan2(b, a)
                        sqrt(b  + a ) %e

```
""" polarform

@doc """
    realpart(expr)

Find the real part of a complex expression

## Examples
```julia
julia> realpart(MExpr(\"a + %i*b\"))

                            a

```
""" realpart

@doc """
    imagpart(expr)

Find the imaginary part of a complex expression.

## Examples
```julia
julia> realpart(MExpr(\"a + %i*b\"))

                            b

```
""" imagpart

@doc """
    demoivre(expr)

Break exponential terms into hyperbolic and trigonometric functions. Roughly the
opposite in function to `exponentialize`

## Examples
```julia
julia> demoivre(m\"exp(a + %i * b)\")

                             a
                           %e  (%i sin(b) + cos(b))

julia> exponentialize(ans)

                         %i b     - %i b     %i b     - %i b
                    a  %e     + %e         %e     - %e
                  %e  (----------------- + -----------------)
                               2                   2

julia> expand(ans)

                                    %i b + a
                                  %e

```
""" demoivre

@doc """
    exponentialize(expr)

Express expression in terms of exponents as much as possible

## Examples
```julia
julia> exponentialize(m\"sin(x)\")

                                   %i x     - %i x
                             %i (%e     - %e      )
                           - ----------------------
                                       2

```
""" exponentialize

@doc """
    float(expr)

Convert rational numbers into floating point numbers in expression.

## Examples
```julia
julia> float(m\"1/3*x\")

                             0.3333333333333333 x

```
""" float

@doc """
    subst(a, b, expr)

Replace `a` with `b` in `expr`.

## Examples
```julia
julia> subst(:b, :a, :(a^2 + b^2))
:(2 * b ^ 2)

```
""" subst
for t in ty
    quote
        function subst(a, b, expr::$t)
        convert($t, subst(a,b,MExpr(expr)))
        end
    end |> eval
end

function subst(a, b, expr::MExpr)
    str = split(expr).str
    for h in 1:length(str)
        if occursin(":=", str[h])
            sp = split(str[h],":=")
            str[h] = String(str[1])*":="*(_subst(a,b,String(sp[2]))).str
        elseif occursin("block([],", str[h])
            rp = replace(str[h],"block([],","") |> chop
            sp = split(rp,",")
            ns = "block([],"
            for u in 1:length(sp)
                ns = ns * (_subst(a,b,String(sp[u]))).str
            end
            ns = ns * ")"
            str[h] = ns
        elseif occursin( ":", str[h])
            sp = split(str[h], ":")
            str[h] = String(sp[1]) * ":" * (_subst(a,b,String(sp[2]))).str
        else
            str[h] = _subst(a, b, str[h])
        end
    end
    MExpr(str)
end


text:
export Point, Arc, Line, Arcs, Polygon, issimple, iscomplex, boundingbox,
isconcave, isconvex, convexhull

"""
    Point(ϕ::Float64, λ::Float64)

Point type with latitude `ϕ` [deg] and longitude `λ` [deg]
"""
struct Point{T<:Float64}
    ϕ::T
    λ::T
end

(-)(x::Point, y::Point) = Point(x.ϕ - y.ϕ, x.λ - y.λ)
(rad2deg)(x::Point) = Point(rad2deg(x.ϕ), rad2deg(x.λ))
(deg2rad)(x::Point) = Point(deg2rad(x.ϕ), deg2rad(x.λ))

"""
    Arc(point₁::Point, point₂::Point)

Arc type with start position `point₁` [deg] and end position `point₂`
[deg]. The Arc is defined along the shortest great circle distance between
the two points.
"""
struct Arc{T<:Float64}
    point₁::Point{T}
    point₂::Point{T}
end

"""
    Line(point::Point, azimuth::Float64)

Line type with `point` [deg] and `azimuth` [deg]. The line type of a great
circle line that completely circles the sphere.
"""
struct Line{T<:Float64}
    point::Point{T}
    azimuth::Float64
end

"""
    Arcs(points::Vector{Point})

Arcs type with a Vector of points [deg]. A minimum of two points are
necessary for the continuous string of line sections.
"""
struct Arcs{T<:Float64}
    points::Vector{Point{T}}
    function Arcs(points::Vector{Point{T}}) where T<:Float64
        if length(points) > 1
            new{T}(points)
        else
            throw(BoundsError("Invalid Arcs: # Points < 2"))
        end
    end
end

#TODO Ensure that the inside_point is not on any line crossing adjacent points of the polygon
"""
    Polygon(inside_point::Point, points::Vector{Point})

Polygon type with points [deg]. A minimum of three different points are
necessary for a polygon. Also an inside point must be given.
"""
struct Polygon{T<:Float64}
    inside_point::Point{T}
    points::Vector{Point{T}}
    function Polygon(inside_point::Point, points::Vector{Point{T}}) where T<:Float64
        if points[1] != points[end]
            points = vcat(points, points[1])
        end
        if length(points) > 3
            new{T}(inside_point, points)
        else
            throw(BoundsError("Invalid Polygon: # individual Points < 3"))
        end
    end
end

"""
    Polygon(inside_point::Point, arcs::Arcs)

Generate Polygon type based on Arcs. Also an inside point must be given.
"""
Polygon(inside_point::Point, arcs::Arcs) = Polygon(inside_point,
arcs.points)

"""
    issimple(polygon::Polygon)

True if the polygon is 'simple' meaning that it does not self-intersect.
"""
function issimple(polygon::Polygon)
    return !isselfintersecting(polygon)
end

"""
    iscomplex(polygon::Polygon)

True if the polygon is 'complex' meaning that it self-intersects.
"""
iscomplex(polygon::Polygon) = !issimple(polygon)

"""
    isconcave(polygon::Polygon)

Tests if a polygon is concave.

This is true when the polygon is simple and not convex.
"""
isconcave(polygon::Polygon) = issimple(polygon) && !isconvex(polygon)

"""
    isconvex(polygon::Polygon)

Tests if a polygon is convex.

The test is setup by checking exhaustively all possible arcs between all
points of the polygon to intersect with the polygon. If false the polygon
is convex.
"""
function isconvex(polygon::Polygon)
    for (id,pt1) in enumerate(polygon.points[1:end-2])
        for pt2 in polygon.points[id+1:end]
            if length(intersection_points(Arc(pt1,pt2), polygon)) > 2
                return false
            end
        end
    end
    return true
end

"""
    boundingbox(arcs::Arcs)

    Output: [[ϕ_south,ϕ_north],[λ_west,λ_east]]

Generate bounding box values from the given arcs. All arcs stay
within the lat/lon-limits of the bounding box.

Source: Chamberlain & Duquette - Some Algorithms for Polygons on a Sphere
"""
function boundingbox(arcs::Arcs)
    λw = minimum([point.λ for point in arcs.points])
    λe = maximum([point.λ for point in arcs.points])
    pt1 = arcs.points[1]
    ϕmax = -90.0
    ϕmin = 90.0
    for pt2 in arcs.points[2:end]
        ϕmax = max(ϕmax, highest_latitude_point(pt1, pt2).ϕ)
        ϕmin = min(ϕmin, lowest_latitude_point(pt1, pt2).ϕ)
        pt1=pt2
    end
    return [[ϕmin,ϕmax],[λw,λe]]
end

"""
    boundingbox(polygon::Polygon)

    Output: [[ϕ_south,ϕ_north],[λ_west,λ_east]]

Generate bounding box values from the given polygon. This bounding box can be
used for quick point-in-polygon tests. Any arc of the polygon stays within
the lat/lon-limits of the bounding box.

Source: Chamberlain & Duquette - Some Algorithms for Polygons on a Sphere
"""
boundingbox(polygon::Polygon) = boundingbox(Arcs(polygon.points))

"""
    convexhull(polygon::Polygon)

    output: Polygon

Generate a convex hull for the polygon.

Source: stackoverflow.com/questions/9678624/convex-hull-of-longitude-latitude-points-on-the-surface-of-a-sphere
"""
function convexhull(polygon::Polygon)
    n_pts = length(polygon.points) - 1

    inside_point = point_cartesian(polygon.inside_point)
    points = point_cartesian.(polygon.points[1:end-1])
    # move points that are on a great circle arc section to inside_point
    for idA in 1:n_pts
        for idB in 1:n_pts
            if idA == idB
                #
            else
                for id in 1:n_pts
                    if id ≠ idA && id ≠ idB && ison(polygon.points[id], Arc(polygon.points[idA],
                        polygon.points[idB]))
                        points[id] = inside_point
                    end
                end
            end
        end
    end

    # add great circle lines that have all points on or inside the convex hull
    hull_points = []
    for (idA, pointA) in enumerate(points)
        for (idB, pointB) in enumerate(points)
            if pointA == pointB
                #
            else
                great_circle = cross(pointA, pointB)
                tst = map(x -> dot(great_circle, x), points)
                tst[idA] = 0.0
                tst[idB] = 0.0
                if all(tst.≥0.0)
                    append!(hull_points, [[idA, idB]])
                end
            end
        end
    end

    # create convex hull by sequencing points
    convex_set = [hull_points[1][1], hull_points[1][2]]
    for _ in 1:length(hull_points)-1
        id = findfirst(x -> x[1] == convex_set[end], hull_points)
        push!(convex_set, hull_points[id][2])
    end

    return Polygon(polygon.inside_point, polygon.points[convex_set])
end


text:
# TODO: This function should probably be moved to some other file,
#       unless new_package.jl and new_version.jl are merged into one.
function update_status(data::GitHubAutoMergeData; kwargs...)
    if data.read_only
        @info "`read_only` mode; skipping updating the status"
        return nothing
    end
    my_retry(() -> GitHub.create_status(data.api,
                                        data.registry,
                                        data.current_pr_head_commit_sha;
                                        auth = data.auth,
                                        params = Dict(kwargs...)))
end

function pull_request_build(data::GitHubAutoMergeData, ::NewPackage; check_license)::Nothing
    # Rules:
    # 1. Registry consistency tests pass
    # 2.  Only changes a subset of the following files:
    #         - `Registry.toml`,
    #         - `E/Example/Compat.toml`
    #         - `E/Example/Deps.toml`
    #         - `E/Example/Package.toml`
    #         - `E/Example/Versions.toml`
    # 3.  TODO: implement this check. When implemented, this check will make sure that the changes to `Registry.toml` only modify the specified package.
    # 4.  Normal capitalization
    #         - name should match r"^[A-Z]\w*[a-z]\w*[0-9]?$"
    #         - i.e. starts with a capital letter, ASCII alphanumerics only, contains at least 1 lowercase letter
    # 5.  Not too short
    #         - at least five letters
    #         - you can register names shorter than this, but doing so requires someone to approve
    # 6.  Meets julia name check
    #         - does not include the string "julia" with any case
    #         - does not start with "Ju"
    # 7.  DISABLED. Standard initial version number - one of 0.0.1, 0.1.0, 1.0.0, X.0.0
    #         - does not apply to JLL packages
    # 8.  Repo URL ends with /$name.jl.git where name is the package name. We only apply this check if the package is not a subdirectory package.
    # 9.  Compat for Julia
    #         - there should be a [compat] entry for Julia
    #         - the Julia [compat] should have an upper bound
    #         - the Julia [compat] upper bound must have major version 1
    # 10.  Compat for all dependencies
    #         - there should be a [compat] entry for Julia
    #         - all [deps] should also have [compat] entries
    #         - all [compat] entries should have upper bounds
    #         - dependencies that are standard libraries do not need [compat] entries
    #         - dependencies that are JLL packages do not need [compat] entries
    # 11. (only applies to JLL packages) The only dependencies of the package are:
    #         - Pkg
    #         - Libdl
    #         - other JLL packages
    # 12. Package's name has only ASCII characters
    # 13. Version can be installed
    #         - given the proposed changes to the registry, can we resolve and install the new version of the package?
    #         - i.e. can we run `Pkg.add("Foo")`
    # 14. Package code can be downloaded. Can we `git clone` the repo and get the files corresponding to the tree hash in the PR?
    #       - and does that tree hash match what we get by looking at the subdir + commit the registration was triggered from?
    # 15. Package repository contains only OSI-approved license(s) in the license file at toplevel in the version being registered
    # 16. Version can be loaded
    #         - once it's been installed (and built?), can we load the code?
    #         - i.e. can we run `import Foo`
    # 17. Dependency confusion check
    #         - Package UUID doesn't conflict with an UUID in the provided
    #           list of package registries. The exception is if also the
    #           package name *and* package URL matches those in the other
    #           registry, in which case this is a valid co-registration.
    # 18. Package's name is sufficiently far from existing package names in the registry
    #         - We always run this check last.
    #         - We exclude JLL packages from the "existing names"
    #         - We use three checks:
    #             i. that the lowercased name is at least 1 away in Damerau Levenshtein distance from any other lowercased name
    #             ii. that the name is at least 2 away in Damerau Levenshtein distance from any other name
    #             iii. that the name is sufficiently far in a visual distance from any other name
    this_is_jll_package = is_jll_name(data.pkg)
    @info("This is a new package pull request",
          pkg = data.pkg,
          version = data.version,
          this_is_jll_package)

    update_status(data;
                  state = "pending",
                  context = "automerge/decision",
                  description = "New package. Pending.")

    this_pr_can_use_special_jll_exceptions =
        this_is_jll_package && data.authorization == :jll

    # Each element is a tuple of a guideline and whether it's
    # applicable. Instead of a guideline there can be the symbol
    # `:update_status` in which case the PR status will be updated
    # with the results so far before continuing to the following
    # guidelines.
    guidelines =
        [
            (guideline_registry_consistency_tests_pass, true),  # 1
            (guideline_pr_only_changes_allowed_files, true),    # 2
            # (guideline_only_changes_specified_package, true), # 3 (unimplemented)
            (guideline_normal_capitalization,                   # 4
                !this_pr_can_use_special_jll_exceptions),
            (guideline_name_length,                             # 5
                !this_pr_can_use_special_jll_exceptions),
            (guideline_julia_name_check, true),                 # 6
            # (guideline_standard_initial_version_number,       # 7 (deactivated)
            #     !this_pr_can_use_special_jll_exceptions),
            (guideline_repo_url_requirement, true),             # 8
            (guideline_compat_for_julia, true),                 # 9
            (guideline_compat_for_all_deps, true),              # 10
            (guideline_allowed_jll_nonrecursive_dependencies,   # 11
                this_is_jll_package),
            (guideline_name_ascii, true),                       # 12
            (:update_status, true),
            (guideline_version_can_be_pkg_added, true),         # 13
            (guideline_code_can_be_downloaded, true),           # 14
            # `guideline_version_has_osi_license` must be run
            # after `guideline_code_can_be_downloaded` so
            # that it can use the downloaded code!
            (guideline_version_has_osi_license, check_license), # 15
            (guideline_version_can_be_imported, true),          # 16
            (:update_status, true),
            (guideline_dependency_confusion, true),             # 17
            # We always run the `guideline_distance_check`
            # check last, because if the check fails, it
            # prints the list of similar package names in
            # the automerge comment. To make the comment easy
            # to read, we want this list to be at the end.
            (guideline_distance_check, true),                   # 18
         ]

    checked_guidelines = Guideline[]

    for (guideline, applicable) in guidelines
        applicable || continue
        if guideline == :update_status
            if !all(passed, checked_guidelines)
                update_status(data;
                              state = "failure",
                              context = "automerge/decision",
                              description = "New version. Failed.")
            end
        else
            check!(guideline, data)
            @info(guideline.info,
                  meets_this_guideline = passed(guideline),
                  message = message(guideline))
            push!(checked_guidelines, guideline)
        end
    end

    if all(passed, checked_guidelines) # success
        description = "New package. Approved. name=\"$(data.pkg)\". sha=\"$(data.current_pr_head_commit_sha)\""
        update_status(data;
                      state = "success",
                      context = "automerge/decision",
                      description = description)
        this_pr_comment_pass = comment_text_pass(data.registration_type,
                                                 data.suggest_onepointzero,
                                                 data.version,
                                                 this_pr_can_use_special_jll_exceptions)
        my_retry(() -> update_automerge_comment!(data, this_pr_comment_pass))
    else # failure
        update_status(data;
                      state = "failure",
                      context = "automerge/decision",
                      description = "New package. Failed.")
        failing_messages = message.(filter(!passed, checked_guidelines))
        this_pr_comment_fail = comment_text_fail(data.registration_type,
                                                 failing_messages,
                                                 data.suggest_onepointzero,
                                                 data.version)
        my_retry(() -> update_automerge_comment!(data, this_pr_comment_fail))
        throw(AutoMergeGuidelinesNotMet("The automerge guidelines were not met."))
    end
    return nothing
end