text:
################################################################################
#  Copyright 2022, Tom Van Acker (BASF)                                        #
################################################################################
# TFT.jl                                                                       #
# A Julia package for Taylor-Fourier Transform.                                #
# See http://github.com/timmyfaraday/TFT.jl                                    #
################################################################################

@testset "FTFT" begin

    # fundamental frequency and angular frequency
    F   = 50.0
    ω   = 2 * pi * F

    # tft input
    D   = 2
    K   = 9

    # discrete time
    t   = 0.0:0.001:1.0
    tm  = 0.536
    idm = findfirst(x -> x == tm, t)

    @testset "Zeroth Harmonic Constant Signal" begin
        # input - amplitude
        A(t)    = 10.0

        # derivatives - amplitude and anti-rotating phase
        dA(t)   = _FD.derivative(A,t)
        d2A(t)  = _FD.derivative(dA,t)

        # derived input
        Ξ(t)    = A.(t)
        Ψ(t)    = A.(t)
        S(t)    = A.(t)

        # perform taylor-fourier transform
        sol = ftft(S.(t), collect(t), [0,10], D, F, K)

        # tests
        ## amplitude
        @test isapprox(TFT.a(sol,0,0)[idm], A(tm), atol=atol)
        @test isapprox(TFT.a(sol,1,0)[idm], dA(tm), atol=atol)
        @test isapprox(TFT.a(sol,2,0)[idm], d2A(tm), atol=atol)
        ### --- amplitude of non-present harmonic should be zero
        @test isapprox(TFT.a(sol,0,10)[idm], 0.0, atol=atol)
        ## phase
        @test isapprox(TFT.ϕ(sol,0,0)[idm], 0.0, atol=atol)
        ## anti-rotating phase
        @test isapprox(TFT.φ(sol,0,0)[idm], 0.0, atol=atol)
        ## dynamic phasor
        @test isapprox(TFT.ξ(sol,0,0)[idm], Ξ(tm), atol=atol)
        ## anti-rotating dynamic phasor
        @test isapprox(TFT.ψ(sol,0,0)[idm], Ψ(tm), atol=atol)
        @test isapprox(TFT.ψ(sol,1,0)[idm], TFT.ξ(sol,1,0)[idm], atol=atol)
        @test isapprox(TFT.ψ(sol,2,0)[idm], TFT.ξ(sol,2,0)[idm], atol=atol)
        ## signal
        @test isapprox(TFT.signal(sol)[idm], S(tm), atol=atol)
        @test isapprox(TFT.signal(sol,0)[idm], S(tm), atol=atol)
        ### --- signal of non-present harmonic should be zero
        @test isapprox(TFT.signal(sol,10)[idm], 0.0, atol=atol)
        ## error
        @test isapprox(TFT.error(sol)[idm], 0.0, atol=atol)

    end

    @testset "Zeroth Harmonic Linear Signal" begin
        # input - amplitude
        A(t)    = 10.0 .- t

        # derivatives - amplitude and anti-rotating phase
        dA(t)   = _FD.derivative(A,t)
        d2A(t)  = _FD.derivative(dA,t)

        # derived input
        Ξ(t)    = A.(t)
        Ψ(t)    = A.(t)
        S(t)    = A.(t)

        # perform taylor-fourier transform
        sol = ftft(S.(t), collect(t), [0,10], D, F, K)

        # tests
        ## amplitude
        @test isapprox(TFT.a(sol,0,0)[idm], A(tm), atol=atol)
        @test isapprox(TFT.a(sol,1,0)[idm], dA(tm), atol=atol)
        @test isapprox(TFT.a(sol,2,0)[idm], d2A(tm), atol=atol)
        ### --- amplitude of non-present harmonic should be zero
        @test isapprox(TFT.a(sol,0,10)[idm], 0.0, atol=atol)
        ## phase
        @test isapprox(TFT.ϕ(sol,0,0)[idm], 0.0, atol=atol)
        ## anti-rotating phase
        @test isapprox(TFT.φ(sol,0,0)[idm], 0.0, atol=atol)
        ## dynamic phasor
        @test isapprox(TFT.ξ(sol,0,0)[idm], Ξ(tm), atol=atol)
        ## anti-rotating dynamic phasor
        @test isapprox(TFT.ψ(sol,0,0)[idm], Ψ(tm), atol=atol)
        @test isapprox(TFT.ψ(sol,1,0)[idm], TFT.ξ(sol,1,0)[idm], atol=atol)
        @test isapprox(TFT.ψ(sol,2,0)[idm], TFT.ξ(sol,2,0)[idm], atol=atol)
        ## signal
        @test isapprox(TFT.signal(sol)[idm], S(tm), atol=atol)
        @test isapprox(TFT.signal(sol,0)[idm], S(tm), atol=atol)
        ### --- signal of non-present harmonic should be zero
        @test isapprox(TFT.signal(sol,10)[idm], 0.0, atol=atol)
        ## error
        @test isapprox(TFT.error(sol)[idm], 0.0, atol=atol)

    end

    @testset "Zeroth Harmonic Quadratic Signal" begin
        # input - amplitude
        A(t)    = 10.0 .- t .+ 2.0 .* t.^2

        # derivatives - amplitude and anti-rotating phase
        dA(t)   = _FD.derivative(A,t)
        d2A(t)  = _FD.derivative(dA,t)

        # derived input
        Ξ(t)    = A.(t)
        Ψ(t)    = A.(t)
        S(t)    = A.(t)

        # perform taylor-fourier transform
        sol = ftft(S.(t), collect(t), [0,10], D, F, K)

        # tests
        ## amplitude
        @test isapprox(TFT.a(sol,0,0)[idm], A(tm), atol=atol)
        @test isapprox(TFT.a(sol,1,0)[idm], dA(tm), atol=atol)
        @test isapprox(TFT.a(sol,2,0)[idm], d2A(tm), atol=atol)
        ### --- amplitude of non-present harmonic should be zero
        @test isapprox(TFT.a(sol,0,10)[idm], 0.0, atol=atol)
        ## phase
        @test isapprox(TFT.ϕ(sol,0,0)[idm], 0.0, atol=atol)
        ## anti-rotating phase
        @test isapprox(TFT.φ(sol,0,0)[idm], 0.0, atol=atol)
        ## dynamic phasor
        @test isapprox(TFT.ξ(sol,0,0)[idm], Ξ(tm), atol=atol)
        ## anti-rotating dynamic phasor
        @test isapprox(TFT.ψ(sol,0,0)[idm], Ψ(tm), atol=atol)
        @test isapprox(TFT.ψ(sol,1,0)[idm], TFT.ξ(sol,1,0)[idm], atol=atol)
        @test isapprox(TFT.ψ(sol,2,0)[idm], TFT.ξ(sol,2,0)[idm], atol=atol)
        ## signal
        @test isapprox(TFT.signal(sol)[idm], S(tm), atol=atol)
        @test isapprox(TFT.signal(sol,0)[idm], S(tm), atol=atol)
        ### --- signal of non-present harmonic should be zero
        @test isapprox(TFT.signal(sol,10)[idm], 0.0, atol=atol)
        ## error
        @test isapprox(TFT.error(sol)[idm], 0.0, atol=atol)

    end

    @testset "Fundamental Periodic Signal" begin
        # input - amplitude and anti-rotating phase
        A(t)    = 10.0
        Φ(t)    = 2 * pi / 3

        # derivatives - amplitude and anti-rotating phase
        dA(t)   = _FD.derivative(A,t)
        d2A(t)  = _FD.derivative(dA,t)
        dΦ(t)   = _FD.derivative(Φ,t)
        d2Φ(t)  = _FD.derivative(dΦ,t)

        # derived input
        Φr(t)   = rem.(ω .*t .+ Φ.(t) .+ pi, 2 * pi) .- pi
        Fr(t)   = F + dΦ(t) / (2 * pi)
        Rr(t)   = d2Φ(t) / (2 * pi)^2
        Ξ(t)    = A.(t) ./ 2 .* exp.(im .* Φ.(t)) .* exp.(im .* ω .* t)
        Ψ(t)    = A.(t) ./ 2 .* exp.(im .* Φ.(t))
        S(t)    = real.(A.(t) ./ 2 .* exp.(im .* Φ.(t)) .* exp.(im .* ω .* t) .+
                  conj.(A.(t) ./ 2 .* exp.(im .* Φ.(t)) .* exp.(im .* ω .* t)))

        # perform taylor-fourier transform
        sol = ftft(S.(t), collect(t), [1,10], D, F, K)

        # tests
        ## amplitude
        @test isapprox(TFT.a(sol,0,1)[idm], A(tm), atol=atol)
        @test isapprox(TFT.a(sol,1,1)[idm], dA(tm), atol=atol)
        @test isapprox(TFT.a(sol,2,1)[idm], d2A(tm), atol=atol)
        ### --- amplitude of derivative degree higher than two not supported
        @test TFT.a(sol,3,1) === nothing
        ### --- amplitude of non-present harmonic should be zero
        @test isapprox(TFT.a(sol,0,10)[idm], 0.0, atol=atol)
        ## phase
        @test isapprox(TFT.ϕ(sol,0,1)[idm], Φr(tm), atol=atol)
        @test isapprox(TFT.ϕ(sol,1,1)[idm], dΦ(tm), atol=atol)
        @test isapprox(TFT.ϕ(sol,2,1)[idm], d2Φ(tm), atol=atol)
        ### --- phase of derivative degree higher than two not supported
        @test TFT.ϕ(sol,3,1) === nothing
        ## anti-rotating phase
        @test isapprox(TFT.φ(sol,0,1)[idm], Φ(tm), atol=atol)
        @test isapprox(TFT.φ(sol,1,1)[idm], dΦ(tm), atol=atol)
        @test isapprox(TFT.φ(sol,2,1)[idm], d2Φ(tm), atol=atol)
        ### --- ar phase of derivative degree higher than two not supported
        @test TFT.φ(sol,3,1) === nothing
        ## frequency
        @test isapprox(TFT.f(sol,1)[idm], Fr(tm), atol=atol)
        ## rocof
        @test isapprox(TFT.r(sol,1)[idm], Rr(tm), atol=atol)
        ## dynamic phasor
        @test isapprox(TFT.ξ(sol,0,1)[idm], Ξ(tm), atol=atol)
        ## anti-rotating dynamic phasor
        @test isapprox(TFT.ψ(sol,0,1)[idm], Ψ(tm), atol=atol)
        @test isapprox(TFT.ψ(sol,1,1)[idm], TFT.ξ(sol,1,1)[idm] * exp(-im * ω * tm), atol=atol)
        @test isapprox(TFT.ψ(sol,2,1)[idm], TFT.ξ(sol,2,1)[idm] * exp(-im * ω * tm), atol=atol)
        ## signal
        @test isapprox(TFT.signal(sol)[idm], S(tm), atol=atol)
        @test isapprox(TFT.signal(sol,1)[idm], S(tm), atol=atol)
        ### --- signal of non-present harmonic should be zero
        @test isapprox(TFT.signal(sol,10)[idm], 0.0, atol=atol)
        ## error
        @test isapprox(TFT.error(sol)[idm], 0.0, atol=atol)
    end

    @testset "Fundamental Aperiodic Linear Signal" begin
        # input - amplitude and anti-rotating phase
        A(t)    = 10.0 .- t
        Φ(t)    = pi / 2 .* t

        # derivatives - amplitude and anti-rotating phase
        dA(t)   = _FD.derivative(A,t)
        d2A(t)  = _FD.derivative(dA,t)
        dΦ(t)   = _FD.derivative(Φ,t)
        d2Φ(t)  = _FD.derivative(dΦ,t)

        # derived input
        Φr(t)   = rem.(ω .*t .+ Φ.(t) .+ pi, 2 * pi) .- pi
        Fr(t)   = F + dΦ(t) / (2 * pi)
        Rr(t)   = d2Φ(t) / (2 * pi)^2
        Ξ(t)    = A.(t) ./ 2 .* exp.(im .* Φ.(t)) .* exp.(im .* ω .* t)
        Ψ(t)    = A.(t) ./ 2 .* exp.(im .* Φ.(t))
        S(t)    = real.(A.(t) ./ 2 .* exp.(im .* Φ.(t)) .* exp.(im .* ω .* t) .+
                  conj.(A.(t) ./ 2 .* exp.(im .* Φ.(t)) .* exp.(im .* ω .* t)))

        # perform taylor-fourier transform
        sol = ftft(S.(t), collect(t), [1,10], D, F, K)

        # tests
        ## amplitude
        @test isapprox(TFT.a(sol,0,1)[idm], A(tm), atol=atol)
        @test isapprox(TFT.a(sol,1,1)[idm], dA(tm), atol=atol)
        @test isapprox(TFT.a(sol,2,1)[idm], d2A(tm), atol=atol)
        ### --- amplitude of derivative degree higher than two not supported
        @test TFT.a(sol,3,1) === nothing
        ### --- amplitude of non-present harmonic should be zero
        @test isapprox(TFT.a(sol,0,10)[idm], 0.0, atol=atol)
        ## phase
        @test isapprox(TFT.ϕ(sol,0,1)[idm], Φr(tm), atol=atol)
        @test isapprox(TFT.ϕ(sol,1,1)[idm], dΦ(tm), atol=atol)
        @test isapprox(TFT.ϕ(sol,2,1)[idm], d2Φ(tm), atol=atol)
        ### --- phase of derivative degree higher than two not supported
        @test TFT.ϕ(sol,3,1) === nothing
        ## anti-rotating phase
        @test isapprox(TFT.φ(sol,0,1)[idm], Φ(tm), atol=atol)
        @test isapprox(TFT.φ(sol,1,1)[idm], dΦ(tm), atol=atol)
        @test isapprox(TFT.φ(sol,2,1)[idm], d2Φ(tm), atol=atol)
        ### --- ar phase of derivative degree higher than two not supported
        @test TFT.φ(sol,3,1) === nothing
        ## frequency
        @test isapprox(TFT.f(sol,1)[idm], Fr(tm), atol=atol)
        ## rocof
        @test isapprox(TFT.r(sol,1)[idm], Rr(tm), atol=atol)
        ## dynamic phasor
        @test isapprox(TFT.ξ(sol,0,1)[idm], Ξ(tm), atol=atol)
        ## anti-rotating dynamic phasor
        @test isapprox(TFT.ψ(sol,0,1)[idm], Ψ(tm), atol=atol)
        @test isapprox(TFT.ψ(sol,1,1)[idm], TFT.ξ(sol,1,1)[idm] * exp(-im * ω * tm), atol=atol)
        @test isapprox(TFT.ψ(sol,2,1)[idm], TFT.ξ(sol,2,1)[idm] * exp(-im * ω * tm), atol=atol)
        ## signal
        @test isapprox(TFT.signal(sol)[idm], S(tm), atol=atol)
        @test isapprox(TFT.signal(sol,1)[idm], S(tm), atol=atol)
        ### --- signal of non-present harmonic should be zero
        @test isapprox(TFT.signal(sol,10)[idm], 0.0, atol=atol)
        ## error
        @test isapprox(TFT.error(sol)[idm], 0.0, atol=atol)
    end

    @testset "Fundamental Aperiodic Quadratic Signal" begin
        # input - amplitude and anti-rotating phase
        A(t)    = 10.0 .- t .+ 2.0 .* t.^2
        Φ(t)    = pi / 2 .* t.^2

        # derivatives - amplitude and anti-rotating phase
        dA(t)   = _FD.derivative(A,t)
        d2A(t)  = _FD.derivative(dA,t)
        dΦ(t)   = _FD.derivative(Φ,t)
        d2Φ(t)  = _FD.derivative(dΦ,t)

        # derived input
        Φr(t)   = rem.(ω .* t .+ Φ.(t) .+ pi, 2 * pi) .- pi
        Fr(t)   = F + dΦ(t) / (2 * pi)
        Rr(t)   = d2Φ(t) / (2 * pi)^2
        Ξ(t)    = A.(t) ./ 2 .* exp.(im .* Φ.(t)) .* exp.(im .* ω .* t)
        Ψ(t)    = A.(t) ./ 2 .* exp.(im .* Φ.(t))
        S(t)    = real.(A.(t) ./ 2 .* exp.(im .* Φ.(t)) .* exp.(im .* ω .* t) .+
                  conj.(A.(t) ./ 2 .* exp.(im .* Φ.(t)) .* exp.(im .* ω .* t)))

        # perform taylor-fourier transform -- not passing specific harmonics
        # to use the alternative problem_builder function.
        sol = ftft(S.(t), collect(t), D, F, K)

        # tests
        ## amplitude
        @test isapprox(TFT.a(sol,0,1)[idm], A(tm), atol=atol)
        @test isapprox(TFT.a(sol,1,1)[idm], dA(tm), atol=atol)
        @test isapprox(TFT.a(sol,2,1)[idm], d2A(tm), atol=atol)
        ### --- amplitude of derivative degree higher than two not supported
        @test TFT.a(sol,3,1) === nothing
        ### --- amplitude of non-present harmonic should be zero
        @test isapprox(TFT.a(sol,0,10)[idm], 0.0, atol=atol)
        ## phase
        @test isapprox(TFT.ϕ(sol,0,1)[idm], Φr(tm), atol=atol)
        @test isapprox(TFT.ϕ(sol,1,1)[idm], dΦ(tm), atol=atol)
        @test isapprox(TFT.ϕ(sol,2,1)[idm], d2Φ(tm), atol=atol)
        ### --- amplitude of derivative degree higher than two not supported
        @test TFT.ϕ(sol,3,1) === nothing
        ## anti-rotating phase
        @test isapprox(TFT.φ(sol,0,1)[idm], Φ(tm), atol=atol)
        @test isapprox(TFT.φ(sol,1,1)[idm], dΦ(tm), atol=atol)
        @test isapprox(TFT.φ(sol,2,1)[idm], d2Φ(tm), atol=atol)
        ### --- ar phase of derivative degree higher than two not supported
        @test TFT.φ(sol,3,1) === nothing
        ## frequency
        @test isapprox(TFT.f(sol,1)[idm], Fr(tm), atol=atol)
        ## rocof
        @test isapprox(TFT.r(sol,1)[idm], Rr(tm), atol=atol)
        ## dynamic phasor
        @test isapprox(TFT.ξ(sol,0,1)[idm], Ξ(tm), atol=atol)
        ## anti-rotating dynamic phasor
        @test isapprox(TFT.ψ(sol,0,1)[idm], Ψ(tm), atol=atol)
        @test isapprox(TFT.ψ(sol,1,1)[idm], TFT.ξ(sol,1,1)[idm] * exp(-im * ω * tm), atol=atol)
        @test isapprox(TFT.ψ(sol,2,1)[idm], TFT.ξ(sol,2,1)[idm] * exp(-im * ω * tm), atol=atol)
        ## signal
        @test isapprox(TFT.signal(sol)[idm], S(tm), atol=atol)
        @test isapprox(TFT.signal(sol,1)[idm], S(tm), atol=atol)
        ### --- signal of non-present harmonic should be zero
        @test isapprox(TFT.signal(sol,10)[idm], 0.0, atol=atol)
        ## error
        @test isapprox(TFT.error(sol)[idm], 0.0, atol=atol)
    end

end


text:
function coef_H_uniform!(r::Integer, N::Integer, M::Integer,
                         store::Store{T}, version=2) where T <: AbstractFloat
    H = OffsetVector{Matrix{T}}(undef, 0:N-1)
    for ℓbar = 0:N-1
        H[ℓbar] = Array{T}(undef, r, r)
    end
    coef_H_uniform!(H, M, store, version)
    return H
end

function coef_H_uniform!(H::OffsetVector{Matrix{T}}, M::Integer,
                         store::Store{T}, version=2) where T <: AbstractFloat
    N = length(H)
    coef_H0_uniform!(H[0], store)
    if N ≥ 2
        coef_H1_uniform!(H[1], M, store)
    end
    if N ≥ 3
        if version == 1
            coef_H_uniform_ver1!(H, 2:N-1, M, store)
        elseif version == 2
            coef_H_uniform_ver2!(H, 2:N-1, M, store)
        else
            throw(ArgumentError("version must be 1 or 2"))
        end
    end
end

# H^0
function coef_H0_uniform!(r::Integer, store::Store{T}) where T <: AbstractFloat
    H0 =Array{T}(undef, r, r)
    coef_H0_uniform!(H0, store)
    return H0
end

function coef_H0_uniform!(H0::Matrix{T},
                          store::Store{T}) where T <: AbstractFloat
    r = size(H0, 1)
    α = store.α
    @argcheck r ≤ store.rmax
    @argcheck size(H0, 2) == r
    y, wy = rules(store.jacobi5)
    z, wz = rules(store.legendre)
    Φ  = view(store.C, 1:2r)
    for j = 1:r, i = 1:r
        Ψ  = view(store.Ψ, 1:j)
        dΨ = view(store.dΨ, 1:i)
        My = ceil(Integer, (i+j-1)/2)
        Mz = ceil(Integer, (i+j)/2)
        for my = 1:My
            Φ[my] = zero(T)
            y_ = y[My][my]
            for mz = 1:Mz
                z_ = z[Mz][mz]
                τ = ( y_ - z_ * y_ + 1 + z_ ) / 2
                σ = τ - 1 - y_
                legendre_polys!(Ψ, σ)
                deriv_legendre_polys!(dΨ, τ)
                Φ[my] += wz[Mz][mz] * Ψ[j] * dΨ[i]
            end
            Φ[my] /= 2
        end
        s = zero(T)
        for my = 1:My
            s += wy[My][my] * Φ[my]
        end
        H0[i,j] = s
    end
    σ, w = rules(store.jacobi1)
    c = 1 / ( 2^α * Γ(α) )
    for j = 1:r
        Ψ = view(store.Ψ, 1:j)
        M = ceil(Integer, j/2)
        s = zero(T)
        for m = 1:M
            σ_ = σ[M][m]
            legendre_polys!(Ψ, σ_)
            s += w[M][m] * Ψ[j]
        end
        for i = 1:r
            H0[i,j] = c * ( s - H0[i,j] )
        end
    end
end

# H^1
function coef_H1_uniform!(r::Integer, M::Integer,
                          store::Store{T}) where T <: AbstractFloat
    H1 = Array{T}(undef, r, r)
    coef_H1_uniform!(H1, M, store)
    return H1
end

function coef_H1_uniform!(H1::Matrix{T}, M::Integer,
                          store::Store{T}) where T <: AbstractFloat
    r = size(H1, 1)
    α = store.α
    @argcheck size(H1, 2) == r
    @argcheck r ≤ store.rmax

    σ, wσ = rule(store.legendre[M])
    Ψ = view(store.Ψ, 1:r, 1:M)
    legendre_polys!(Ψ, σ)
    A = view(store.A, 1:r)
    c = 2^(1-α)
    for j = 1:r
        s = zero(T)
        for m = 1:M
            s += wσ[m] * ( 3 - σ[m] )^(α-1) * Ψ[j,m]
        end
        A[j] = c * s
    end

    Mσ = ceil(Integer, r/2)
    σ, wσ = rule(store.jacobi1[Mσ])
    Ψ = view(store.Ψ, 1:r, 1:Mσ)
    legendre_polys!(Ψ, σ)
    B = view(store.B, 1:r)
    for j = 1:r
        s = zero(T)
        for m = 1:Mσ
            s+= wσ[m] * Ψ[j,m]
        end
        B[j] = 2^(1-α) * s
    end

    Mτ = max(1, ceil(Integer, r-1))
    Mσ = Mτ
    Mz = M
    z, wz = rule(store.unitlegendre[Mz])
    σ, wσ = rule(store.jacobi3[Mσ])
    τ, wτ = rule(store.jacobi4[Mτ])
    dΨ = view(store.dΨ, 1:r, 1:Mτ)
    deriv_legendre_polys!(dΨ, τ)
    C = view(store.C, 1:r, 1:r)
    for j = 1:r
        Ψ = view(store.Ψ, 1:j)
        for i = 1:r
            outer = zero(T)
            for mτ = 1:Mτ
                τ_ = τ[mτ]
                inner = zero(T)
                for mz = 1:Mz
                    z_ = z[mz]
                    legendre_polys!(Ψ, 1 - z_*(1+τ_))
                    inner += wz[mz] * (1+z_)^(α-1) * Ψ[j]
                end
                outer += wτ[mτ] * dΨ[i,mτ] * inner
            end
            C[i,j] = 2^(1-α) * outer
        end
    end
    Ψ = view(store.Ψ, 1:r, 1:Mσ)
    legendre_polys!(Ψ, σ)
    for j = 1:r
        for i = 1:r
            dΨ = view(store.dΨ, 1:i)
            outer = zero(T)
            for mσ = 1:Mσ
                σ_ = σ[mσ]
                inner = zero(T)
                for mz = 1:Mz
                    z_ = z[mz]
                    deriv_legendre_polys!(dΨ, z_*(1-σ_)-1)
                    inner += wz[mz] * (z_+1)^(α-1) * dΨ[i]
                end
                outer += wσ[mσ] * Ψ[j,mσ] * inner
            end
            C[i,j] += 2^(1-α) * outer
        end
    end

    c = 1 / ( 2Γ(α) )
    for j = 1:r
        pow = one(T)
        for i = 1:r
            pow = -pow
            H1[i,j] = c * ( A[j] + pow * B[j] - C[i,j] )
        end
    end
end

# H^ℓbar for ℓbar in rng, first version (rng must not include ℓbar ≤ 1)
function coef_H_uniform_ver1!(H::OffsetVector{Matrix{T}}, rng::UnitRange,
                              M::Integer, store::Store{T}
                             ) where T <: AbstractFloat
    @argcheck rng.start ≥ 2
    r = size(H[rng.start], 1)
    A = view(store.A, 1:r)
    B = view(store.B, 1:r)
    C = view(store.C, 1:r, 1:r)

    σ, wσ = rule(store.legendre[M])
    τ, wτ = σ, wσ
    Ψ  = view(store.Ψ, 1:r, 1:M)
    dΨ = view(store.dΨ, 1:r, 1:M)
    legendre_polys!(Ψ, σ)
    deriv_legendre_polys!(dΨ, τ)
    α = store.α
    c = 1 / ( 2Γ(α) )
    for ℓbar in rng

        for j = 1:r
            s = zero(T)
            for m = 1:M
                Δ = ( 1 - σ[m] ) / (2ℓbar)
                s += wσ[m] * ( 1 + Δ )^(α-1) * Ψ[j,m]
            end
            A[j] = s
        end

        for j = 1:r
            s = zero(T)
            for m = 1:M
                 Δ = ( 1 + σ[m] ) / ( 2ℓbar )
                 s += wσ[m] * (1-Δ)^(α-1) * Ψ[j,m]
            end
            B[j] = s
        end

        τ, wτ = σ, wσ
        for j = 1:r, i = 1:r
            outer = zero(T)
            for mτ = 1:M
                inner = zero(T)
                for mσ = 1:M
                    Δ = ( τ[mτ] - σ[mσ] ) / ( 2ℓbar )
                    inner += wσ[mσ] * (1+Δ)^(α-1) * Ψ[j,mσ]
                end
                outer += wτ[mτ] * dΨ[i,mτ] * inner
            end
            C[i,j] = outer
        end

        for j = 1:r
            pow = one(T)
            for i = 1:r
                pow = -pow
                H[ℓbar][i,j] = c * ℓbar^(α-1) * ( A[j] + pow * B[j] - C[i,j] )
            end
        end
    end
end

# H^ℓbar for ℓbar in rng, second version (rng must not include ℓbar ≤ 1)
function coef_H_uniform_ver2!(H::OffsetVector{Matrix{T}}, rng::UnitRange,
                              M::Integer, store::Store{T}
                             ) where T <: AbstractFloat
    @argcheck rng.start ≥ 2
    r = size(H[rng.start], 1)
    σ, wσ = rule(store.legendre[M])
    τ, wτ = σ, wσ
    Ψ  = view(store.Ψ, 1:r, 1:M)
    legendre_polys!(Ψ, σ)
    α = store.α
    c = - (1-α) / ( 4Γ(α) )
    for ℓbar in rng
        for j = 1:r, i = 1:r
            outer = zero(T)
            for mτ = 1:M
                inner = zero(T)
                for mσ = 1:M
                    Δ = ( τ[mτ] - σ[mσ] ) / ( 2ℓbar )
                    inner += wσ[mσ] * (1+Δ)^(α-2) * Ψ[j,mσ]
                end
                outer += wτ[mτ] * Ψ[i,mτ] * inner
            end
            H[ℓbar][i,j] = c * ℓbar^(α-2) * outer
        end
    end
end

"""
    coef_H(r)

Returns diagonal matrix `H` for the classical case `α=1`.
"""
function coef_H(::Type{T}, r::Integer) where T <: AbstractFloat
    d = T[ one(T) / (2j+1) for j = 0:r-1 ]
    return LinearAlgebra.Diagonal(d)
end

ext: jl


text:
"""
Case 9:
This case study a three bus system with 1 machine (One d- One q-: 4th order model), a VSM of 19 states and an infinite source.
The perturbation increase the reference power (analogy for mechanical power) of the inverter from 1.0 to 1.2.
"""

##################################################
############### LOAD DATA ########################
##################################################

include(joinpath(TEST_FILES_DIR, "data_tests/test09.jl"))

##################################################
############### SOLVE PROBLEM ####################
##################################################

@testset "Test 09 VSM Inverter and OneDoneQ ResidualModel" begin
    path = (joinpath(pwd(), "test-09"))
    !isdir(path) && mkdir(path)
    #time span
    tspan = (0.0, 20.0)
    case_inv = collect(PSY.get_components(PSY.DynamicInverter, threebus_sys))[1]
    case_gen = collect(PSY.get_components(PSY.DynamicGenerator, threebus_sys))[1]

    #Define Fault using Callbacks
    Pref_change = ControlReferenceChange(1.0, case_inv, :P_ref, 1.2)
    gen_trip = GeneratorTrip(1.5, case_gen)
    try
        #Define Simulation Problem
        sim = Simulation(
            ResidualModel,
            threebus_sys, # system
            path,
            tspan,
            [Pref_change, gen_trip],
        )

        # Test Initial Conditions
        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test09_x0_init
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end
        @test (diff[1] < 1e-3)

        # Obtain small signal results for initial conditions
        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        # Test Eigenvalues
        @test LinearAlgebra.norm(eigs - test09_eigvals) < 1e-3

        #Solve problem
        @test execute!(sim, IDA(), dtmax = 0.02) == PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        #Obtain data for angles
        t, θ_oc = get_state_series(results, ("generator-103-1", :θ_oc))
        @test length(t) == length(θ_oc)
        t, ω_oc = get_state_series(results, ("generator-103-1", :ω_oc))
        @test length(t) == length(ω_oc)
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end

@testset "Test 09 VSM Inverter and OneDoneQ MassMatrixModel" begin
    path = (joinpath(pwd(), "test-09"))
    !isdir(path) && mkdir(path)
    #time span
    tspan = (0.0, 20.0)
    case_inv = collect(PSY.get_components(PSY.DynamicInverter, threebus_sys))[1]
    case_gen = collect(PSY.get_components(PSY.DynamicGenerator, threebus_sys))[1]

    #Define Fault using Callbacks
    Pref_change = ControlReferenceChange(1.0, case_inv, :P_ref, 1.2)
    gen_trip = GeneratorTrip(1.5, case_gen)
    try
        #Define Simulation Problem
        sim = Simulation(
            MassMatrixModel,
            threebus_sys, # system
            path,
            tspan,
            [Pref_change, gen_trip],
        )

        # Test Initial Conditions
        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test09_x0_init
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end
        @test (diff[1] < 1e-3)

        # Obtain small signal results for initial conditions
        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        # Test Eigenvalues
        @test LinearAlgebra.norm(eigs - test09_eigvals) < 1e-3

        #Solve problem
        @test execute!(sim, Rodas4(), dtmax = 0.02) == PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        #Obtain data for angles
        t, θ_oc = get_state_series(results, ("generator-103-1", :θ_oc))
        @test length(t) == length(θ_oc)
        t, ω_oc = get_state_series(results, ("generator-103-1", :ω_oc))
        @test length(t) == length(ω_oc)

        power = PSID.get_activepower_series(results, "generator-103-1")
        rpower = PSID.get_reactivepower_series(results, "generator-103-1")
        @test isa(power, Tuple{Vector{Float64}, Vector{Float64}})
        @test isa(rpower, Tuple{Vector{Float64}, Vector{Float64}})

    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end


text:
@testset "layers" begin
    df = (x = rand(1000), y = rand(1000), c = rand(["a", "b", "c"], 1000))
    d = mapping(:x, :y, color=:c)
    s = visual(color=:red) + mapping(markersize=:c)
    layers = data(df) * d * s
    @test layers[1].transformation isa AlgebraOfGraphics.Visual
    @test layers[1].transformation.attributes[:color] == :red
    @test layers[1].positional == Any[:x, :y]
    @test layers[1].named == NamedArguments((color=:c,))
    @test layers[1].data == df
end

@testset "process_mappings" begin
    df = (x=rand(1000), y=rand(1000), z=rand(1000), c=rand(["a", "b", "c"], 1000))
    d = mapping(:x => exp, [:y, :z], color=:c, marker = dims(1) => renamer(["a", "b"]))
    layer = data(df) * d
    processedlayer = AlgebraOfGraphics.process_mappings(layer)
    @test processedlayer.positional[1] == fill(map(exp, df.x))
    @test processedlayer.positional[2] == [df.y, df.z]
    @test processedlayer.primary[:color] == fill(df.c)
    @test processedlayer.primary[:marker] == [fill(Sorted(1, "a")), fill(Sorted(2, "b"))]
    @test processedlayer.named == NamedArguments((;))
    @test processedlayer.labels[1] == fill("x")
    @test processedlayer.labels[2] ==  ["y", "z"]
    @test processedlayer.labels[:color] == fill("c")
    @test processedlayer.labels[:marker] == ""
end

@testset "grouping" begin
    df = (x=rand(1000), y=rand(1000), z=rand(1000), w=rand(1000), c=rand(["a", "b", "c"], 1000))
    df.c[1:3] .= ["a", "b", "c"] # ensure all three values exist
    d = mapping(:x => exp, [:y, :z], color=:c, marker=dims(1) => t -> ["1", "2"][t], markersize=:w)
    layer = data(df) * d
    processedlayer = AlgebraOfGraphics.ProcessedLayer(layer)
    processedlayers = map(CartesianIndices(AlgebraOfGraphics.shape(processedlayer))) do c
        primary, positional, named = map((processedlayer.primary, processedlayer.positional, processedlayer.named)) do tup
            return map(v -> v[c], tup)
        end
        labels = map(l -> AlgebraOfGraphics.getnewindex(l, c), processedlayer.labels)
        return ProcessedLayer(processedlayer; primary, positional, named, labels)
    end
    @test length(processedlayers) == 6
    for i in 1: 6
        @test processedlayers[i].plottype === Any
        @test isempty(processedlayers[i].attributes)
        @test processedlayers[i].labels[1] == "x"
        @test processedlayers[i].labels[2] == (i ≤ 3 ? "y" : "z")
        @test processedlayers[i].labels[:color] == "c"
        @test processedlayers[i].labels[:marker] == ""
        @test processedlayers[i].labels[:markersize] == "w"
    end

    @test processedlayers[1].primary[:color] == "a"
    @test processedlayers[1].primary[:marker] == "1"
    @test processedlayers[1].positional[1] == exp.(df.x[df.c .== "a"])
    @test processedlayers[1].positional[2] == df.y[df.c .== "a"]
    @test processedlayers[1].named[:markersize] == df.w[df.c .== "a"]

    @test processedlayers[2].primary[:color] == "b"
    @test processedlayers[2].primary[:marker] == "1"
    @test processedlayers[2].positional[1] == exp.(df.x[df.c .== "b"])
    @test processedlayers[2].positional[2] == df.y[df.c .== "b"]
    @test processedlayers[2].named[:markersize] == df.w[df.c .== "b"]

    @test processedlayers[3].primary[:color] == "c"
    @test processedlayers[3].primary[:marker] == "1"
    @test processedlayers[3].positional[1] == exp.(df.x[df.c .== "c"])
    @test processedlayers[3].positional[2] == df.y[df.c .== "c"]
    @test processedlayers[3].named[:markersize] == df.w[df.c .== "c"]

    @test processedlayers[4].primary[:color] == "a"
    @test processedlayers[4].primary[:marker] == "2"
    @test processedlayers[4].positional[1] == exp.(df.x[df.c .== "a"])
    @test processedlayers[4].positional[2] == df.z[df.c .== "a"]
    @test processedlayers[4].named[:markersize] == df.w[df.c .== "a"]

    @test processedlayers[5].primary[:color] == "b"
    @test processedlayers[5].primary[:marker] == "2"
    @test processedlayers[5].positional[1] == exp.(df.x[df.c .== "b"])
    @test processedlayers[5].positional[2] == df.z[df.c .== "b"]
    @test processedlayers[5].named[:markersize] == df.w[df.c .== "b"]

    @test processedlayers[6].primary[:color] == "c"
    @test processedlayers[6].primary[:marker] == "2"
    @test processedlayers[6].positional[1] == exp.(df.x[df.c .== "c"])
    @test processedlayers[6].positional[2] == df.z[df.c .== "c"]
    @test processedlayers[6].named[:markersize] == df.w[df.c .== "c"]
end


text:
using ProgressMeter, MultivariateStats

function getdatapatch(data, pos, imsize; patchsize = 10, colorimage = false)
    if colorimage
        data_all_colors = reshape(reshape(data[1:imsize^2,:], imsize, imsize, size(data)[2])[pos[1]:pos[1] + patchsize - 1,
            pos[2]:pos[2] + patchsize - 1,:], patchsize^2, size(data)[2])
        for i in 2:3
            data_all_colors = vcat(data_all_colors, reshape(reshape(data[(i-1) * imsize^2 + 1:i * imsize^2,:], imsize, imsize, size(data)[2])[pos[1]:pos[1] + patchsize - 1,
                    pos[2]:pos[2] + patchsize - 1,:], patchsize^2, size(data)[2]))
        end
        data_all_colors
    else
        reshape(reshape(data, imsize, imsize, size(data)[2])[pos[1]:pos[1] + patchsize - 1,
            pos[2]:pos[2] + patchsize - 1,:], patchsize^2, size(data)[2])
    end
end

function getPCAscores(data, datatest, n_hidden_per_patch)
    (n_hidden_per_patch >= size(data)[1]) && error("cannot extract more PCs than (patch) input dimension")
    PCA_fit = fit(PCA, data, pratio = 1., maxoutdim = n_hidden_per_patch)
    return transform(PCA_fit, data), transform(PCA_fit, datatest)
end

function getICAscores(data, datatest, n_hidden_per_patch)
    (n_hidden_per_patch >= size(data)[1]) && error("cannot extract more ICs than (patch) input dimension")
    ICA_fit = fit(ICA, hcat(data,datatest), n_hidden_per_patch, tol = 2.5 * sqrt(n_hidden_per_patch), do_whiten=true)#false) #
    return transform(ICA_fit, data), transform(ICA_fit, datatest)
end

# output: smallimgs-array for shallow net!
function getpatchyhiddenreps(data, datatest, n_patches, n_hidden_per_patch;
    featurefunction = getPCAscores, patchsize = 10, randpos = true, colorimage = false)

    smallimgs = zeros(n_patches * n_hidden_per_patch, size(data)[2])
    smallimgstest = zeros(n_patches * n_hidden_per_patch, size(datatest)[2])

    if colorimage
        input_size = Int(size(data)[1]/3)
    else
        input_size = Int(size(data)[1])
    end
    imsize = Int(sqrt(input_size))
    print(string("calculate unsupervised hidden representations with function: ",string(featurefunction),"\n"))
    @showprogress for i in 1:n_patches
        if !randpos
            error("randpos must be true (in this state of the code)")
        elseif randpos
            if patchsize == imsize
                pos = [1,1]
            else
                pos = rand(1:imsize - patchsize,2)
            end
            patchdata = getdatapatch(data, pos, imsize; patchsize = patchsize, colorimage = colorimage)
            patchdatatest = getdatapatch(datatest, pos, imsize; patchsize = patchsize, colorimage = colorimage)
            smallimgs[(i-1) * n_hidden_per_patch + 1:i * n_hidden_per_patch,:],
             smallimgstest[(i-1) * n_hidden_per_patch + 1:i * n_hidden_per_patch,:] =
                        featurefunction(patchdata, patchdatatest, n_hidden_per_patch)
        end
    end
    return smallimgs ./ maximum(smallimgs), smallimgstest ./ maximum(smallimgs)
end


text:
import Printf: @sprintf

const srcpath = joinpath(splitdir(Base.source_path())[1], "src", "blend2d", "src")
const outpath = joinpath(splitdir(Base.source_path())[1], "..", "src", "gen")
const header = "# This file has been auto-generated, do not edit!\n"

const enumsignore = ["BLByteOrder"]
const funcsignore = [
	v -> !(v in [
		# implemented manually
		"blImageCodecBuiltInCodecs",
		# BLGlyphRun
		"blContextFillGlyphRunI",
		"blContextFillGlyphRunD",
		"blContextStrokeGlyphRunI",
		"blContextStrokeGlyphRunD",
		# va_list
		"blRuntimeMessageVFmt",
		"blStringApplyOpFormatV",
		# BLFormatInfo
		"blFormatInfoSanitize",
		# BLStrokeOptionsCore
		"blContextGetStrokeOptions",
		"blContextSetStrokeOptions",
		"blPathAddStrokedPath",
		# BLImageScaleOptions
		"blImageScale",
		# BLContextHints
		"blContextSetHints",
		# unused
		"blRuntimeAssertionFailure"
	]),
	v -> !startswith(v, "blFont"),
	v -> !startswith(v, "blGlyph"),
	v -> !startswith(v, "blGradient"),
	v -> !startswith(v, "blPixelConverter"),
	v -> !startswith(v, "blStroke"),
	v -> !startswith(v, "blFile"),
	v -> !endswith(v, "FromData"),
]
const structsignore = [
	# implemented manually
	"BLMatrix2D"
]

function sanitize_name(s)
	if s == "end"
		return "end_"
	end
	return s
end

function parse_header(filename)
	content = join(readlines(filename), "\n")

	# enums
	enums = collect(eachmatch(r"^BL_DEFINE_ENUM\((?<name>\w+)\)\s*\{(?<values>.*?)\}"sm, content))
	enums = filter(v -> !(v[:name] in enumsignore), enums)
	enums = map(function(enum)
		vals = collect(eachmatch(r"^\s*(?<item>BL_\w+)(\s*?=\s*?(?<const>[x0-9a-fA-F]+)u?)?"sm, enum[:values]))
		vals = map(v -> Dict(:item => v[:item], :const => v[:const] == nothing ? "" : String(v[:const])), vals)
		for (index, item) in enumerate(vals)
			if item[:const] == ""
				item[:const] = index == 1 ? "0" : @sprintf("0x%08x", parse(Int, vals[index - 1][:const]) + 1)
			else
				item[:const] = @sprintf("0x%08x", parse(Int, item[:const]))
			end
		end
		Dict(:name => enum[:name], :values => vals)
	end, enums)

	# funcs
	# keep closing ) in 'params' so we can match the end of the params or comma later on
	funcs = collect(eachmatch(r"^BL_API_C\s+(?<return>.+?)\s+BL_CDECL\s+(?<name>\w+)\s*\((?<params>.+?\)).*?;$"sm, content))
	funcs = filter(v -> all(t -> t(v[:name]), funcsignore), funcs)
	funcs = map(function(func)
		params = collect(eachmatch(r"\s*(?<type>.+?)\s+(?<name>\w+?)\s*[,\)]", func[:params]))
		params = map(p -> Dict(:name => p[:name], :type => p[:type]), params)
		Dict(:name => func[:name], :return => func[:return], :params => params)
	end, funcs)

	# structs
	structs = collect(eachmatch(r"^struct\s+(?<name>\w+)\s*\{(?<fields>.*?)\}\s*;"sm, content))
	structs = map(s -> Dict(:name => s[:name], :fields => s[:fields]), structs)
	structs = map(function(s)
		fields = s[:fields]
		fields = replace(fields, r"#if.+?#endif"sm => "")
		fields = replace(fields, r"^\s*//.*?$"sm => "")
		fields = collect(eachmatch(r"\s*(?<type>.+?)\s+(?<name>\w+?)\s*(?<array>\[\d*\])?\s*;"sm, fields))
		fields = map(p -> Dict(:name => p[:name], :type => p[:type]), fields)
		Dict(:name => s[:name], :fields => fields)
	end, structs)

	Dict(:enums => enums, :funcs => funcs, :structs => structs)
end

function ctype_match(ctype)
	return match(r"(const\s+)?(?<basetype>[^*]+)(?<ptr>\*)?", ctype)
end

function is_mutable(structname)
	if any(v -> startswith(structname, v), ["BLBox", "BLRect", "BLPoint", "BLSize"])
		return false
	end
	return true
end

function default_value(type)
	if type == "Ptr{Cvoid}"
		return "Ptr{Cvoid}(0)"
	elseif startswith(type, "C")
		return "0"
	else
		return "$type()"
	end
end

function convert_type(structnames, enumnames, ctype, isccall)
	t = ctype_match(ctype)
	if t[:basetype] in structnames
		if endswith(t[:basetype], "Impl")
			return "Ptr{Cvoid}"
		elseif is_mutable(t[:basetype])
			return isccall ? "Ref{$(t[:basetype])}" : t[:basetype]
		else
			return t[:basetype]
		end
	elseif t[:basetype] in enumnames
		return t[:basetype]
	elseif t[:basetype] == "char"
		if t[:ptr] != nothing
			return isccall ? "Cstring" : "String"
		else
			return "Cchar"
		end
	elseif t[:basetype] == "uint8_t"
		return "Cuchar"
	elseif t[:basetype] == "uint16_t"
		return isccall ? "Cushort" : "Integer"
	elseif t[:basetype] == "uint32_t"
		return isccall ? "Cuint" : "Integer"
	elseif t[:basetype] == "int64_t"
		return isccall ? "Clonglong" : "Integer"
	elseif t[:basetype] == "uint64_t"
		return isccall ? "Culonglong" : "Integer"
	elseif t[:basetype] == "float"
		return "Cfloat"
	elseif t[:basetype] == "double"
		return "Cdouble"
	elseif t[:basetype] == "bool"
		return "Cuchar"
	elseif t[:basetype] == "char" && t[:ptr] != nothing
		return "Cstring"
	elseif t[:basetype] == "int"
		return isccall ? "Cint" : "Integer"
	elseif t[:basetype] == "size_t"
		return isccall ? "Csize_t" : "Integer"
	elseif t[:basetype] == "intptr_t"
		return "Culonglong"
	elseif t[:basetype] == "BLResult"
		return "BLResultCode"
	elseif t[:basetype] == "BLBitWord"
		return "Ptr{Cvoid}"
	elseif endswith(t[:basetype], "Impl")
		return "Ptr{Cvoid}"
	elseif t[:basetype] == "void"
		if t[:ptr] != nothing
			return "Ptr{Cvoid}"
		else
			return "Cvoid"
		end
	end
	error("unknown type [$(ctype)]")
end

function entry_point_parse_headers()
	mainheader = joinpath(srcpath, "blend2d.h")
	content = join(readlines(mainheader), "\n")
	matches = collect(eachmatch(r"#include \"(?<include>.+?)\""sm, content))
	includes = map(m -> m[:include], matches)
	results = Dict()
	for inc in includes
		incpath = joinpath(srcpath, inc)
		result = parse_header(incpath)
		merge!(append!, results, result)
	end

	exportnames = Set()
	enumnames = Set(map(v -> v[:name], results[:enums]))
	structnames = Set(map(v -> v[:name], results[:structs]))

	mkpath(outpath)
	open(joinpath(outpath, "enums.jl"), "w") do f
		write(f, header)
		for enum in results[:enums]
			push!(exportnames, enum[:name])
			write(f, "# $(enum[:name])\n")
			write(f, "const $(enum[:name]) = UInt32\n")
			for value in enum[:values]
				push!(exportnames, value[:item])
				write(f, "const $(value[:item]) = $(value[:const])\n")
			end
		end
	end

	mkpath(outpath)
	open(joinpath(outpath, "funcs.jl"), "w") do f
		write(f, header)
		for func in results[:funcs]
			push!(exportnames, func[:name])
			ret = convert_type(structnames, enumnames, func[:return], true)
			fparams = join(map(p -> "$(p[:name])::$(convert_type(structnames, enumnames, p[:type], false))", func[:params]), ", ")
			cparams = join(map(p -> "$(convert_type(structnames, enumnames, p[:type], true))", func[:params]), ", ")
			cparamnames = join(map(p -> p[:name], func[:params]), ", ")
			write(f, "function $(func[:name])($(fparams))\n")
			if size(func[:params], 1) > 0
				write(f, "\tccall((:$(func[:name]), libblend2d), $(ret), ($(cparams),), $(cparamnames))\n")
			else
				write(f, "\tccall((:$(func[:name]), libblend2d), $(ret), ())\n")
			end
			write(f, "end\n")
		end
	end

	structstokeep = Set()
	for func in results[:funcs]
		push!(structstokeep, ctype_match(func[:return])[:basetype])
		for p in func[:params]
			push!(structstokeep, ctype_match(p[:type])[:basetype])
		end
	end
	for stru in results[:structs]
		if stru[:name] in structstokeep && !endswith(stru[:name], "Core")
			for f in stru[:fields]
				push!(structstokeep, ctype_match(f[:type])[:basetype])
			end
		end
	end

	mkpath(outpath)
	open(joinpath(outpath, "structs.jl"), "w") do f
		write(f, header)
		for stru in results[:structs]
			if stru[:name] in structstokeep
				if stru[:name] in structsignore
					continue
				end
				push!(exportnames, stru[:name])
				try
					fields = join(map(f -> "\t$(sanitize_name(f[:name]))::$(convert_type(structnames, [], f[:type], true))", stru[:fields]), "\n")
					defaultvalues = join(map(f -> "$(default_value(convert_type(structnames, [], f[:type], true)))", stru[:fields]), ", ")
					cons = "$(stru[:name])() = new($defaultvalues)"
					if is_mutable(stru[:name])
						write(f, "mutable struct $(stru[:name])\n$fields\n\t$cons\nend\n\n")
					else
						defaultconsnames = join(map(f -> f[:name], stru[:fields]), ", ")
						cons = "$cons\n\t$(stru[:name])($defaultconsnames) = new($defaultconsnames)"
						write(f, "struct $(stru[:name])\n$fields\n\t$cons\nend\n\n")
					end
				catch e
					rethrow([ErrorException(stru[:name]), e])
				end
			end
		end
	end

	mkpath(outpath)
	open(joinpath(outpath, "export.jl"), "w") do f
		write(f, header)
		for exportname in exportnames
			write(f, "export $exportname\n")
		end
	end
end

entry_point_parse_headers()


text:
#Adding all the packages we need for this bot
import Pkg
Pkg.add(["OAuth", "HTTP", "JSON"])
using OAuth, HTTP, JSON

#Creating a function to filter out the tweets we want
function filter_tweets(tweets::Array)

    #List of tweets to be returned
    filtered_tweets = String[]

    #List of keywords we want the tweet to have
    keywords = ["brexit", "european union", "referendum", "conservative", "labour", "boris johnson", "united kingdom", "scotland", "ukip", "london", "manchester", "birmingham", "edinburgh", "british", "britain"]

    #Loop checking whether each tweet has contains any of the keywords, and if so, adds the tweet's id to the list of tweets to be returned
    for i in keywords
        for k = 1:length(tweets)
            if occursin(i, lowercase(tweets[k]["text"]))
                push!(filtered_tweets, tweets[k]["id_str"])
            end
        end
    end

    #Returns the list of the id of the tweets containing the keywords
    filtered_tweets
end

#Main function containing the body of the bot
function main()

    #Returns the JSON from the OAuth Request to get the most recent <number of tweets specified in the variable below> tweets
    tweets = JSON.parse(String(oauth_get("https://api.twitter.com/1.1/statuses/home_timeline.json", Dict("count" => "$no_of_tweets")).body))

    #Passes the tweets to be filtered and returns the tweets we want
    filtered_tweets = filter_tweets(tweets)

    #If there are no new tweets, we will give the user a message to let them know
    if length(filtered_tweets) == 0
        println("sorry, no new british news :(")

    #If there are new tweets, we want to retweet each and one of them
    else
        for id in filtered_tweets
            try
                oauth_post("https://api.twitter.com/1.1/statuses/retweet/$id.json", Dict("id" => "$id")) #runs the retweet request

            #Catches the error given when it tries to retweet the same tweet for the 2nd time
            catch
                continue #Skips the rest of the current loop iteration to move onto the next retweet
            end
        end
        println("added ", length(filtered_tweets), " tweets") #tells the user how many tweets were added
    end
end

#Variables containing the keys and tokens required to power the account the bot is running on
consumer_key = ARGS[1]
secret_consumer_key = ARGS[2]
access_token = ARGS[3]
secret_access_token = ARGS[4]

#Rewriting GET and POST functions to make life easier: we don't need to call the keys and tokens every time!
oauth_get(endpoint::String, options::Dict) = oauth_request_resource(endpoint, "GET", options, consumer_key, secret_consumer_key, access_token, secret_access_token)
oauth_post(endpoint::String, options::Dict) = oauth_request_resource(endpoint, "POST", options, consumer_key, secret_consumer_key, access_token, secret_access_token)

#The number of the most recent tweets the bot looks through on its home timeline every 8 minutes
no_of_tweets = 15

#Runs the main function
main()

max_stars_repo_name: kim-fung/UKNewsTwitterBot

max_issues_repo_name: kim-fung/UKNewsTwitterBot

max_forks_repo_name: kim-fung/UKNewsTwitterBot


text:
using LinearAlgebra

"""
    iscommuting(A::AbstractMatrix, B::AbstractMatrix)

returns true if [`AbstractMatrix`](@ref) `A` and `B` commute.
"""
function iscommuting(A::AbstractMatrix, B::AbstractMatrix)
    # require compatible dimensions
    if size(A, 2) == size(B, 1) && size(A, 1) == size(B, 2)
        return (A * B ≈ B * A)
    else
        return false
    end
end

# gates of same type do not necessarily commute, see, e.g., MatrixGate

iscommuting(::XGate, ::XGate) = true
iscommuting(::YGate, ::YGate) = true
iscommuting(::ZGate, ::ZGate) = true

iscommuting(::XGate, ::YGate) = false
iscommuting(::YGate, ::XGate) = false
iscommuting(::XGate, ::ZGate) = false
iscommuting(::ZGate, ::XGate) = false
iscommuting(::YGate, ::ZGate) = false
iscommuting(::ZGate, ::YGate) = false

iscommuting(::RxGate, ::RxGate) = true
iscommuting(::RyGate, ::RyGate) = true
iscommuting(::RzGate, ::RzGate) = true

iscommuting(::XGate,  ::RxGate) = true
iscommuting(::RxGate, ::XGate)  = true

iscommuting(::RyGate, ::YGate) = true
iscommuting(::YGate, ::RyGate) = true

iscommuting(::RzGate, ::ZGate) = true
iscommuting(::ZGate, ::RzGate) = true

# entanglement gates commute with each other
iscommuting(::EntanglementXXGate, ::EntanglementXXGate) = true
iscommuting(::EntanglementXXGate, ::EntanglementYYGate) = true
iscommuting(::EntanglementXXGate, ::EntanglementZZGate) = true
iscommuting(::EntanglementYYGate, ::EntanglementXXGate) = true
iscommuting(::EntanglementYYGate, ::EntanglementYYGate) = true
iscommuting(::EntanglementYYGate, ::EntanglementZZGate) = true
iscommuting(::EntanglementZZGate, ::EntanglementXXGate) = true
iscommuting(::EntanglementZZGate, ::EntanglementYYGate) = true
iscommuting(::EntanglementZZGate, ::EntanglementZZGate) = true

# cover MatrixGate, too
iscommuting(A::AbstractGate, B::AbstractGate) = iscommuting(Qaintessent.sparse_matrix(A), Qaintessent.sparse_matrix(B))

"""
    iscommuting(A::CircuitGate{L,ControlledGate{G}}, B::CircuitGate{M,ControlledGate{G}}) where {L,M,G<:AbstractGate}

returns true if controlled [`CircuitGate`](@ref) `A` and `B` commute.
"""
function iscommuting(A::CircuitGate{L,ControlledGate{G}}, B::CircuitGate{M,ControlledGate{G}}) where {L,M,G<:AbstractGate}
    N = maximum((maximum(B.iwire), maximum(A.iwire)))
    # check whether only control wires overlap, if at all
    if (length(intersect(A.iwire[end - L + 1:end], B.iwire)) == 0 &&
        length(intersect(B.iwire[end - M + 1:end], A.iwire)) == 0)
        return true
    end

    if A.iwire == B.iwire
        return iscommuting(A.gate, B.gate)
    end

    # TODO: more fine-grained case switches for partially overlapping 'iwire'
    iscommuting(Qaintessent.sparse_matrix(A, N), Qaintessent.sparse_matrix(B, N))
end

"""
    iscommuting(A::CircuitGate{L,G}, B::CircuitGate{M,H}) where {L,M,G,H}

returns true if general [`CircuitGate`](@ref) `A` and `B` commute.
"""
function iscommuting(A::CircuitGate{L,G}, B::CircuitGate{M,H}) where {L,M,G,H}
    N = maximum((maximum(B.iwire), maximum(A.iwire)))
    if length(intersect(A.iwire, B.iwire)) == 0
        return true
    end

    if A.iwire == B.iwire
        return iscommuting(A.gate, B.gate)
    end

    # TODO: more fine-grained case switches for partially overlapping 'iwire'
    iscommuting(Qaintessent.sparse_matrix(A,N), Qaintessent.sparse_matrix(B,N))
end

"""
    iscommuting(A::MeasurementOperator{L,G}, B::MeasurementOperator{M,H}) where {L,M,G,H}

returns true if [`MeasurementOperator`](@ref) `A` and `B` commute.
"""
function iscommuting(A::MeasurementOperator{L,G}, B::MeasurementOperator{M,H}) where {L,M,G,H}
    N = maximum((maximum(A.iwire), maximum(B.iwire)))
    if length(intersect(A.iwire, B.iwire)) == 0
        return true
    end

    if A.iwire == B.iwire
        return iscommuting(A.operator, B.operator)
    end

    iscommuting(sparse_matrix(A, N), sparse_matrix(B, N))
end

function iscommuting(A::MeasurementOperator{M,G}, B::MeasurementOperator{L,H}) where {L,M,G<:AbstractGate,H<:AbstractMatrix}
    if length(intersect(A.iwire, B.iwire)) == 0
        return true
    end

    N = maximum((maximum(A.iwire), maximum(B.iwire)))
    iscommuting(sparse_matrix(A, N), sparse_matrix(B, N))
end

iscommuting(A::MeasurementOperator{M,G}, B::MeasurementOperator{L,H}) where {L,M,G<:AbstractMatrix,H<:AbstractGate} = iscommuting(B, A)


text:
########## Parameters for scripts ##########

###### bspline-MLE ######

#### cross-time only ####
maxscl = 60

parms = (
    estimator = "bspline_MLE_time",
    method = :optim,  # method of optimizer

    # Parameters for B-Spline wavelet transform
    sclrng = [maxscl],  # index of working scale
    vm = 1,  # vanishing moment of the wavelet
    fmode = :causal,  # mode of filtration
    ρ = 0,  # threshold for truncation of small eigen values in the covariance matrix
    partial = false,  # true for whittle estimator

    dfct = maxscl,  # down-sampling factor
    ssize = 5,  # size of sub window, set ssize=wsize/dfct to disable sub window
    dlen = 1,  # length of decorrelation of sub window

    wsize = maxscl*15,  # size of moving window
    pov = maxscl,  # step of moving window
    boundary = :hard,  # boundary condition
)

base_estimator = x -> FracFin.bspline_MLE_estim(x, 2*parms.sclrng, parms.vm, parms.ssize, parms.dlen, parms.dlen; ρ=parms.ρ, partial=parms.partial)

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
function create_drawing(img::Array{UInt32,2})
    return Cairo.CairoRGBSurface(img)
end

function create_contex(surface::Cairo.CairoSurface)
    return Cairo.CairoContext(surface)
end

function draw_vectors(ctx::Cairo.CairoContext, vectors, color, factor = 1)
    Cairo.save(ctx)
    Cairo.set_source_rgb(ctx, color...)
    for line in vectors
        diff = line[2] - line[1]
        Cairo.move_to(ctx, line[1]...)
        Cairo.line_to(ctx, line[1] + diff * factor...)
    end
    Cairo.set_line_width(ctx, 3.0)
    Cairo.stroke(ctx)
    return ctx
end

function draw_line(ctx::Cairo.CairoContext, line, color = (0, 0, 0), thickness = 2.0)
    Cairo.save(ctx)
    Cairo.set_source_rgb(ctx, color...)
    Cairo.move_to(ctx, line[1]...)
    Cairo.line_to(ctx, line[2]...)
    Cairo.set_line_width(ctx, thickness)
    Cairo.stroke(ctx)
    return ctx
end

function draw_points(ctx::Cairo.CairoContext, points, color = (0, 0, 0), factor = 1)
    d = 5
    Cairo.save(ctx)
    Cairo.set_source_rgb(ctx, color...)
    for point in points
        Cairo.move_to(ctx, point + [-d, -d]...)
        Cairo.line_to(ctx, point + [d, d]...)
        Cairo.move_to(ctx, point + [-d, d]...)
        Cairo.line_to(ctx, point + [d, -d]...)
    end
    Cairo.set_line_width(ctx, 3.0)
    Cairo.stroke(ctx)
    return ctx
end

function draw_point(ctx::Cairo.CairoContext, point, color, factor = 1)
    d = 5
    Cairo.save(ctx)
    Cairo.set_source_rgb(ctx, color...)
    Cairo.move_to(ctx, point + [-d, -d]...)
    Cairo.line_to(ctx, point + [d, d]...)
    Cairo.move_to(ctx, point + [-d, d]...)
    Cairo.line_to(ctx, point + [d, -d]...)
    Cairo.set_line_width(ctx, 2.0)
    Cairo.stroke(ctx)
    return ctx
end

function draw_text(ctx::Cairo.CairoContext, text, point, offset, fontsize, color)
    fill_color = color
    line_color = [1, 1, 1] - color
    Cairo.save(ctx)
    Cairo.select_font_face(ctx, "Sans", Cairo.FONT_SLANT_NORMAL, Cairo.FONT_WEIGHT_BOLD)
    Cairo.set_font_size(ctx, fontsize)
    Cairo.move_to(ctx, point + offset...)
    Cairo.rotate(ctx, -pi / 2)
    Cairo.scale(ctx, -1, 1)
    Cairo.text_path(ctx, text)
    Cairo.set_source_rgb(ctx, fill_color...)
    Cairo.fill_preserve(ctx)
    Cairo.set_source_rgb(ctx, line_color...)
    Cairo.set_line_width(ctx, 1.0)
    Cairo.scale(ctx, -1, 1)
    Cairo.rotate(ctx, pi / 2)
    Cairo.stroke(ctx)
    return ctx
end

function draw_indices(ctx::Cairo.CairoContext, points, offset, fontsize, color)
    fill_color = color
    line_color = [1, 1, 1] - color
    Cairo.save(ctx)
    Cairo.select_font_face(ctx, "Sans", Cairo.FONT_SLANT_NORMAL, Cairo.FONT_WEIGHT_BOLD)
    Cairo.set_font_size(ctx, fontsize)
    for (k, point) in enumerate(points)
        Cairo.move_to(ctx, point + offset...)
        Cairo.rotate(ctx, -pi / 2)
        Cairo.scale(ctx, -1, 1)
        Cairo.text_path(ctx, "$k")
        Cairo.set_source_rgb(ctx, fill_color...)
        Cairo.fill_preserve(ctx)
        Cairo.set_source_rgb(ctx, line_color...)
        Cairo.set_line_width(ctx, 1.0)
        Cairo.scale(ctx, -1, 1)
        Cairo.rotate(ctx, pi / 2)
    end
    Cairo.stroke(ctx)
    return ctx
end

function get_drawing(surface::Cairo.CairoSurface)
    return convert(Array{RGB24}, surface.data)
end

function draw_reference(ctx, h, w, factor)
    reference = 1
    margin = 10
    a = [h - margin, w - margin - reference * factor]
    b = [h - margin, w - margin]
    color = [0.9, 0.9, 0.9]
    draw_line(ctx, (a, b), color, 3.0)
    draw_text(ctx, "1", a + [-20, reference * factor / 2], [0, 0], 24.0, color)
    return ctx
end


text:
using Distributions

function betarnd(a,b,n,m)
	d = Beta(a,b)
	r = rand(d,int(n*m))
	reshape(r, n,m)
end

function randVonMisesFisher(m, n, kappa, mu = None)
	#% randVonMISESFISHERM Random number generation from von Mises Fisher
	#% distribution.
	#% X = randVonMisesFisherm(m, n, kappa) returns n samples of random unit
	#% directions in m dimensional space, with concentration parameter kappa,
	#% and the direction parameter mu = e_m
	#% X = randVonMisesFisherm(m, n, kappa, mu) with direction parameter mu
	#% (m-dimensional column unit vector)
	#%
	#% Sungkyu Jung, Feb 3, 2010.
	#
	# converted into Julia, Daniel Perry, April 2015,
	#
	# see http://www.stat.pitt.edu/sungkyu/MiscPage.html
	#     http://www.stat.pitt.edu/sungkyu/software/randVonMisesFisherm.zip
	#     http://www.stat.pitt.edu/sungkyu/software/randVonMisesFisher3.pdf

	if mu == None
		muflag = false;
	else
		muflag = true;
	end

	if m < 2;
			println("Message from randVonMisesFisherm.m: dimension m must be > 2")
			println("Message from randVonMisesFisherm.m: Set m to be 2")
			m = 2;
	end

	if kappa < 0;
			println("Message from randVonMisesFisherm.m: kappa must be >= 0");
			println("Message from randVonMisesFisherm.m: Set kappa to be 0");
			kappa = 0;
	end

	#%
	#% the following algorithm is following the modified Ulrich's algorithm
	#% discussed by Andrew T.A. Wood in "SIMULATION OF THE VON MISES FISHER
	#% DISTRIBUTION", COMMUN. STATIST 23(1), 1994.

	#% step 0 : initialize
	b = (-2*kappa + sqrt(4*kappa^2 + (m-1)^2))/(m-1);
	x0 = (1-b)/(1+b);
	c = kappa*x0 + (m-1)*log(1-x0^2);

	#% step 1 & step 2
	nnow = n;
	w = [];
	#%cnt = 0;
	while true
			ntrial = int(max(round(nnow*1.2),nnow+10))
			Z = betarnd((m-1)/2,(m-1)/2,ntrial,1);
			U = rand(ntrial,1);
			W = (1-(1+b)*Z)./(1-(1-b)*Z);

			indicator = (kappa*W + (m-1)*log(1-x0*W) - c) .>= log(U)
			if sum(indicator) >= nnow
					w1 = W[find(indicator)]
					w = [w ;w1[1:nnow]];
					break;
			else
					w = [w ; W[find(indicator)]];
					nnow = nnow-sum(indicator);
					#%cnt = cnt+1;disp(['retrial' num2str(cnt) '.' num2str(sum(indicator))]);
			end
	end

	#% step 3
	V = UNIFORMdirections(m-1,n);
	X = [repmat(sqrt(1-w'.^2),m-1,1).*V ;w'];

	if muflag
			mu = mu / norm(mu);
			X = rotMat(mu)'*X;
	end
	return X
end

function UNIFORMdirections(m,n)
	#% generate n uniformly distributed m dim'l random directions
	#% Using the logic: "directions of Normal distribution are uniform on sphere"

	V = zeros(m,n);
	nr = randn(m,n); #%Normal random
	for i=1:n
			while true
					ni=(nr[:,i]'*nr[:,i])[1]; #% length of ith vector
					#% exclude too small values to avoid numerical discretization
					if ni < 1e-10
							#% so repeat random generation
							nr[:,i]=randn(m,1);
					else
							 V[:,i]=nr[:,i]/sqrt(ni);
							break;
					end
			end
	end

	return V
end

function rotMat(b, a = None, alpha = None)
	#% ROTMAT returns a rotation matrix that rotates unit vector b to a
	#%
	#%   rot = rotMat(b) returns a d x d rotation matrix that rotate
	#%   unit vector b to the north pole (0,0,...,0,1)
	#%
	#%   rot = rotMat(b,a ) returns a d x d rotation matrix that rotate
	#%   unit vector b to a
	#%
	#%   rot = rotMat(b,a,alpha) returns a d x d rotation matrix that rotate
	#%   unit vector b towards a by alpha (in radian)
	#%
	#%    See also .
	#
	#% Last updated Nov 7, 2009
	#% Sungkyu Jung

	s1 = size(b,1)
	s2 = size(b,2)
	d = max(s1,s2);
	b= b/norm(b);

	if s1<=s2
		b = b'
	end

	if a == None
			a = [zeros(d-1,1); 1];
	end
	if alpha == None
			alpha = acos((a'*b)[1]);
	end

	if abs((a'*b)[1] - 1) < 1e-15
		rot = eye(d)
		return rot
	end
	if abs((a'*b)[1] + 1) < 1e-15
		rot = -eye(d)
		return rot
	end

	c = b - a * (a'*b); c = c / norm(c);
	A = a*c' - c*a' ;

	rot = eye(d) + sin(alpha)*A + (cos(alpha) - 1)*(a*a' +c*c');
end


text:
using Printf

function abc(str::AbstractString, list)
    isempty(str) && return true
    for i in eachindex(list)
        str[end] in list[i] &&
            any([abc(str[1:end-1], deleteat!(copy(list), i))]) &&
            return true
    end
    return false
end

let test = ["A", "BARK","BOOK","TREAT","COMMON","SQUAD","CONFUSE"],
    list = ["BO","XK","DQ","CP","NA","GT","RE","TG","QD","FS",
            "JW","HU","VI","AN","OB","ER","FS","LY","PC","ZM"]
    for str in test
        @printf("%-8s |  %s\n", str, abc(str, list))
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@coeff function K_H(reactor::Reactor)
  cur_K = 2.272

  cur_K *= K_L * K_beta ^ 1.34 * K_F ^ 1.17

  cur_K /= A ^ 0.19 * epsilon ^ 0.58 * kappa ^ 0.78 * K_G ^ 0.41 * P_F ^ 0.48

  cur_K *= ( Q / ( Q + 5 ) ) ^ 0.31

  cur_K
end


text:
"""
Base module for Yao.
"""
module YaoBase

using LinearAlgebra, LuxurySparse, SparseArrays, Random
using Reexport

@reexport using YaoAPI
using YaoAPI

import YaoAPI:
    isunitary,
    isreflexive,
    iscommute,
    AbstractRegister,
    AdjointRegister,
    AbstractBlock,
    PostProcess,
    NotImplementedError,
    LocationConflictError,
    QubitMismatchError,
    instruct!,
    focus!,
    relax!,
    nqubits,
    nremain,
    nactive,
    nbatch,
    viewbatch,
    addbits!,
    insert_qubits!,
    measure,
    measure!,
    occupied_locs,
    invorder!,
    partial_tr,
    select!,
    ρ,
    reorder!
export basis, ishermitian

include("utils/ast_tools.jl")

include("utils/constants.jl")
include("utils/math.jl")

include("error.jl")
include("abstract_register.jl")
include("adjoint_register.jl")

include("inspect.jl")
include("instruct.jl")

# compat with older version of dependencies
include("compat.jl")

# deprecation warns
include("deprecations.jl")

end # module


text:
# Run tests on functions in src/by_reference_doc/tracks.jl

using Spotify.Tracks

@testset verbose = true "GET-request endpoints for tracks" begin

    track_id = SpId()

    @test ~isempty(tracks_get_audio_analysis(track_id)[1])

    @test ~isempty(tracks_get_audio_features(track_id)[1])

    # Cycle through different input keywords for testing
    markets = ["US", "NL", "DE"]

    @testset "For market = $(market_id)" for market_id in markets

        @test ~isempty(tracks_get_single(track_id, market = market_id)[1])
        @test ~isempty(tracks_get_multiple("4iV5W9uYEdYUVa79Axb7Rh,1301WleyT98MSxVHPZCA6",
                                                    market = market_id)[1])
        @test ~isempty(tracks_get_saved(market = market_id)[1])

    end

    @test ~isempty(tracks_get_contains("$track_id, 4VqPOruhp5EdPBeR92t6lQ")[1])

end


text:
# defines AB_SERVER_BIN
include(joinpath(dirname(dirname(pathof(PLCTag))), "deps", "ab_server.jl"))

const TAG_TYPES = Dict(
	:SINT  => (Int8,    :int8),
	:INT   => (Int16,   :int16),
	:DINT  => (Int32,   :int32),
	:LINT  => (Int64,   :int64),
	:USINT => (UInt8,   :uint8),
	:UINT  => (UInt16,  :uint16),
	:UDINT => (UInt32,  :uint32),
	:ULINT => (UInt64,  :uint64),
	:REAL  => (Float32, :float32),
	:LREAL => (Float64, :float64),
)

const TAGS = [
	# (:SINT,  1),  # WARNING: using SINT results in read/write errors from server
	(:INT,   1),
	(:DINT,  1),
	(:LINT,  1),
	(:REAL,  1),
	(:LREAL, 1),
	# (:SINT,  4),  # WARNING: arrays seem problematic too
	# (:INT,   4),
	# (:DINT,  4),
	# (:LINT,  4),
	# (:REAL,  4),
	# (:LREAL, 4),
]

@testset "PLCTag" begin
	server = nothing
	try
		rng = MersenneTwister(123456)

		ip = get(ENV, "AB_SERVER_IP", "127.0.0.1")
		if ip == "127.0.0.1"
			server = open(Cmd([AB_SERVER_BIN, "--plc=ControlLogix", "--path=1,0", map(((t, c),) -> "--tag=tag$(t)$(c):$(t)[$(c)]", TAGS)...]))
			timedwait(() -> try ; close(connect(44818)) ; true ; catch ; false ; end, 3.0, pollint = 0.001) === :ok || error("Failed to start libplctag test server")
		end

		vals = []
		@testset "Low-level API" begin
			@test unsafe_string(plc_tag_decode_error(c"PLCTAG_STATUS_OK")) == "PLCTAG_STATUS_OK"

			tags = map(TAGS) do (t, c)
				(x, y) = TAG_TYPES[t]
				push!(vals, rand(rng, x, c))

				tag = plc_tag_create("protocol=ab-eip&gateway=$(ip)&path=1,0&cpu=compactlogix&elem_size=$(sizeof(x))&elem_count=$(c)&name=tag$(t)$(c)&debug=0", 1000)
				@test tag > 0
				@test plc_tag_status(tag) == c"PLCTAG_STATUS_OK"
				@test plc_tag_get_size(tag) == sizeof(x)*c
				@test plc_tag_get_int_attribute(tag, "size", -1) == sizeof(x)*c
				@test plc_tag_get_int_attribute(tag, "elem_size", -1) == sizeof(x)
				@test plc_tag_get_int_attribute(tag, "elem_count", -1) == c
				return (tag, c, x, y)
			end

			@testset "Sync Usage" begin
				for (ind, (tag, c, x, y)) in enumerate(tags)
					vals[ind] = rand(rng, x, c)
					for i in eachindex(vals[ind])
						@test getproperty(PLCTag, Symbol(:plc_tag_set_, y))(tag, i-1, vals[ind][i]) == c"PLCTAG_STATUS_OK"
					end
					@test plc_tag_write(tag, 1000) == c"PLCTAG_STATUS_OK"

					@test plc_tag_read(tag, 1000) == c"PLCTAG_STATUS_OK"
					for i in eachindex(vals[ind])
						@test getproperty(PLCTag, Symbol(:plc_tag_get_, y))(tag, i-1) == vals[ind][i]
					end
				end
			end

			@testset "Async Usage" begin
				for (ind, (tag, c, x, y)) in enumerate(tags)
					vals[ind] = rand(rng, x, c)
					for i in eachindex(vals[ind])
						@test getproperty(PLCTag, Symbol(:plc_tag_set_, y))(tag, i-1, vals[ind][i]) == c"PLCTAG_STATUS_OK"
					end
					@test plc_tag_write(tag, 0) in (c"PLCTAG_STATUS_OK", c"PLCTAG_STATUS_PENDING")
					@test timedwait(() -> plc_tag_status(tag) == c"PLCTAG_STATUS_OK", 1.0, pollint = 0.001) === :ok
				end

				# WARNING: truly async is unsupported by test server?
				# for (ind, (tag, c, x, y)) in enumerate(tags)
				# 	@test timedwait(() -> plc_tag_status(tag) == c"PLCTAG_STATUS_OK", 1.0, pollint = 0.001) === :ok
				# end

				# for (ind, (tag, c, x, y)) in enumerate(tags)
				# 	@test plc_tag_read(tag, 0) in (c"PLCTAG_STATUS_OK", c"PLCTAG_STATUS_PENDING")
				# end

				for (ind, (tag, c, x, y)) in enumerate(tags)
					@test plc_tag_read(tag, 0) in (c"PLCTAG_STATUS_OK", c"PLCTAG_STATUS_PENDING")
					@test timedwait(() -> plc_tag_status(tag) == c"PLCTAG_STATUS_OK", 1.0, pollint = 0.001) === :ok
					for i in eachindex(vals[ind])
						@test getproperty(PLCTag, Symbol(:plc_tag_get_, y))(tag, i-1) == vals[ind][i]
					end
				end
			end

			for (tag, c, x, y) in tags
				tag > 0 && @test plc_tag_destroy(tag) == c"PLCTAG_STATUS_OK"
			end
		end

		@testset "High-level API" begin
			plc = PLC(
				protocol = "ab-eip",
				gateway = "127.0.0.1",
				path = "1,0",
				cpu = "compactlogix",
				debug = 0,
			)

			tags = map(TAGS) do (t, c)
				(x, y) = TAG_TYPES[t]
				tag = PLCRef{x}(plc, "tag$(t)$(c)", c)
				@test length(tag) == c
				return (tag, c, x, y)
			end

			@testset "Sync Usage" begin
				for (ind, (tag, c, x, y)) in enumerate(tags)
					vals[ind] = rand(rng, x, c)

					for i in eachindex(vals[ind])
						tag[i] = vals[ind][i]
						@test tag[i] == vals[ind][i]
					end
				end
			end

			@testset "Async Usage" begin
				for (ind, (tag, c, x, y)) in enumerate(tags)
					vals[ind] = rand(rng, x, c)

					for i in eachindex(vals[ind])
						write(tag, vals[ind][i], i)
					end
					flush(tag)
				end

				# WARNING: truly async is unsupported by test server?
				# for (ind, (tag, c, x, y)) in enumerate(tags)
				# 	flush(tag)
				# end

				# for (ind, (tag, c, x, y)) in enumerate(tags)
				# 	read(tag)
				# end

				for (ind, (tag, c, x, y)) in enumerate(tags)
					read(tag)
					for i in eachindex(vals[ind])
						@test fetch(tag, i) == vals[ind][i]
					end
				end
			end
		end

	finally
		isnothing(server) || (kill(server); wait(server))
	end
end


text:
#= using Pkg =#
#= Pkg.add("Plots") =#
#= Pkg.add("CSV") =#
#= Pkg.add("DataFrames") =#
#= Pkg.add("NPZ") =#
using Plots
using CSV
using DataFrames
using NPZ

function stretch_to_500k(input_data)
    stretched = Float64[]
    for x in input_data
        for i in range(1, stop=x)
            append!(stretched, x)
        end
    end
    return stretched
end

println("Loading Data")
darlaDqn = npzread("darla_dqn_history.npz")["episode_reward"]
myVaeDqn = npzread("myvae_dqn_history.npz")["episode_reward"]
#= stylegan = npzread("stylegan_dqn_training_history_500k_again.npz")["episode_reward"] =#
#= temporalVae = npzread("vae_dqn_training_history_500k_modified.npz")["episode_reward"] =#

smoothedDarlaDqn = float(copy(darlaDqn))
smoothedmyVaeDqn = float(copy(myVaeDqn))
#= smoothedstylegan = float(copy(stylegan)) =#
#= smoothedtemporalVae = float(copy(temporalVae)) =#

#= Smooth the array =#
println("Smoothing Data")
alpha = 0.01
for i in 2:length(smoothedDarlaDqn)
    smoothedDarlaDqn[i] = alpha * darlaDqn[i] + (1. - alpha) * smoothedDarlaDqn[i-1]
end
for i in 2:length(smoothedmyVaeDqn)
    smoothedDarlaDqn[i] = alpha * darlaDqn[i] + (1. - alpha) * smoothedDarlaDqn[i-1]
end
#= for i in 2:length(smoothedstylegan) =#
#=     smoothedstylegan[i] = alpha * stylegan[i] + (1. - alpha) * smoothedstylegan[i-1] =#
#= end =#
#= for i in 2:length(smoothedtemporalVae) =#
#=     smoothedtemporalVae[i] = alpha * temporalVae[i] + (1. - alpha) * smoothedtemporalVae[i-1] =#
#= end =#

println("Stretching Data")
darlaDqn = stretch_to_500k(darlaDqn)
smoothedDarlaDqn = stretch_to_500k(smoothedDarlaDqn)
#= smoothedstylegan = stretch_to_500k(smoothedstylegan) =#
#= smoothedtemporalVae = stretch_to_500k(smoothedtemporalVae) =#

#= minLength = minimum([length(smoothedvanillaJaco) length(smoothedDarlaDqn) length(smoothedstylegan) length(smoothedtemporalVae)]) =#
minLength = minimum([length(darlaDqn) length(smoothedDarlaDqn)])
println("minLength is ", minLength)

darlaDqn = darlaDqn[1:minLength]
smoothedDarlaDqn = smoothedDarlaDqn[1:minLength]
#= smoothedstylegan = smoothedstylegan[1:minLength] =#
#= smoothedtemporalVae = smoothedtemporalVae[1:minLength] =#

println("Plotting Data")
x_data = range(1, stop=minLength)
#= plot(x_data, [smoothedvanillaJaco smoothedDarlaDqn smoothedstylegan smoothedtemporalVae], label=["Smoothed Vanilla Jaco" "Smoothed Visual Jaco" "Smoothed Stylegan" "Smoothed Temporal VAE"], xlabel="Number of Timesteps", ylabel="Reward", title="Reward vs Training Timestep", legend=:bottomright) =#
plot(x_data, [darlaDqn smoothedDarlaDqn], label=["Visual Jaco" "Smoothed Visual Jaco"], linealpha=[0.5 1], xlabel="Number of Timesteps", ylabel="Reward", title="Reward vs Training Timestep", legend=:bottomright)

savefig("losses.png")


text:
# This file was generated, do not modify it.

using CategoricalArrays

v = categorical(["AA", "BB", "CC", "AA", "BB", "CC"])

levels(v)

v = categorical([1, 2, 3, 1, 2, 3, 1, 2, 3], ordered=true)

levels(v)

v[1] < v[2]

v = categorical(["high", "med", "low", "high", "med", "low"], ordered=true)

levels(v)

v[1] < v[2]

levels!(v, ["low", "med", "high"])

v[1] < v[2]

v = categorical(["AA", "BB", missing, "AA", "BB", "CC"]);

levels(v)


text:
using MyLib

function count_to_ten()
    count = zero(Int32)
    while count < 10
        count = increment32(count)
    end
end

count_to_ten()


text:
using ComplexMixtures, PDBTools
const CM = ComplexMixtures

# PDB file of the system simulated
pdb = readPDB("../Data/system.pdb")

# Load results of a ComplexMixtures run
R = CM.load("../Data/results_glyc50.json")

# Inform which is the solute
protein = select(pdb,"protein")
solute = CM.Selection(protein,nmols=1)

# Compute the 3D density grid and output it to the PDB file
# here we use dmax=3.5 such that the the output file is not too large
grid = CM.grid3D(solute=solute,
                 solute_atoms=protein,
                 mddf_result=R,
                 output_file="grid.pdb",
                 dmin=1.5,dmax=3.5)


text:
using MultivariateFunctions
using DataFrames
using Dates
using Random
using GLM

tol = 10*eps()

Random.seed!(1)
obs = 1000
X = rand(obs)
y = X .+ rand(Normal(),obs) .+ 7

# Basic use case with 2 degrees
lm1 = fit(LinearModel, hcat(ones(obs), X, X .^ 2), y)
glm_preds = predict(lm1, hcat(ones(obs), X, X .^ 2))
package_approximation, reg1 = create_ols_approximation(y, X, 2; intercept = true, base_x = 0.0)
package_predictions = evaluate.(package_approximation, X)
sum(abs.(glm_preds .- package_predictions)) < 1e-10

# Degree of 1 with no intercept
lm1 = fit(LinearModel, hcat(X), y)
glm_preds = predict(lm1, hcat(X))
package_approximation, reg2 = create_ols_approximation(y, X, 1; intercept = false)
package_predictions = evaluate.(package_approximation, X)
sum(abs.(glm_preds .- package_predictions)) < 1e-10

# Degree of 1 with no intercept
lm1 = fit(LinearModel, hcat(ones(obs)), y)
glm_preds = predict(lm1, hcat(ones(obs)))
package_approximation, reg3 = create_ols_approximation(y, X, 0; intercept = true)
package_predictions = evaluate.(package_approximation, X)
sum(abs.(glm_preds .- package_predictions)) < 1e-10

# With applying a base
Xmin = X .- 5
lm1 = fit(LinearModel, hcat(ones(obs), Xmin, Xmin .^ 2), y)
glm_preds = predict(lm1, hcat(ones(obs), Xmin, Xmin .^ 2))
package_approximation, reg4 = create_ols_approximation(y, X, 2; intercept = true, base_x = 5.0)
package_predictions = evaluate.(package_approximation, X)
sum(abs.(glm_preds .- package_predictions)) < 1e-10

# With x as dates and a base.
X = Array{Date}(undef, obs)
baseDate =  Date(2016, 7, 21)
StartDate = Date(2018, 7, 21)
for i in 1:obs
    X[i] = StartDate +Dates.Day(2* (i-1))
end
XConverted = years_between.(X, baseDate)
lm1 = fit(LinearModel, hcat(ones(obs), XConverted, XConverted .^ 2), y)
glm_preds = predict(lm1, hcat(ones(obs), XConverted, XConverted .^ 2))
package_approximation, reg5 = create_ols_approximation(y, X, 2; intercept = true, base_date = baseDate)
package_predictions = evaluate.(package_approximation, X)
sum(abs.(glm_preds .- package_predictions)) < 1e-10

##### Multivariate Tests

Random.seed!(1992)
nObs = 1000
dd = DataFrame()
dd[:x] = rand( Normal(),nObs) + 0.1 .* rand( Normal(),nObs)
dd[:z] = rand( Normal(),nObs) + 0.1 .* rand( Normal(),nObs)
dd[:w] = (0.5 .* rand( Normal(),nObs)) .+ 0.7.*(dd[:z] .- dd[:x]) + 0.1 .* rand( Normal(),nObs)
dd[:y] = (dd[:x] .*dd[:w] ) .* (dd[:z] .- dd[:w]) .+ dd[:x] + rand( Normal(),nObs)
dd[7,:y] = 1.0

# Creating basic PE_Units
linear    = PE_Unit(0.0,0.0,1)
quadratic = PE_Unit(0.0,0.0,2)
# Creating Model
constant_term = PE_Function(1.0, Dict{Symbol,PE_Unit}())
x_lin         = PE_Function(1.0, Dict{Symbol,PE_Unit}(:x => linear))
z_lin         = PE_Function(1.0, Dict{Symbol,PE_Unit}(:z => linear))
w_lin         = PE_Function(1.0, Dict{Symbol,PE_Unit}(:w => linear))
w_quad        = PE_Function(1.0, Dict{Symbol,PE_Unit}(:w => quadratic))
x_lin_z_quad  = PE_Function(1.0, Dict{Symbol,PE_Unit}([:x, :z] .=> [linear, quadratic]))
model = constant_term + x_lin + z_lin + w_lin + w_quad + x_lin_z_quad
# Regression and testing
mod_1, reg_1 = create_ols_approximation(dd, :y, model)
dd[:predicted_y] = evaluate(mod_1, dd)
sum(abs.(dd[:predicted_y] - reg_1.rr.mu)) < 1e-12

# Now we look at a saturated model.
y = :y
x_variables = [:x,:z]
degree = 2
intercept = true
univariate_dim_name = :default
mod_2, reg_2 = create_saturated_ols_approximation(dd, y, x_variables, degree)
dd[:predicted_y_2] = evaluate(mod_2, dd)
sum(abs.(dd[:predicted_y_2] - reg_2.rr.mu)) < 1e-12

# Now we look at a saturated model.
y = :y
x_variables = [:x,:y,:z]
degree = 3
intercept = true
univariate_dim_name = :default
mod_3, reg_3 = create_saturated_ols_approximation(dd, y, x_variables, degree)
dd[:predicted_y_3] = evaluate(mod_3, dd)
sum(abs.(dd[:predicted_y_3] - reg_3.rr.mu)) < 1e-11

# Now the second model should have higher RSS than the first because it has fewer variables and terms. We can test this:
sum((dd[:predicted_y_2] .- dd[:y]).^2) > sum((dd[:predicted_y] .- dd[:y]).^2)
# The third model should have lower RSS than the first because it has more terms. Testing this:
sum((dd[:predicted_y_3] .- dd[:y]).^2) < sum((dd[:predicted_y] .- dd[:y]).^2)


text:
#= Including this file

julia> include("code/julia/cmocean4plotly.jl")

generates file "ColormapsOcean.jl" which then has the code for
the Julia module ColormapsOcean. It makes available the ocean colormaps

https://matplotlib.org/cmocean/

with the PlotlyJS.jl plotting package. This file is included
in the ColormapsOcean.jl package, but it is only needed if the original
color maps get updated or enhanced, or the Julia code has to be regenerated
for some reason. Please edit then paths according to your local setup.

=#
ppath = "code/cmocean/cmocean/rgb"
fid = open(joinpath("code/julia", "ColormapsOcean.jl"), "w")
@printf(fid, "__precompile__()\n\nmodule ColormapsOcean\n\n")
@printf(fid, "export cmocean\n\n")
@printf(fid, "const cmocean = Dict{Symbol, Vector{Tuple{Float64,String}}}()\n\n")
for (_, __, files) in walkdir(ppath)
    for fnm in files
        if !ismatch(r"txt$", fnm)
            continue
        end
        @printf(fid, "cmocean[:%s] = [\n", fnm[1:search(fnm, '-')-1])
        lns = open(readdlm, joinpath(ppath, fnm))
        for k in 1:size(lns, 1)
            rgb = view(lns, k, :)
            if k>1
                println(fid, ",")
            end
            print(fid, "($((k-1)/255), \"rgb($(round(UInt8, rgb[1]*255)),$(round(UInt8, rgb[2]*255)),$(round(UInt8, rgb[3]*255)))\")")
        end
        println(fid, "]\n")
    end
end
println(fid, "end   # module ColormapsOcean")
close(fid)


text:
module BarFrameworks

# Write your package code here.

using HomotopyContinuation

include("myStructs.jl")
include("myFunctions.jl")

end


text:
using ..Phylo
using Diversity
using Diversity.API
using Statistics
using AxisArrays

abstract type AbstractPhyloTypes{Tree <: AbstractTree} <:
    Diversity.API.AbstractTypes
end

import Diversity.API: _addedoutputcols
function _addedoutputcols(::AbstractPhyloTypes{TS}) where
    {LABEL, NL, BL, TS <: TreeSet{LABEL, NL, BL, <: AbstractTree}}
    Dict{Symbol, Type}(:treename => LABEL)
end

import Diversity.API: _getaddedoutput
function _getaddedoutput(pt::AbstractPhyloTypes{TS}) where
    {LABEL, NL, BL, TS <: TreeSet{LABEL, NL, BL, <: AbstractTree}}
    Dict{Symbol, LABEL}(:treename => first(treenameiter(pt.tree)))
end

struct PhyloBranches{Tree} <: AbstractPhyloTypes{Tree}
    tree::Tree
    nleaf::Int64
    nancestral::Int64
    leafnames::Vector{String}
    ancestralnames::Vector{String}
    ancestralmatrix::Matrix{Float64}
    Zmatrix::Matrix{Float64}
end

function PhyloBranches(tree::Tree) where Tree <: AbstractTree
    leafnames = [getnodename(tree, node)
                 for node in traversal(tree, preorder) if isleaf(tree, node)]
    nleaf = length(leafnames)
    nleaf > 0 || error("Too few species")
    leafinfo = Dict{String, Float64}()
    speciesinfo = Dict{String, Tuple{String, String, Float64}}()
    ancestralnames = String[]

    for leaf in leafnames
        branches = branchhistory(tree, leaf)
        leafinfo["$leaf"] = heighttoroot(tree, leaf)
        for branch in branches
            branchname = getbranchname(tree, branch)
            name = "$leaf : $branchname"
            push!(ancestralnames, name)
            speciesinfo[name] =
                tuple("$leaf", "$branchname",
                      getlength(tree, branch) / leafinfo["$leaf"])
        end
    end

    nancestral = length(ancestralnames)
    Lbar = Statistics.mean(collect(values(leafinfo)))
    ancestralmatrix = Matrix{Float64}(undef, nancestral, nleaf)
    fill!(ancestralmatrix, 0.0)
    for i in 1:nancestral
        for j in 1:nleaf
            if speciesinfo[ancestralnames[i]][1] == leafnames[j]
                ancestralmatrix[i, j] =
                    speciesinfo[ancestralnames[i]][3] * leafinfo[leafnames[j]]
            end
        end
    end

    Zmatrix = Matrix{Float64}(undef, nancestral, nancestral)
    fill!(Zmatrix, 0.0)
    for i in 1:nancestral
        for j in 1:nancestral
            if haskey(speciesinfo,
                      "$(speciesinfo[ancestralnames[j]][1]) : " *
                      "$(speciesinfo[ancestralnames[i]][2])")
                Zmatrix[i, j] = 1.0 /
                    leafinfo[speciesinfo[ancestralnames[j]][1]]
            end
        end
    end

    return PhyloBranches{Tree}(tree, nleaf, nancestral, leafnames,
                               ancestralnames, ancestralmatrix, Zmatrix)
end

function PhyloBranches(treeset::TS) where TS <: TreeSet
    ntrees(treeset) == 1 ||
        error("Can currently only handle one tree in a PhyloSet")
    pt = PhyloBranches(first(treeset))
    return PhyloBranches{TS}(treeset, pt.nleaf, pt.nancestral,
                             pt.leafnames, pt.ancestralnames,
                             pt.ancestralmatrix, pt.Zmatrix)
end

import Diversity.API: _gettypenames
function _gettypenames(phy::PhyloBranches, raw::Bool)
    return raw ? phy.leafnames : phy.ancestralnames
end

import Diversity.API: _counttypes
function _counttypes(phy::PhyloBranches, raw::Bool)
    return raw ? phy.nleaf : phy.nancestral
end

import Diversity.API: _calcabundance
function _calcabundance(phy::PhyloBranches,
                        raw::AbstractMatrix{<: AbstractFloat})
    processed = phy.ancestralmatrix * raw
    scale = sum(processed)
    processed ./= scale
    return processed, scale
end

import Diversity.API: _calcsimilarity
function _calcsimilarity(phy::PhyloBranches, scale::Real)
    return phy.Zmatrix .* scale
end

import Diversity.API: _calcordinariness
function _calcordinariness(phy::PhyloBranches,
                           processed::AbstractMatrix{<: AbstractFloat},
                           scale::Real)
    return (phy.Zmatrix * processed) .* scale
end

import Diversity.API: _getdiversityname
_getdiversityname(::PhyloBranches) = "Phylogenetic Branch"


text:
# Initial tests to see if deleteLeaf works
# as expected with QuartetNetwork type
# Claudia November 2014

include("../examples/bad_triangle_example.jl")
include("../examples/case_f_example.jl")
qnet = QuartetNetwork(net);
printEdges(qnet)
printNodes(qnet)
printEdges(net)
printNodes(net)
qnet.hasEdge

# bad triangle
deleteLeaf!(qnet,qnet.node[8])
deleteLeaf!(qnet,qnet.node[7])
deleteLeaf!(qnet,qnet.node[5])
deleteLeaf!(qnet,qnet.node[6])
deleteLeaf!(qnet,qnet.node[1])

# bad diamond
deleteLeaf!(qnet,qnet.node[10])
deleteLeaf!(qnet,qnet.node[8])
deleteLeaf!(qnet,qnet.node[6])
deleteLeaf!(qnet,qnet.node[7])
deleteLeaf!(qnet,qnet.node[4])

# Case G : no simplification in deleteLeaf!
include("../examples/case_g_example.jl")

qnet = QuartetNetwork(net);
printEdges(qnet)
printNodes(qnet)

deleteLeaf!(qnet,qnet.node[10]) # 8
deleteLeaf!(qnet,qnet.node[7]) # 7
deleteLeaf!(qnet,qnet.node[4]) # 4
deleteLeaf!(qnet,qnet.node[2]) # 2
deleteLeaf!(qnet,qnet.node[1]) # 1

# Case C: bad triangle II
include("../examples/case_c_example.jl")

qnet = QuartetNetwork(net);
printEdges(qnet)
printNodes(qnet)

deleteLeaf!(qnet,qnet.node[1]) # 1
deleteLeaf!(qnet,qnet.node[2]) # 2
deleteLeaf!(qnet,qnet.node[8]) # 5
deleteLeaf!(qnet,qnet.node[9]) # 6
deleteLeaf!(qnet,qnet.node[4]) # 3

# test extractQuartet -------------------

include("../examples/bad_triangle_example.jl")
q1 = Quartet(1,["6","7","1","8"],[0.5,0.4,0.1]);
qnet = extractQuartet(net,q1);
printEdges(qnet)
printNodes(qnet)
qnet.hasEdge

include("../examples/bad_triangle_example.jl")
qnet = extractQuartet(net,net.node[6],net.node[7],net.node[1],net.node[8]);
printEdges(qnet)
printNodes(qnet)
qnet.hasEdge

include("../examples/case_f_example.jl")

# bad triangle
qnet=QuartetNetwork(net);
printEdges(qnet)
printNodes(qnet)
deleteLeaf!(qnet,qnet.node[8])
deleteLeaf!(qnet,qnet.node[7])
deleteLeaf!(qnet,qnet.node[5])
deleteLeaf!(qnet,qnet.node[6])
deleteLeaf!(qnet,qnet.node[1])

# bad diamond
deleteLeaf!(qnet,qnet.node[10])
deleteLeaf!(qnet,qnet.node[8])


text:
using JSON
Nproc=4

rmprocs(workers())
addprocs(Nproc)

println("workers: ",length(workers()))

@everywhere using Base.Test
@everywhere import Base.convert
@everywhere using MAT
@everywhere using ForwardDiff
@everywhere using Polynomials
@everywhere import ForwardDiff.DiffBase
@everywhere using Optim
@everywhere gc(true)

# cluster idibaps
@everywhere path="/home/emma/github/fit_behaviour_julia/"
# cluster
# @everywhere path="/home/hcli64/hcli64751/rats/fit_behaviour_julia/"
#
@everywhere include(path*"/functions/functions_julia.jl")
@everywhere include(path*"/functions/parameter_functions.jl")
@everywhere include(path*"/functions/utils.jl")

rat_num = 4
path_load=path*"scripts_rats/"
filename="/data/processed_data.json"

data = Dict()
open(path_load*filename, "r") do f
    global data
    dicttxt = readstring(f)  # file information to string
    data=JSON.parse(dicttxt)  # parse and transform data
end

stimulus = data[1][rat_num]             # Data frames obtained from a gaussian distribution of a determined coherence
# coherences = data[2][rat_num]           # Target coherences used for each trial
# rewards = data[3][rat_num]              # Correct choice side
decisions = data[4][rat_num]            # Actual decision made by the rat each trial
# performances = data[5][rat_num]         # 0 if the rat chose the correct choice, 1 if the rat chose the incorrect choice
# target_sigmas = data[6][rat_num]        # Target sigmas used for each trial

stim_10 = []
stim_20 = []
choices_10 = []
choices_20 = []

# Separar els estímuls de 10 frames dels de 20 frames
for session in 1:size(stimulus)[1]
    for trial in 1:size(stimulus[session])[1]
        stim_1 = stimulus[session][trial]
        if size(stim_1)[1] == 10
            stim_10 = [stim_10;stim_1]
            choices_10 = [choices_10;decisions[session][trial]]
        end
        if size(stim_1)[1] == 20
            stim_20 = [stim_20;stim_1]
            choices_20 = [choices_20;decisions[session][trial]]
        end
    end
end
n_stim_10 = convert(Int64, size(stim_10)[1]/10)
n_stim_20 = convert(Int64, size(stim_20)[1]/20)

# if n_stim_10 != 0
#     stim_10 = permutedims(reshape(stim_10,(10,n_stim_10)))
# end
# if n_stim_20 != 0
#     stim_20 = permutedims(reshape(stim_20,(20,n_stim_20)))
# end

results = Dict("stim_10"=>stim_10,"stim_20"=>stim_20,"choices_10"=>choices_10,"choices_20"=>choices_20,"n_stim_10"=>n_stim_10,"n_stim_20"=>n_stim_20)

filename_save = path_load*"data/processed_data_rat"*string(rat_num)*".json"
# pass data as a json string (how it shall be displayed in a file)
stringdata = JSON.json(results)

# write the file with the stringdata variable information
open(filename_save, "w") do f
        write(f, stringdata)
     end

max_stars_repo_name: emmafdz/FinalDegreeProject

max_issues_repo_name: emmafdz/FinalDegreeProject

max_forks_repo_name: emmafdz/FinalDegreeProject