text:
using GeometricBase.Config
using GeometricBase.Utils
using GeometricIntegrators.Integrators
using GeometricIntegrators.Integrators.SPARK
using GeometricIntegrators.Tableaus
using GeometricProblems.LotkaVolterra2d
using SimpleSolvers
using Test

SimpleSolvers.set_config(:nls_atol, 8eps())
SimpleSolvers.set_config(:nls_rtol, 2eps())

SimpleSolvers.set_config(:nls_atol_break, Inf)
SimpleSolvers.set_config(:nls_rtol_break, Inf)
SimpleSolvers.set_config(:nls_stol_break, Inf)

const Δt = 0.01
const nt = 10
const q₀ = [1.0, 1.0]
const parameters = (a₁=1.0, a₂=1.0, b₁=-1.0, b₂=-2.0)

ode  = lotka_volterra_2d_ode(q₀; params=parameters)
idae = lotka_volterra_2d_idae(q₀; params=parameters)
pdae = lotka_volterra_2d_pdae(q₀; params=parameters)
ldae = lotka_volterra_2d_ldae(q₀; params=parameters)
ldae_slrk = lotka_volterra_2d_slrk(q₀; params=parameters)

int  = IntegratorFIRK(ode, TableauGauss(8), Δt)
sol  = integrate(ode, int, nt)
refx = sol.q[end]

@testset "$(rpad("SLRK integrators",80))" begin

    int = IntegratorSLRK(ldae_slrk, TableauSLRKLobattoIIIAB(2), Δt)
    sol = integrate(ldae_slrk, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = IntegratorSLRK(ldae_slrk, TableauSLRKLobattoIIIAB(3), Δt)
    sol = integrate(ldae_slrk, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-11

    int = IntegratorSLRK(ldae_slrk, TableauSLRKLobattoIIIAB(4), Δt)
    sol = integrate(ldae_slrk, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-15

    int = IntegratorSLRK(ldae_slrk, TableauSLRKLobattoIIIBA(2), Δt)
    sol = integrate(ldae_slrk, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = IntegratorSLRK(ldae_slrk, TableauSLRKLobattoIIIBA(3), Δt)
    sol = integrate(ldae_slrk, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-11

    int = IntegratorSLRK(ldae_slrk, TableauSLRKLobattoIIIBA(4), Δt)
    sol = integrate(ldae_slrk, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-15

    int = IntegratorSLRK(ldae_slrk, TableauSLRKLobattoIIICC̄(2), Δt)
    sol = integrate(ldae_slrk, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = IntegratorSLRK(ldae_slrk, TableauSLRKLobattoIIICC̄(3), Δt)
    sol = integrate(ldae_slrk, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-11

    int = IntegratorSLRK(ldae_slrk, TableauSLRKLobattoIIICC̄(4), Δt)
    sol = integrate(ldae_slrk, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-15

    int = IntegratorSLRK(ldae_slrk, TableauSLRKLobattoIIIC̄C(2), Δt)
    sol = integrate(ldae_slrk, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = IntegratorSLRK(ldae_slrk, TableauSLRKLobattoIIIC̄C(3), Δt)
    sol = integrate(ldae_slrk, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-11

    int = IntegratorSLRK(ldae_slrk, TableauSLRKLobattoIIIC̄C(4), Δt)
    sol = integrate(ldae_slrk, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-15

    int = IntegratorSLRK(ldae_slrk, TableauSLRKLobattoIIID(2), Δt)
    sol = integrate(ldae_slrk, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = IntegratorSLRK(ldae_slrk, TableauSLRKLobattoIIID(3), Δt)
    sol = integrate(ldae_slrk, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-11

    int = IntegratorSLRK(ldae_slrk, TableauSLRKLobattoIIID(4), Δt)
    sol = integrate(ldae_slrk, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-15

    int = IntegratorSLRK(ldae_slrk, TableauSLRKLobattoIIIE(2), Δt)
    sol = integrate(ldae_slrk, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-6

    int = IntegratorSLRK(ldae_slrk, TableauSLRKLobattoIIIE(3), Δt)
    sol = integrate(ldae_slrk, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = IntegratorSLRK(ldae_slrk, TableauSLRKLobattoIIIE(4), Δt)
    sol = integrate(ldae_slrk, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-15

end

@testset "$(rpad("SPARK integrators",80))" begin

    int = Integrator(idae, TableauSPARKGLRK(1), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-6

    int = Integrator(idae, TableauSPARKGLRK(2), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = Integrator(idae, TableauSPARKLobattoIIIAIIIB(3), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-6

    int = Integrator(idae, TableauSPARKLobattoIIIAIIIB(4), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-10

    int = Integrator(idae, TableauSPARKGLRKLobattoIIIAIIIB(1), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 5E-4

    int = Integrator(idae, TableauSPARKGLRKLobattoIIIAIIIB(2), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 3E-4

    int = Integrator(idae, TableauSPARKGLRKLobattoIIIAIIIB(3), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-4

    int = Integrator(idae, TableauSPARKGLRKLobattoIIIBIIIA(1), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 5E-4

    int = Integrator(idae, TableauSPARKGLRKLobattoIIIBIIIA(2), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 3E-4

    int = Integrator(idae, TableauSPARKGLRKLobattoIIIBIIIA(3), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-4

end

@testset "$(rpad("VPARK integrators",80))" begin

    int = Integrator(idae, TableauSymplecticProjection(:glrk1ps, TableauGauss(1), TableauGauss(1)), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-6

    int = Integrator(idae, TableauSymplecticProjection(:glrk2ps, TableauGauss(2), TableauGauss(2)), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = Integrator(idae, TableauGausspSymplectic(1), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-6

    int = Integrator(idae, TableauGausspSymplectic(2), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = Integrator(idae, TableauLobattoIIIAIIIBpSymplectic(2), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-6

    int = Integrator(idae, TableauLobattoIIIAIIIBpSymplectic(3), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-11

    int = Integrator(idae, TableauLobattoIIIAIIIBpSymplectic(4), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-15

    int = Integrator(idae, TableauLobattoIIIBIIIApSymplectic(2), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-6

    int = Integrator(idae, TableauLobattoIIIBIIIApSymplectic(3), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-11

    int = Integrator(idae, TableauLobattoIIIBIIIApSymplectic(4), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-16

end

@testset "$(rpad("VSPARK integrators",80))" begin

    ### VSPARK Integrators ###

    int = IntegratorVSPARK(idae, TableauSPARKGLRKLobattoIIIAIIIB(1), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-6

    int = IntegratorVSPARK(idae, TableauSPARKGLRKLobattoIIIAIIIB(2), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = IntegratorVSPARK(idae, TableauSPARKGLRKLobattoIIIAIIIB(3), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-15

    int = IntegratorVSPARK(idae, TableauSPARKGLRKLobattoIIIBIIIA(1), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-6

    int = IntegratorVSPARK(idae, TableauSPARKGLRKLobattoIIIBIIIA(2), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = IntegratorVSPARK(idae, TableauSPARKGLRKLobattoIIIBIIIA(3), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-15

    ### VSPARKprimary Integrators ###

    int = Integrator(idae, TableauVSPARKGLRKpMidpoint(1), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-6

    int = Integrator(idae, TableauVSPARKGLRKpMidpoint(2), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = Integrator(idae, TableauVSPARKGLRKpSymplectic(1), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-6

    int = Integrator(idae, TableauVSPARKGLRKpSymplectic(2), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = Integrator(idae, TableauVSPARKGLRKpSymmetric(1), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-6

    int = Integrator(idae, TableauVSPARKGLRKpSymmetric(2), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = IntegratorVSPARKprimary(idae, TableauVSPARKLobattoIIIAIIIBpSymmetric(2), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-6

    int = IntegratorVSPARKprimary(idae, TableauVSPARKLobattoIIIAIIIBpSymmetric(3), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 5E-11

    int = IntegratorVSPARKprimary(idae, TableauVSPARKLobattoIIIAIIIBpSymmetric(4), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-15

    int = IntegratorVSPARKprimary(idae, TableauVSPARKLobattoIIIBIIIApSymmetric(2), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-6

    # int = IntegratorVSPARKprimary(idae, TableauVSPARKLobattoIIIBIIIApSymmetric(3), Δt)
    # sol = integrate(idae, int, nt)
    # @test relative_maximum_error(sol.q, refx) < 5E-11

    int = IntegratorVSPARKprimary(idae, TableauVSPARKLobattoIIIBIIIApSymmetric(4), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-15

    ### VSPARKsecondary Integrators ###

    int = Integrator(ldae, TableauVSPARKLobattoIIIAB(2), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = Integrator(ldae, TableauVSPARKLobattoIIIAB(3), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-11

    int = Integrator(ldae, TableauVSPARKLobattoIIIAB(4), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-15

    int = Integrator(ldae, TableauVSPARKLobattoIIIBA(2), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = Integrator(ldae, TableauVSPARKLobattoIIIBA(3), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-11

    int = Integrator(ldae, TableauVSPARKLobattoIIIBA(4), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-15

    int = Integrator(ldae, TableauVSPARKLobattoIIICC̄(2), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = Integrator(ldae, TableauVSPARKLobattoIIICC̄(3), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-11

    int = Integrator(ldae, TableauVSPARKLobattoIIICC̄(4), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-15

    int = Integrator(ldae, TableauVSPARKLobattoIIIC̄C(2), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = Integrator(ldae, TableauVSPARKLobattoIIIC̄C(3), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-11

    int = Integrator(ldae, TableauVSPARKLobattoIIIC̄C(4), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-15

    int = Integrator(ldae, TableauVSPARKLobattoIIID(2), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = Integrator(ldae, TableauVSPARKLobattoIIID(3), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-11

    int = Integrator(ldae, TableauVSPARKLobattoIIID(4), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-15

    int = Integrator(ldae, TableauVSPARKLobattoIIIE(2), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-6

    int = Integrator(ldae, TableauVSPARKLobattoIIIE(3), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = Integrator(ldae, TableauVSPARKLobattoIIIE(4), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-15

    int = Integrator(ldae, TableauVSPARKGLRKLobattoIIIAB(1), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = Integrator(ldae, TableauVSPARKGLRKLobattoIIIAB(2), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = Integrator(ldae, TableauVSPARKGLRKLobattoIIIAB(3), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-15

    int = Integrator(ldae, TableauVSPARKGLRKLobattoIIIBA(1), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = Integrator(ldae, TableauVSPARKGLRKLobattoIIIBA(2), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = Integrator(ldae, TableauVSPARKGLRKLobattoIIIBA(3), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-15

    int = Integrator(ldae, TableauVSPARKGLRKLobattoIIICC̄(1), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = Integrator(ldae, TableauVSPARKGLRKLobattoIIICC̄(2), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = Integrator(ldae, TableauVSPARKGLRKLobattoIIICC̄(3), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-15

    int = Integrator(ldae, TableauVSPARKGLRKLobattoIIIC̄C(1), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = Integrator(ldae, TableauVSPARKGLRKLobattoIIIC̄C(2), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = Integrator(ldae, TableauVSPARKGLRKLobattoIIIC̄C(3), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-15

    int = Integrator(ldae, TableauVSPARKGLRKLobattoIIID(1), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = Integrator(ldae, TableauVSPARKGLRKLobattoIIID(2), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = Integrator(ldae, TableauVSPARKGLRKLobattoIIID(3), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-15

    int = Integrator(ldae, TableauVSPARKGLRKLobattoIIIE(1), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-6

    int = Integrator(ldae, TableauVSPARKGLRKLobattoIIIE(2), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = Integrator(ldae, TableauVSPARKGLRKLobattoIIIE(3), Δt)
    sol = integrate(ldae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-15
end

@testset "$(rpad("HPARK integrators",80))" begin

    int = Integrator(pdae, TableauHPARKGLRK(1), Δt)
    sol = integrate(pdae, int, nt)
    # println(relative_maximum_error(sol.q, refx))
    @test relative_maximum_error(sol.q, refx) < 2E-6

    int = Integrator(pdae, TableauHPARKGLRK(2), Δt)
    sol = integrate(pdae, int, nt)
    # println(relative_maximum_error(sol.q, refx))
    @test relative_maximum_error(sol.q, refx) < 8E-7

    # int = Integrator(pdae, TableauHPARKLobattoIIIAIIIB(2), Δt)
    # sol = integrate(pdae, int, nt)
    # TODO
    # println(relative_maximum_error(sol.q, refx))
    # @test relative_maximum_error(sol.q, refx) < 2E-2

    int = Integrator(pdae, TableauHPARKLobattoIIIAIIIB(3), Δt)
    sol = integrate(pdae, int, nt)
    # TODO
    # println(relative_maximum_error(sol.q, refx))
    @test relative_maximum_error(sol.q, refx) < 8E-2

    int = Integrator(pdae, TableauHPARKLobattoIIIAIIIB(4), Δt)
    sol = integrate(pdae, int, nt)
    # TODO
    # println(relative_maximum_error(sol.q, refx))
    @test relative_maximum_error(sol.q, refx) < 8E-4

    # int = Integrator(pdae, TableauHPARKLobattoIIIBIIIA(2), Δt)
    # sol = integrate(pdae, int, nt)
    # TODO
    # println(relative_maximum_error(sol.q, refx))
    # @test relative_maximum_error(sol.q, refx) < 2E-2

    int = Integrator(pdae, TableauHPARKLobattoIIIBIIIA(3), Δt)
    sol = integrate(pdae, int, nt)
    # TODO
    # println(relative_maximum_error(sol.q, refx))
    @test relative_maximum_error(sol.q, refx) < 4E-3

    int = Integrator(pdae, TableauHPARKLobattoIIIBIIIA(4), Δt)
    sol = integrate(pdae, int, nt)
    # TODO
    # println(relative_maximum_error(sol.q, refx))
    @test relative_maximum_error(sol.q, refx) < 2E-2

end

@testset "$(rpad("HSPARK integrators",80))" begin
    ### HSPARK Integrators ###

    int = IntegratorHSPARK(pdae, TableauSPARKGLRK(1), Δt)
    sol = integrate(pdae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-6

    int = IntegratorHSPARK(pdae, TableauSPARKGLRK(2), Δt)
    sol = integrate(pdae, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    ### HSPARKprimary Integrators ###

    int = Integrator(pdae, TableauHSPARKGLRKpSymmetric(1), Δt)
    sol = integrate(pdae, int, nt)
    # TODO
    # println(relative_maximum_error(sol.q, refx))
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = Integrator(pdae, TableauHSPARKGLRKpSymmetric(2), Δt)
    sol = integrate(pdae, int, nt)
    # TODO
    # println(relative_maximum_error(sol.q, refx))
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = Integrator(pdae, TableauHSPARKLobattoIIIAIIIBpSymmetric(2), Δt)
    sol = integrate(pdae, int, nt)
    # TODO
    # println(relative_maximum_error(sol.q, refx))
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = Integrator(pdae, TableauHSPARKLobattoIIIAIIIBpSymmetric(3), Δt)
    sol = integrate(pdae, int, nt)
    # TODO
    # println(relative_maximum_error(sol.q, refx))
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = Integrator(pdae, TableauHSPARKLobattoIIIAIIIBpSymmetric(4), Δt)
    sol = integrate(pdae, int, nt)
    # TODO
    # println(relative_maximum_error(sol.q, refx))
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = Integrator(pdae, TableauHSPARKLobattoIIIBIIIApSymmetric(2), Δt)
    sol = integrate(pdae, int, nt)
    # TODO
    # println(relative_maximum_error(sol.q, refx))
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = Integrator(pdae, TableauHSPARKLobattoIIIBIIIApSymmetric(3), Δt)
    sol = integrate(pdae, int, nt)
    # TODO
    # println(relative_maximum_error(sol.q, refx))
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = Integrator(pdae, TableauHSPARKLobattoIIIBIIIApSymmetric(4), Δt)
    sol = integrate(pdae, int, nt)
    # TODO
    # println(relative_maximum_error(sol.q, refx))
    @test relative_maximum_error(sol.q, refx) < 4E-6

    ### HSPARKsecondary Integrators ###

    # int = IntegratorHSPARKsecondary(ldae, TableauHSPARKLobattoIIIAB(2), Δt)
    # sol = integrate(ldae, int, nt)
    # println(relative_maximum_error(sol.q, refx))
    # @test relative_maximum_error(sol.q, refx) < 4E-6

    # int = IntegratorHSPARKsecondary(ldae, TableauHSPARKLobattoIIIAB(3), Δt)
    # sol = integrate(ldae, int, nt)
    # println(relative_maximum_error(sol.q, refx))
    # @test relative_maximum_error(sol.q, refx) < 2E-11

    # int = IntegratorHSPARKsecondary(ldae, TableauHSPARKLobattoIIIAB(4), Δt)
    # sol = integrate(ldae, int, nt)
    # println(relative_maximum_error(sol.q, refx))
    # @test relative_maximum_error(sol.q, refx) < 1E-15

end


text:
#--------------Replicating Gertler and Karadi 2015-----------------------

#---------------------Load Packages--------------------------------------
using VectorAutoregressions
using Plots
plotly()
using Random, LinearAlgebra
using Statistics, GrowableArrays
using DelimitedFiles: readdlm

#-----------Set base-path------------------------------------------
path = joinpath(dirname(pathof(VectorAutoregressions)), "..")

#--------------------Load Data-------------------------------------------
X,names = readdlm(joinpath(path,"test","gertler_karadi_var_data.csv"),',',header=true)

#-------------------Setting-up the variables-----------------------------
Y = convert(Array{Float64},X[:,3:6])
Z = convert(Array{Float64},X[findlast(λ-> λ =="",X[:,7])+1:end,8:8])

#------------------Hyperparameters---------------------------------------
α = [.05, .10, .15, .25] # confidence interval
intercept = true
K = size(Y,2)
nrep = 1000
H = 50
p = 12

#---------------Cholesky VAR---------------------------------------------
V = VAR(Y,p,intercept)
mIRFa = IRFs_a(V,H,intercept)

pIRF_asy = plot(layout = grid(K,1));
let counter = 0
    for i in collect(3:K:15)
        counter +=1
        plot!(pIRF_asy, [mIRFa.IRF[i,:] mIRFa.CI.CIl[i,:] mIRFa.CI.CIh[i,:]], color = ["red" "red" "red"],
        line = [:solid :dash :dash], legend = false, label = ["cholesky"], subplot = counter)
    end
    gui(pIRF_asy)
end
#-------------External instrument VAR-------------------------------------
Y_z = Y[:,[3, 1, 2, 4]] # fedfund first
V = VAR(Y_z,p,intercept)

IRF = IRFs_ext_instrument(V, Z, H, nrep, α, intercept)

pIRF = plot(layout = grid(4,1));
for i in [3 1 2 4]
    plot!(pIRF,[IRF.IRF[:,i] IRF.CI.CIl[4,:,i] IRF.CI.CIh[4,:,i]], subplot = i, legend = false, color = ["blue" "blue" "blue"], line = [:solid :dash :dash])
end
gui(pIRF)

ext: jl

max_stars_repo_name: lucabrugnolini/VAR

max_issues_repo_name: lucabrugnolini/VAR

max_forks_repo_name: lucabrugnolini/VAR


text:
# HelperRoutines.jl

function extract_options(args, ndof::Int64 = 1)
    o = Dict{Symbol, Any}(args)
    if haskey(o, :reltol)
        reltol = o[:reltol]
    else
        reltol = 1e-8
    end
    if haskey(o, :abstol)
        abstol = o[:abstol]
    else
        abstol = 1e-6
    end
    if haskey(o, :dtmax)
        dtmax = o[:dtmax]
    else
        dtmax = 1e-3
    end
    if haskey(o, :lower)
        lb = o[:lower]
    else
        lb = -Inf * ones(ndof)
    end
    if haskey(o, :upper)
        ub = o[:upper]
    else
        ub = Inf * ones(ndof)
    end
    return (
        reltol = reltol,
        abstol = abstol,
        dtmax = dtmax,
        lower = lb,
        upper = ub
    )
end


text:
"""
    lagrange(xⁱ, m, x)

Calculate the Lagrange interpolating polynomial Lⁱₘ(x), given the roots
`xⁱ`.

Lⁱₘ(x) = ∏(j≠m) (x-xⁱⱼ)/(xⁱₘ-xⁱⱼ)

Eq. (8) Rescigno2000
"""
function lagrange(xⁱ::AbstractVector, m::Integer, x::AbstractVector)
    Lₘ = fill(1.0, length(x))
    for j in eachindex(xⁱ)
        j == m && continue
        Lₘ .*= (x .- xⁱ[j])/(xⁱ[m]-xⁱ[j])
    end
    Lₘ
end
lagrange(xⁱ::AbstractVector, m::Integer, x::Number) =
    lagrange(xⁱ, m, [x])[1]

"""
    δ(a,b)

Kronecker δ function"""
δ(a,b) = a == b ? 1 : 0

"""
    lagrangeder!(xⁱ, m, L′)

Calculate the derivative of the Lagrange interpolating polynomial
Lⁱₘ(x), given the roots `xⁱ`, *at* the roots, and storing the result
in `L′`.

Lⁱₘ′(x) = (xⁱₘ-xⁱₘ,)⁻¹ ∏(k≠m,m′) (xⁱₘ,-xⁱₖ)/(xⁱₘ-xⁱₖ), m≠m′,
          [δ(m,n) - δ(m,1)]/2wⁱₘ,                    m=m′

Eq. (20) Rescigno2000
"""
function lagrangeder!(xⁱ::AbstractVector, wⁱ::AbstractVector,
                      L′::AbstractMatrix)
    n = length(xⁱ)
    for m in 1:n
        L′[m,m] = (δ(m,n)-δ(m,1))/2wⁱ[m]
        for m′ in 1:n
            m′ == m && continue
            f = 1
            for k = 1:n
                k in [m,m′] && continue
                f *= (xⁱ[m′]-xⁱ[k])/(xⁱ[m]-xⁱ[k])
            end
            L′[m,m′] = f/(xⁱ[m]-xⁱ[m′])
        end
    end
end

function lagrangeder(grid::Grid)
    n = order(grid)

    L′ = zeros(elcount(grid), n, n)
    for i = elems(grid)
        sel = (1:n) .+ (i-1)*(n-1)
        lagrangeder!(grid.X[i,:], grid.W[i,:],
                     view(L′, i, :, :))
    end
    L′
end


text:
using Dates
using GitHub
using Pkg
using Printf
using RegistryCI
using Test
using TimeZones

const AutoMerge = RegistryCI.AutoMerge

@testset "RegistryCI.jl" begin
    @testset "RegistryCI.jl unit tests" begin
        @info("Running the RegistryCI.jl unit tests")
        include("registryci.jl")
    end

    @testset "AutoMerge.jl unit tests" begin
        @info("Running the AutoMerge.jl unit tests")
        include("automerge-unit.jl")
    end

    AUTOMERGE_RUN_INTEGRATION_TESTS = get(ENV, "AUTOMERGE_RUN_INTEGRATION_TESTS", "")::String
    if AUTOMERGE_RUN_INTEGRATION_TESTS == "true"
        @testset "AutoMerge.jl integration tests" begin
            @info("Running the AutoMerge.jl integration tests")
            include("automerge-integration.jl")
        end
    end
end


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script to compile data for some descriptive statistics (VALUES)
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

using DataFrames, CSV, XLSX, LinearAlgebra, Statistics, StatsBase, StatsPlots, Dates, Pipe

# other scripts
dir_home = "x:/VIVES/1-Personal/Florian/git/2022_01_comext/src/"
include(dir_home * "functions.jl")

# location of data input/output (io)
dir_io = "C:/Users/u0148308/data/comext/"
dir_dropbox = "C:/Users/u0148308/Dropbox/BREXIT/"

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Raw data
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

df_VLAIO = CSV.read(dir_dropbox * "rawdata/" * "df_VLAIO" * ".csv", DataFrame)
transform!(df_VLAIO, ["PRODUCT_NC", "PRODUCT_NC_LAB1", "PRODUCT_NC_LAB2", "DECLARANT_ISO", "TRADE_TYPE", "PARTNER_ISO", "FLOW"] .=> ByRow(string), renamecols=false)
transform!(df_VLAIO, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> convert(Union{Missing, Float64}, x)), renamecols=false)
df_VLAIO.PRODUCT_NC = lpad.(string.(df_VLAIO.PRODUCT_NC), 8, '0') # needs to be done again

# remove missing, only for VALUES
#   - drop missing [VALUE, QUANTITY]: lose 1-16027346/21667553 ~26% observations
#   - drop missing VALUE: lose 1-21666543/21667553 ~ 0.1% observations
subset!(df_VLAIO, :VALUE_IN_EUROS .=> ByRow(x -> !ismissing(x)))

# not needed in case of VALUES!!
# MAD adjustment
#   - cutoff = 3: lose further 1-15177297/21666543 ~30% observations (~30% of total)
#   - cutoff = 4: lose further 1-15812254/21666543 ~27% observations (~27% of total)
# outlier_cutoff = 4
# gdf = groupby(df_VLAIO, "PRODUCT_NC")
# df_VLAIO = transform(gdf, :VALUE_IN_EUROS => MAD_method => :MAD)
# subset!(df_VLAIO, :MAD => ByRow(x -> x < outlier_cutoff))
# df_VLAIO = df_VLAIO[:, Not(:MAD)]

# add EU/WORLD again
EU27 = ["Roemenië", "Griekenland", "Oostenrijk", "Polen", "Duitsland", "Spanje", "Hongarije", "Slovakije", "Italië", "Nederland",
       "Frankrijk", "Letland", "Kroatië", "Cyprus", "Malta", "Litouwen", "Slovenië", "Estland", "Portugal", "Finland", "Tsjechië",
       "Luxemburg", "Zweden", "Denemarken", "Bulgarije", "Ierland", "België"]

# slightly modified functions
function append_WORLD(df::DataFrame)

    # sum over PARTNER
    cols_grouping = ["DECLARANT_ISO", "PRODUCT_NC", "PRODUCT_NC_LAB1", "PRODUCT_NC_LAB2", "TRADE_TYPE", "FLOW", "PERIOD"]
    gdf = groupby(df, cols_grouping)
    df_WORLD = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)

    df_WORLD.TRADE_TYPE .= "total"
    df_WORLD.PARTNER_ISO .= "WORLD"

    df = vcat(df, df_WORLD)

    return df
end

function append_EU(df::DataFrame, EU::Vector{String})

    # subset and sum over EU ctrys
    df_EU = subset(df, :PARTNER_ISO => ByRow(x -> x in EU))
    cols_grouping = ["DECLARANT_ISO", "PRODUCT_NC", "PRODUCT_NC_LAB1", "PRODUCT_NC_LAB2", "TRADE_TYPE", "FLOW", "PERIOD"]
    gdf = groupby(df_EU, cols_grouping)
    df_EU = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)

    df_EU.TRADE_TYPE .= "intra"
    df_EU.PARTNER_ISO .= "EU"

    df = vcat(df, df_EU)

    return df
end

df_VLAIO = append_WORLD(df_VLAIO)
df_VLAIO = append_EU(df_VLAIO, EU27)

# double check if no zeros/missing introduced by WORLD/EU
transform!(df_VLAIO, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> ifelse(iszero(x) | ismissing(x), missing, x)), renamecols=false)
subset!(df_VLAIO, :VALUE_IN_EUROS => ByRow(x -> !ismissing(x))) # lose 0% observations

cols_subset = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD", "PRODUCT_NC", "VALUE_IN_EUROS", "QUANTITY_IN_KG"]
df_VLAIO = df_VLAIO[:, cols_subset]

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Figure 1
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# prepare data
df = copy(df_VLAIO)
partners = ["Verenigd Koninkrijk", "EU", "Duitsland", "Nederland", "Frankrijk", "Italië"]
subset!(df, :DECLARANT_ISO => ByRow(x -> x == "Vlaanderen"), :PARTNER_ISO => ByRow(x -> x in partners))
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)

# formatting
transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)

# need to sort before computing HP and 3MMA
sort!(df)

# HP filter
λ = 20
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, :VALUE_IN_EUROS => (x -> HP(x, λ)) => :VALUE_IN_EUROS_HP)

# create 3MMA as observations are very volatile
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, :VALUE_IN_EUROS => (x -> movingaverage(x,3)) => :VALUE_3MMA)

# rolling standard deviation
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, [:VALUE_IN_EUROS, :VALUE_3MMA] .=> (x -> rolling_std(x, 6)) .=> [:STD_VALUE, :STD_VALUE_3MMA])

# plotting
#sort!(df, :DATE)

for flow in ["imports", "exports"]

    # values
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x != "EU")) plot(:DATE, :VALUE_IN_EUROS/1e9,
    group=:PARTNER_ISO, lw=2, legend=:topleft, ylabel="euros (billion)", title="Flemish "* flow*": values")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/" * "fig1_" * flow * "_values" * ".png") # export image dropbox

    # HP
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x != "EU")) plot(:DATE, :VALUE_IN_EUROS_HP/1e9,
            group=:PARTNER_ISO, lw=2, legend=:topleft, ylabel="euros (billion)", title="Flemish "*flow*": values (HP, λ=$λ)")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/" * "fig1_" * flow * "_values" * "_HP" * ".png") # export image dropbox

    # 3MMA
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x != "EU")) plot(:DATE, :VALUE_3MMA/1e9,
            group=:PARTNER_ISO, lw=2, legend=:topleft, ylabel="euros (billion)", title="Flemish "*flow*": values (3MMA)")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/" * "fig1_" * flow * "_values" * "_3MMA" * ".png") # export image dropbox

    # STD
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x != "EU")) plot(:DATE, :STD_VALUE/1e9,
            group=:PARTNER_ISO, lw=2, legend=:topleft, ylabel="euros (billion)", title="Flemish "*flow*": 6 months STD \n (values, rolling window)")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/" * "fig4_" * flow * "_values" * "_STD" * ".png") # export image dropbox

end

# STD
p = @df subset(df, :PARTNER_ISO => ByRow(x -> x == "Verenigd Koninkrijk")) plot(:DATE, :STD_VALUE/1e9,
        group=:FLOW, lw=2, legend=:bottomleft, ylabel="euros (billion)", title="Flemish: 6 months STD \n (values, rolling window)")
vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/" * "fig4_" * "_values" * "_STD" * ".png") # export image dropbox

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Figure 2
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# prepare data
df = copy(df_VLAIO)
partners = ["Verenigd Koninkrijk"]
subset!(df, :DECLARANT_ISO => ByRow(x -> x == "Vlaanderen"), :PARTNER_ISO => ByRow(x -> x in partners))

# CPA classification
df_CN_PROD = DataFrame(XLSX.readtable(dir_dropbox* "results/correspondence/" * "table_CN_CPA21" * ".xlsx", "Sheet1")...)
transform!(df_CN_PROD, names(df_CN_PROD) .=> ByRow(string), renamecols=false)
#rename!(df_CN_PROD, :CN_code => :PRODUCT_NC)

df_join = leftjoin(df, df_CN_PROD, on=:PRODUCT_NC)
subset!(df_join, :CPA21_lab => ByRow(x -> !ismissing(x))) # lose 1-775329/801316 ~ 3%

cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD", "CPA21_2digits", "CPA21_lab", "CPA21_lab_custom"]
gdf = groupby(df_join, cols_grouping)
df = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)

# add DATE
transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)
sort!(df)

# HP filter
λ = 20
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "CPA21_2digits", "CPA21_lab", "CPA21_lab_custom"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, x -> nrow(x) < 3 ? DataFrame() : x) # remove groups with only two observations otherwise we cannot apply HP filter
gdf = groupby(df, cols_grouping)
df = transform(gdf, :VALUE_IN_EUROS => (x -> HP(x, λ)) => :VALUE_IN_EUROS_HP)

# plotting
sort!(df)
products = sort(unique(df.CPA21_2digits))
transform!(df, [:CPA21_2digits, :CPA21_lab_custom] => ByRow((x,y) -> x * " " * y) => :legend)

# -----------------------
# correspondance table with country names
# df_ind = DataFrame(XLSX.readtable(dir_dropbox * "results/correspondence/" * "table_NACE_custom" * ".xlsx", "Sheet1")...)
# transform!(df_ind, names(df_ind) .=> ByRow(string), renamecols=false)
# rename!(df_ind, :NACE_code => :PRODUCT_NC_digits)
# df = leftjoin(df, df_ind, on=:PRODUCT_NC_digits)

# df_ind = CSV.read(dir_dropbox * "rawdata/correspondence/CN/" * "CN_2022" * ".csv", DataFrame)
# transform!(df_ind, names(df_ind) .=> ByRow(string), renamecols=false)
# df_new = subset(df_ind, :Code_1 => ByRow(x -> length(x) == 2))
# subset!(df_new, :Parent_1 => ByRow(x -> x != "missing"))
# df_ind = df_new[:, ["Level", "Code", "Code_1", "Parent_1", "Description", "Self-explanatory texts in English"]]
# rename!(df_ind, names(df_ind) .=> ["level", "code_12dig", "code_2dig", "code_1dig", "lab", "lab_custom"])
# transform!(df_ind, :code_12dig => ByRow(x -> lpad(x, 12, '0')), renamecols=false)
# #XLSX.writetable(dir_dropbox * "results/correspondence/" * "table_HS2_correspondence" * ".xlsx", df_ind, overwrite=true)

# rename!(df_ind, :code_2dig => :PRODUCT_NC_digits)
# df = leftjoin(df, df_ind, on=:PRODUCT_NC_digits)
# transform!(df, :lab_custom => ByRow(string), renamecols=false)
# transform!(df, :lab_custom => ByRow(x -> first(x, 25)) => :lab_custom2)
# transform!(df, [:PRODUCT_NC_digits, :lab_custom2] => ByRow((x,y) -> x * " " * y) => :lab_custom3)

# -----------------------

for flow in ["imports", "exports"]
    for i in 1:3:length(products)

    product_range = products[i:i+2]

    # values
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :CPA21_2digits => ByRow(x -> x in product_range)) plot(:DATE, :VALUE_IN_EUROS/1e6,
            group=:legend, lw=2, legend=:bottomleft, ylabel="euros (million)", title="Flemish "*flow*": values")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/products/values/" * flow * "/" * "fig2_" * flow * "_" * product_range[1] * "_" * product_range[end] * "_values" * ".png") # export image

    # HP
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :CPA21_2digits => ByRow(x -> x in product_range)) plot(:DATE, :VALUE_IN_EUROS_HP/1e6,
            group=:legend, lw=2, legend=:bottomleft, ylabel="euros (million)", title="Flemish "*flow*": values (HP, λ=$λ)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/products/values/HP/" * flow * "/" * "fig2_" * flow * "_" * product_range[1] * "_" * product_range[end] * "_values" * "_HP" * ".png") # export image

    end
end

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Figure 3
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# prepare data
df = copy(df_VLAIO)
partners = ["Verenigd Koninkrijk", "EU", "Duitsland", "Nederland", "Frankrijk", "WORLD", "Italië"]
subset!(df, :DECLARANT_ISO => ByRow(x -> x == "Vlaanderen"), :PARTNER_ISO => ByRow(x -> x in partners))
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)

df_WORLD = subset(df, :PARTNER_ISO => ByRow(x -> x == "WORLD"))
cols_grouping = ["DECLARANT_ISO", "FLOW", "PERIOD"]
gdf = groupby(df_WORLD, cols_grouping)
df_WORLD = combine(gdf, :VALUE_IN_EUROS => sum => :VALUE_WORLD)

df = leftjoin(subset(df, :PARTNER_ISO => ByRow(x -> x != "WORLD")), df_WORLD, on=[:DECLARANT_ISO, :FLOW, :PERIOD])

transform!(df, [:VALUE_IN_EUROS, :VALUE_WORLD] => ByRow((x,s) -> x/s) => :SHARE_VALUE)

# Date
transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)
sort!(df)

# HP filter
λ = 20
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, :SHARE_VALUE => (x -> HP(x, λ)) => :SHARE_VALUE_HP)

# plotting
partners = ["Verenigd Koninkrijk", "Duitsland", "Nederland", "Frankrijk", "Italië"]

for flow in ["imports", "exports"]

    # values
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x in partners)) plot(:DATE, :SHARE_VALUE*100,
            group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="percentages", title="Flemish "* flow*" share (values)")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/" * "fig3_" * flow * "_value_share" * ".png") # export image dropbox

    # HP
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x in partners)) plot(:DATE, :SHARE_VALUE_HP*100,
        group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="percentages", title="Flemish "* flow*" share (values, λ=$λ)")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/" * "fig3_" * flow * "_value_share" * "_HP" * ".png") # export image dropbox

end

# VLA vis-a-vis GB
p = @df subset(df, :PARTNER_ISO => ByRow(x -> x == "Verenigd Koninkrijk")) plot(:DATE, :SHARE_VALUE_HP*100,
        group=:FLOW, lw=2, legend=:bottomleft, ylabel="percentages", title="Total trade share (values, λ=$λ): \n Flanders vis-a-vis Great Britain")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/" * "fig3_" * "VLA_GB" * "_value_share" * "_HP" * ".png") # export image

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Figure 5
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# prepare data
df = copy(df_VLAIO)
partners = ["Verenigd Koninkrijk", "Duitsland", "Nederland", "Frankrijk", "Italië"]
subset!(df, :DECLARANT_ISO => ByRow(x -> x == "Vlaanderen"), :PARTNER_ISO => ByRow(x -> x in partners))

cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, nrow => :COUNT)

transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)
sort!(df)

# HP filter
λ = 20
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, :COUNT => (x -> HP(x, λ)) => :COUNT_HP)

for flow in ["imports", "exports"]

    p = @df subset(df, :FLOW => ByRow(x -> x == flow)) plot(:DATE, :COUNT,
            group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="number of products", title="Flemish "* flow*": number of products")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/" * "fig5_" * flow * "_product_count" * ".png") # export image dropbox

    p = @df subset(df, :FLOW => ByRow(x -> x == flow)) plot(:DATE, :COUNT_HP,
    group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="number of products", title="Flemish "* flow*": number of products (HP, λ=$λ")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/" * "fig5_" * flow * "_product_count_HP" * ".png") # export image dropbox

end

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Table 1
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

df = copy(df_VLAIO)
partners = ["Verenigd Koninkrijk", "EU", "Duitsland", "Nederland", "Frankrijk"]
subset!(df, :DECLARANT_ISO => ByRow(x -> x == "Vlaanderen"), :PARTNER_ISO => ByRow(x -> x in [partners; "WORLD"]))
transform!(df, :PRODUCT_NC => ByRow(x -> x[1:2]) => :PRODUCT_NC_digits)
transform!(df, :PERIOD => ByRow(x -> string(x)[1:4]) => :YEAR)

cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "YEAR", "PRODUCT_NC_digits"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)

df_join = leftjoin(subset(df, :PARTNER_ISO => ByRow(x -> x != "WORLD")),
                    subset(df, :PARTNER_ISO => ByRow(x -> x == "WORLD")), on=[:DECLARANT_ISO, :FLOW, :PRODUCT_NC_digits, :YEAR], makeunique=true)

rename!(df_join, [:VALUE_IN_EUROS, :VALUE_IN_EUROS_1] .=> [:VALUE_PARTNER, :VALUE_WORLD])
transform!(df_join, [:VALUE_PARTNER, :VALUE_WORLD] => ByRow((x,s) -> x/s) => :SHARE_VALUE)

cols_name = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PRODUCT_NC_digits", "YEAR", "VALUE_PARTNER", "VALUE_WORLD", "SHARE_VALUE"]
df = df_join[:, cols_name]

# ----------------

df_tab1 = subset(df, :PARTNER_ISO => ByRow(x -> x == "Verenigd Koninkrijk"))
cols_name = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PRODUCT_NC_digits", "YEAR", "SHARE_VALUE"]
df_tab1 = df_tab1[:, cols_name]

sort!(df_tab1)
transform!(df_tab1, :SHARE_VALUE => ByRow(x -> round(x*100, digits=2)), renamecols=false)
df_tab1_wide = unstack(df_tab1, :YEAR, :SHARE_VALUE)

# XLSX.writetable(dir_dropbox * "results/images/VLAIO/summary_stats/" * "table1_within_prod_importance" * ".xlsx", df_tab1_wide, overwrite=true)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Table 2
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

df = copy(df_VLAIO)
partners = ["Verenigd Koninkrijk", "EU", "Duitsland", "Nederland", "Frankrijk"]
subset!(df, :DECLARANT_ISO => ByRow(x -> x == "Vlaanderen"), :PARTNER_ISO => ByRow(x -> x in [partners; "WORLD"]))
transform!(df, :PRODUCT_NC => ByRow(x -> x[1:2]) => :PRODUCT_NC_digits)
transform!(df, :PERIOD => ByRow(x -> string(x)[1:4]) => :YEAR)

cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "YEAR", "PRODUCT_NC_digits"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, :VALUE_IN_EUROS => sum => :VALUE_PARTNER)

df_WORLD = subset(df, :PARTNER_ISO => ByRow(x -> x == "WORLD"))
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "YEAR"]
gdf = groupby(df_WORLD, cols_grouping)
df_WORLD = combine(gdf, :VALUE_PARTNER => sum => :VALUE_WORLD)

df_join = leftjoin(subset(df, :PARTNER_ISO => ByRow(x -> x != "WORLD")), df_WORLD[:,Not(:PARTNER_ISO)], on=[:DECLARANT_ISO, :FLOW, :YEAR])
transform!(df_join, [:VALUE_PARTNER, :VALUE_WORLD] => ByRow((x,s) -> x/s) => :SHARE_VALUE)

df_tab2 = subset(df_join, :PARTNER_ISO => ByRow(x -> x == "Verenigd Koninkrijk"))
cols_name = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PRODUCT_NC_digits", "YEAR", "SHARE_VALUE"]
df_tab2 = df_tab2[:, cols_name]

sort!(df_tab2)
transform!(df_tab2, :SHARE_VALUE => ByRow(x -> round(x*100, digits=2)), renamecols=false)

df_tab2_wide = unstack(df_tab2, :YEAR, :SHARE_VALUE)

# XLSX.writetable(dir_dropbox * "results/images/VLAIO/summary_stats/" * "table2_overall_prod_importance" * ".xlsx", df_tab2_wide, overwrite=true)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Figure 6
# y-on-y monthly growth rates
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

df = copy(df_VLAIO)
partners = ["Verenigd Koninkrijk", "EU", "Duitsland", "Nederland", "Frankrijk", "Italië"]
subset!(df, :DECLARANT_ISO => ByRow(x -> x == "Vlaanderen"), :PARTNER_ISO => ByRow(x -> x in [partners; "WORLD"]))
transform!(df, :PRODUCT_NC => ByRow(x -> x[1:2]) => :PRODUCT_NC_digits)

# aggregate to two digits
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD", "PRODUCT_NC_digits"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)

# YoY percentage change
#   - a simple percentage change X(t)-X(t-1)/X(t-1)
#   - computes percentage change if difference between PERIOD == 100 (201512 - 201412 = 100)
function yoy_change(period::AbstractVector, input::AbstractVector)
    M = [period[i]-period[j] == 100 ? (input[i]-input[j])/input[j] : missing for i in eachindex(input), j in eachindex(input)] # matrix
    V = [all(ismissing.(M[i,:])) ? missing : M[i, findfirst(typeof.(M[i,:]) .== Float64)] for i in 1:size(M, 1)] # reduce to vector
    return V
end

# compute YOY monthly percentage change
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "PRODUCT_NC_digits", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, [:PERIOD, :VALUE_IN_EUROS] => yoy_change => :YOY_VALUE)

# ------------

# MAD adjustment
#   - need to drop missing
outlier_cutoff = 4
gdf = groupby(subset(df, :YOY_VALUE => ByRow(x -> !ismissing(x))), "PRODUCT_NC_digits")
df = transform(gdf, :YOY_VALUE => MAD_method => :MAD)
subset!(df, :MAD => ByRow(x -> x < outlier_cutoff))
df = df[:, Not(:MAD)]

# ------------

# Date
transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)
sort!(df)

# HP filter
λ = 20
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PRODUCT_NC_digits"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, :YOY_VALUE => (x -> HP(x, λ)) => :YOY_VALUE_HP)

# plotting
products = 1:8:89

for flow in ["imports", "exports"]
    for i in products

    product_range = lpad.(i:i+6,2,'0')

    # values
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PRODUCT_NC_digits => ByRow(x -> x in product_range), :PARTNER_ISO => ByRow(x -> x == "Verenigd Koninkrijk")) plot(:DATE, :YOY_VALUE*100,
            group=:PRODUCT_NC_digits, lw=2, legend=:bottomleft, ylabel="percentages", title="Flemish "*flow*": monthly y-on-y change")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/products/yoy_change/" * flow * "/" * "fig6_" * flow * "_" * product_range[1] * "_" * product_range[end] * "_yoy_change" * ".png") # export image

    # HP
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PRODUCT_NC_digits => ByRow(x -> x in product_range), :PARTNER_ISO => ByRow(x -> x == "Verenigd Koninkrijk")) plot(:DATE, :YOY_VALUE_HP*100,
            group=:PRODUCT_NC_digits, lw=2, legend=:bottomleft, ylabel="percentages", title="Flemish "*flow*": monthly y-on-y change \n (values, λ=$λ)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/products/yoy_change/HP/" * flow * "/" * "fig2_" * flow * "_" * product_range[1] * "_" * product_range[end] * "_yoy_change" * "_HP" * ".png") # export image

    end
end

df = copy(df_VLAIO)
partners = ["Verenigd Koninkrijk", "EU", "Duitsland", "Nederland", "Frankrijk", "Italië"]
subset!(df, :DECLARANT_ISO => ByRow(x -> x == "Vlaanderen"), :PARTNER_ISO => ByRow(x -> x in [partners; "WORLD"]))

# CPA classification
df_CN_PROD = DataFrame(XLSX.readtable(dir_dropbox* "results/correspondence/" * "table_CN_CPA21" * ".xlsx", "Sheet1")...)
transform!(df_CN_PROD, names(df_CN_PROD) .=> ByRow(string), renamecols=false)
#rename!(df_CN_PROD, :CN_code => :PRODUCT_NC)

df_join = leftjoin(df, df_CN_PROD, on=:PRODUCT_NC)
subset!(df_join, :CPA21_lab => ByRow(x -> !ismissing(x))) # lose 1-775329/801316 ~ 3%

cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD", "CPA21_2digits", "CPA21_lab", "CPA21_lab_custom"]
gdf = groupby(df_join, cols_grouping)
df = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)

# compute YOY monthly percentage change
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "CPA21_2digits", "CPA21_lab", "CPA21_lab_custom"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, [:PERIOD, :VALUE_IN_EUROS] => yoy_change => :YOY_VALUE)

# ------------

# MAD adjustment
#   - need to drop missing
outlier_cutoff = 4
gdf = groupby(subset(df, :YOY_VALUE => ByRow(x -> !ismissing(x))), ["CPA21_2digits", "CPA21_lab", "CPA21_lab_custom"])
df = transform(gdf, :YOY_VALUE => MAD_method => :MAD)
subset!(df, :MAD => ByRow(x -> x < outlier_cutoff))
df = df[:, Not(:MAD)]

# ------------

# Date
transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)
sort!(df)

# HP filter
λ = 20
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "CPA21_2digits", "CPA21_lab", "CPA21_lab_custom"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, x -> nrow(x) < 3 ? DataFrame() : x) # remove groups with only two observations otherwise we cannot apply HP filter
gdf = groupby(df, cols_grouping)
df = transform(gdf, :YOY_VALUE => (x -> HP(x, λ)) => :YOY_VALUE_HP)

# plotting
sort!(df)
products = sort(unique(df.CPA21_2digits))
transform!(df, [:CPA21_2digits, :CPA21_lab_custom] => ByRow((x,y) -> x * " " * y) => :legend)

for flow in ["imports", "exports"]
    for i in 1:4:length(products)

    product_range = products[i:i+3]

    # values
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :CPA21_2digits => ByRow(x -> x in product_range), :PARTNER_ISO => ByRow(x -> x == "Verenigd Koninkrijk")) plot(:DATE, :YOY_VALUE*100,
            group=:legend, lw=2, legend=:bottomleft, ylabel="percentages", title="Flemish "*flow*": monthly y-on-y change")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/products/yoy_change/" * flow * "/" * "fig6_" * flow * "_" * product_range[1] * "_" * product_range[end] * "_yoy_change" * ".png") # export image

    # HP
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :CPA21_2digits => ByRow(x -> x in product_range), :PARTNER_ISO => ByRow(x -> x == "Verenigd Koninkrijk")) plot(:DATE, :YOY_VALUE_HP*100,
            group=:legend, lw=2, legend=:bottomleft, ylabel="percentages", title="Flemish "*flow*": monthly y-on-y change \n (values, λ=$λ)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/products/yoy_change/HP/" * flow * "/" * "fig2_" * flow * "_" * product_range[1] * "_" * product_range[end] * "_yoy_change" * "_HP" * ".png") # export image

    end
end

# ------------
# standard deviation of y-on-y growth
# ------------

# table
referendum = Date(2016, 07, 01)
brexit = Date(2020, 02, 01)
trade = Date(2021, 05, 01)

transform!(df, :DATE => ByRow(x -> ifelse(x in referendum-Month(1)-Year(1):Month(1):referendum-Month(1), "pre",
    ifelse(x in referendum:Month(1):brexit-Month(1), "brexit_1", ifelse(x in brexit:Month(1):trade-Month(1), "brexit_2",
        ifelse(x in trade:Month(1):trade+Month(10), "brexit_3", "no"))))) => :BREXIT)

cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "BREXIT"]
gdf = groupby(df, cols_grouping)
tab_std = combine(gdf, :YOY_VALUE => std => :STD_VALUE)

tab_std_wide = unstack(tab_std, :BREXIT, :STD_VALUE)
sort!(tab_std_wide, :FLOW)

# XLSX.writetable(dir_dropbox * "results/images/VLAIO/summary_stats/" * "table3_std" * ".xlsx", tab_std_wide, overwrite=true)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Figure 7
# STD of y-on-y monthly growth rates
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

df = copy(df_VLAIO)
partners = ["Verenigd Koninkrijk", "Duitsland", "Nederland", "Frankrijk", "Italië"]
subset!(df, :DECLARANT_ISO => ByRow(x -> x == "Vlaanderen"), :PARTNER_ISO => ByRow(x -> x in partners))

# aggregate over products
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)

# compute YOY monthly percentage change
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, [:PERIOD, :VALUE_IN_EUROS] => yoy_change => :YOY_VALUE)

# # ------------

# # MAD adjustment
# #   - need to drop missing
# outlier_cutoff = 4
# gdf = groupby(subset(df, :YOY_VALUE => ByRow(x -> !ismissing(x))), "PRODUCT_NC_digits")
# df = transform(gdf, :YOY_VALUE => MAD_method => :MAD)
# subset!(df, :MAD => ByRow(x -> x < outlier_cutoff))
# df = df[:, Not(:MAD)]

# # ------------

# Date
transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)
sort!(df)

# HP filter
subset!(df, :YOY_VALUE => ByRow(x -> !ismissing(x))) # need to drop missing
λ = 20
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, :YOY_VALUE => (x -> HP(x, λ)) => :YOY_VALUE_HP)

# std
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, [:YOY_VALUE, :YOY_VALUE_HP] .=> (x -> rolling_std(x, 6)) .=> [:STD_YOY, :STD_YOY_HP])

for flow in ["imports", "exports"]

    # YOY HP
    p = @df subset(df, :FLOW => ByRow(x -> x == flow)) plot(:DATE, :STD_YOY_HP,
            group=:PARTNER_ISO, lw=2, legend=:topleft, ylabel="percentages", title="Flemish "*flow*" : 6 months STD \n (YOY change, rolling window, λ=$λ)")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/" * "fig7_" * flow * "_YOY" * "_STD" * ".png") # export image dropbox

end

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Joeps figures from 14/03/2022
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# --------------
# Products per year
# --------------

# prepare data
df = copy(df_VLAIO)
partners = ["Verenigd Koninkrijk", "Duitsland", "Nederland", "Frankrijk", "Italië", "WORLD"]
subset!(df, :DECLARANT_ISO => ByRow(x -> x == "Vlaanderen"), :PARTNER_ISO => ByRow(x -> x in partners))
transform!(df, :PERIOD => ByRow(x -> Date(string(x)[1:4], DateFormat("yyyy"))) => :DATE)

cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "DATE", "PRODUCT_NC"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)

cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "DATE"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, nrow => :COUNT)

sort!(df)
transform!(df, :DATE => ByRow(x -> Dates.year(x)) => :YEAR)

for flow in ["imports", "exports"]

    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x in ["Verenigd Koninkrijk", "WORLD"])) groupedbar(:YEAR, :COUNT,
            group=:PARTNER_ISO, bar_position=:dodge, legend=:bottomleft, ylabel="number of products", title="Flemish "* flow*": number of products")
    #vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    #vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    #vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/" * "fig5_" * flow * "_product_count_year" * ".png") # export image dropbox

end

# absolute difference between successive periods
function diff(period::AbstractVector, input::AbstractVector)
    [i == 1 ? missing : period[i]-period[i-1] != 1 ? missing : input[i]-input[i-1] for i in eachindex(input)]
end

cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, [:YEAR, :COUNT] => diff => :DIFF)
transform!(df, [:COUNT, :DIFF] => ByRow((x,y) -> x - y) => :REMAIN)

for flow in ["imports", "exports"]

    p = @df subset(df, :DIFF => ByRow(x->!ismissing(x)), :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x in ["Verenigd Koninkrijk", "WORLD"])) groupedbar(:YEAR, :DIFF,
            group=:PARTNER_ISO, bar_position=:dodge, legend=:topleft, ylabel="number of products", title="Flemish "* flow*": \n y-on-y change in number of products")
    #vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    #vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    #vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/" * "fig5_" * flow * "_product_diff_year" * ".png") # export image dropbox

end

# --------------
# averages
# --------------

referendum = Date(2016, 07, 01)
trade = Date(2021, 05, 01)
pre_t = referendum-Month(6+1):Month(1):referendum-Month(1)
now_t = trade:Month(1):trade+Month(6)

# prepare data
df = copy(df_VLAIO)
partners = ["Verenigd Koninkrijk", "Duitsland", "Nederland", "Frankrijk"]
subset!(df, :DECLARANT_ISO => ByRow(x -> x == "Vlaanderen"), :PARTNER_ISO => ByRow(x -> x in partners))
transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)
transform!(df, :DATE => ByRow(x -> ifelse(x in pre_t, "pre", ifelse(x in now_t, "now", "delete"))) => :dummy_time)

# average flow per product

# aggregate over time
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PRODUCT_NC", "dummy_time"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)

# averages over products
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "dummy_time"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, nrow => :COUNT, :VALUE_IN_EUROS .=> [sum mean median std])
rename!(df, [:6, :7, :8, :9].=> [:sum, :mean, :median, :std])
subset!(df, :dummy_time => ByRow(x -> x in ["now", "pre"]))

# million euros
transform!(df, [:sum, :mean, :median, :std] .=> ByRow(x -> x/1e6), renamecols=false)
sort!(df, :FLOW)

XLSX.writetable(dir_dropbox * "results/images/VLAIO/summary_stats/" * "table4_averages" * ".xlsx", df, overwrite=true)

# --------------
# follow top 10
# --------------

# # 2 digit codes of top 10 from 2014
# prod_exports = ["87", "30", "27", "39", "84", "29", "85", "71", "90", "64"]
# prod_imports = ["27", "87", "84", "29", "71", "39", "85", "40", "90", "30"]

# 2 digit codes of top 10 from 2021
prod_exports = ["30", "27", "87", "39", "71", "84", "22", "29", "90", "20"]
prod_imports = ["27", "87", "84", "29", "71", "39", "85", "40", "90", "30"]

df = copy(df_VLAIO)
partners = ["Verenigd Koninkrijk", "EU", "Duitsland", "Nederland", "Frankrijk"]
subset!(df, :DECLARANT_ISO => ByRow(x -> x == "Vlaanderen"), :PARTNER_ISO => ByRow(x -> x in [partners; "WORLD"]))
transform!(df, :PRODUCT_NC => ByRow(x -> x[1:2]) => :PRODUCT_NC_digits)
transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)
transform!(df, :PERIOD => ByRow(x -> Date(string(x)[1:4], DateFormat("yyyy"))) => :YEAR)

cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD", "DATE", "PRODUCT_NC_digits"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)

cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD", "DATE"]
gdf = groupby(df, cols_grouping)
df_total_partner = combine(gdf, :VALUE_IN_EUROS => sum => :VALUE_TOTAL)

# subset products and produce share
subset!(df, :PRODUCT_NC_digits => ByRow(x -> x in [prod_imports; prod_exports]))
df = leftjoin(df, df_total_partner, on=[:DECLARANT_ISO, :PARTNER_ISO, :FLOW, :PERIOD, :DATE])
transform!(df, [:VALUE_IN_EUROS, :VALUE_TOTAL] => ByRow((x, y) -> x/y) => :SHARE_VALUE)

sort!(df)

df_ind = CSV.read(dir_dropbox * "rawdata/correspondence/CN/" * "CN_2022" * ".csv", DataFrame)
transform!(df_ind, names(df_ind) .=> ByRow(string), renamecols=false)
df_new = subset(df_ind, :Code_1 => ByRow(x -> length(x) == 2))
subset!(df_new, :Parent_1 => ByRow(x -> x != "missing"))
df_ind = df_new[:, ["Level", "Code", "Code_1", "Parent_1", "Description", "Self-explanatory texts in English"]]
rename!(df_ind, names(df_ind) .=> ["level", "code_12dig", "code_2dig", "code_1dig", "lab", "lab_custom"])
transform!(df_ind, :code_12dig => ByRow(x -> lpad(x, 12, '0')), renamecols=false)
#XLSX.writetable(dir_dropbox * "results/correspondence/" * "table_HS2_correspondence" * ".xlsx", df_ind, overwrite=true)

rename!(df_ind, :code_2dig => :PRODUCT_NC_digits)
df = leftjoin(df, df_ind, on=:PRODUCT_NC_digits)
transform!(df, :lab_custom => ByRow(string), renamecols=false)
transform!(df, :lab_custom => ByRow(x -> first(x, 25)) => :lab_custom2)
transform!(df, [:PRODUCT_NC_digits, :lab_custom2] => ByRow((x,y) -> x * " " * y) => :lab_custom3)

for flow in ["imports", "exports"]

    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x in ["Verenigd Koninkrijk"])) groupedbar(:DATE, :SHARE_VALUE*100,
            group=:lab_custom3, bar_position=:stack, legend=:bottomleft, ylabel="percentages", title="Flemish "* flow*" with GB: \n Top 10 industries")
    #vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    #vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    #vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/" * "fig8_" * flow * "_values_monthly" * "_top10" * ".png") # export image dropbox

end

# -----
# yearly

df = copy(df_VLAIO)
partners = ["Verenigd Koninkrijk", "EU", "Duitsland", "Nederland", "Frankrijk"]
subset!(df, :DECLARANT_ISO => ByRow(x -> x == "Vlaanderen"), :PARTNER_ISO => ByRow(x -> x in [partners; "WORLD"]))
transform!(df, :PRODUCT_NC => ByRow(x -> x[1:2]) => :PRODUCT_NC_digits)
transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)
transform!(df, :PERIOD => ByRow(x -> Date(string(x)[1:4], DateFormat("yyyy"))) => :YEAR)

cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "YEAR", "PRODUCT_NC_digits"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)

cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "YEAR", "FLOW"]
gdf = groupby(df, cols_grouping)
df_total_partner = combine(gdf, :VALUE_IN_EUROS => sum => :VALUE_TOTAL)

# subset products and produce share
subset!(df, :PRODUCT_NC_digits => ByRow(x -> x in [prod_imports; prod_exports]))
df = leftjoin(df, df_total_partner, on=[:DECLARANT_ISO, :PARTNER_ISO, :FLOW, :YEAR])
transform!(df, [:VALUE_IN_EUROS, :VALUE_TOTAL] => ByRow((x, y) -> x/y) => :SHARE_VALUE)

transform!(df, :YEAR => ByRow(x -> string(x)[1:4]), renamecols=false)

sort!(df)

df_ind = CSV.read(dir_dropbox * "rawdata/correspondence/CN/" * "CN_2022" * ".csv", DataFrame)
transform!(df_ind, names(df_ind) .=> ByRow(string), renamecols=false)
df_new = subset(df_ind, :Code_1 => ByRow(x -> length(x) == 2))
subset!(df_new, :Parent_1 => ByRow(x -> x != "missing"))
df_ind = df_new[:, ["Level", "Code", "Code_1", "Parent_1", "Description", "Self-explanatory texts in English"]]
rename!(df_ind, names(df_ind) .=> ["level", "code_12dig", "code_2dig", "code_1dig", "lab", "lab_custom"])
transform!(df_ind, :code_12dig => ByRow(x -> lpad(x, 12, '0')), renamecols=false)
#XLSX.writetable(dir_dropbox * "results/correspondence/" * "table_HS2_correspondence" * ".xlsx", df_ind, overwrite=true)

rename!(df_ind, :code_2dig => :PRODUCT_NC_digits)
df = leftjoin(df, df_ind, on=:PRODUCT_NC_digits)
transform!(df, :lab_custom => ByRow(string), renamecols=false)
transform!(df, :lab_custom => ByRow(x -> first(x, 25)) => :lab_custom2)
transform!(df, [:PRODUCT_NC_digits, :lab_custom2] => ByRow((x,y) -> x * " " * y) => :lab_custom3)

for flow in ["imports", "exports"]

    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x in ["Verenigd Koninkrijk"])) groupedbar(:YEAR, :SHARE_VALUE*100,
            group=:lab_custom3, bar_position=:stack, legend=:bottomleft, ylabel="percentages", title="Flemish "* flow*" with GB: \n Top 10 industries")
    #vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    #vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    #vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/" * "fig8_" * flow * "_values_annual" * "_top10_2021" * ".png") # export image dropbox

end

# --------------- CPA21

df = copy(df_VLAIO)
partners = ["Verenigd Koninkrijk", "EU", "Duitsland", "Nederland", "Frankrijk"]
subset!(df, :DECLARANT_ISO => ByRow(x -> x == "Vlaanderen"), :PARTNER_ISO => ByRow(x -> x in [partners; "WORLD"]))
transform!(df, :PERIOD => ByRow(x -> Date(string(x)[1:4], DateFormat("yyyy"))) => :YEAR)

# CPA classification
df_CN_PROD = DataFrame(XLSX.readtable(dir_dropbox* "results/correspondence/" * "table_CN_CPA21" * ".xlsx", "Sheet1")...)
transform!(df_CN_PROD, names(df_CN_PROD) .=> ByRow(string), renamecols=false)
#rename!(df_CN_PROD, :CN_code => :PRODUCT_NC)

df_join = leftjoin(df, df_CN_PROD, on=:PRODUCT_NC)
subset!(df_join, :CPA21_lab => ByRow(x -> !ismissing(x))) # lose 1-7019070/7260237 ~ 3%

# aggregate over products and country
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "YEAR", "CPA21_2digits" , "CPA21_lab", "CPA21_lab_custom"]
gdf = groupby(df_join, cols_grouping)
df = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)

# aggregate over products
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "YEAR", "FLOW"]
gdf = groupby(df_join, cols_grouping)
df_total_partner = combine(gdf, :VALUE_IN_EUROS => sum => :VALUE_TOTAL)
transform!(df_total_partner, :YEAR => ByRow(x -> string(x)[1:4]), renamecols=false)

# subset products top 10 and produce share
transform!(df, :YEAR => ByRow(x -> string(x)[1:4]), renamecols=false)
top10_exports = first(sort(subset(df, :PARTNER_ISO => ByRow(x -> x == "Verenigd Koninkrijk"), :YEAR => ByRow(x -> x == "2021"),
     :FLOW => ByRow(x -> x == "exports")), order(:VALUE_IN_EUROS, rev=true)).CPA21_2digits, 10)
top10_imports = first(sort(subset(df, :PARTNER_ISO => ByRow(x -> x == "Verenigd Koninkrijk"), :YEAR => ByRow(x -> x == "2021"),
     :FLOW => ByRow(x -> x == "imports")), order(:VALUE_IN_EUROS, rev=true)).CPA21_2digits, 10)

# join
df = leftjoin(df, df_total_partner, on=[:DECLARANT_ISO, :PARTNER_ISO, :FLOW, :YEAR])
transform!(df, [:VALUE_IN_EUROS, :VALUE_TOTAL] => ByRow((x, y) -> x/y) => :SHARE_VALUE)

sort!(df)
transform!(df, [:CPA21_2digits, :CPA21_lab_custom] => ByRow((x,y) -> x * " " * y) => :legend)

flow = "exports"

p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x in ["Verenigd Koninkrijk"]), :CPA21_2digits => ByRow(x -> x in top10_exports)) groupedbar(:YEAR, :SHARE_VALUE*100,
    group=:legend, bar_position=:stack, legend=:outertopleft, ylabel="percentages", title="Flemish "* flow*" with GB: \n Top 10 industries", xrotation=45)
#vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
#vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
#vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/" * "fig8_" * flow * "_values_annual" * "_top10_2021" * ".png") # export image dropbox

flow = "imports"

p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x in ["Verenigd Koninkrijk"]), :CPA21_2digits => ByRow(x -> x in top10_exports)) groupedbar(:YEAR, :SHARE_VALUE*100,
    group=:legend, bar_position=:stack, legend=:outertopleft, ylabel="percentages", title="Flemish "* flow*" with GB: \n Top 10 industries", xrotation=45)
#vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
#vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
#vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
savefig(p, dir_dropbox * "results/images/VLAIO/summary_stats/" * "fig8_" * flow * "_values_annual" * "_top10_2021" * ".png") # export image dropbox


text:
using DelimitedFiles

function solve(input)
    raw = vcat(readdlm(input,Int16)...)
    p1 = process(copy(raw),aggregate)
    @show p1
    p2 = process(copy(raw),findRootValue)
    @show p2

end

function process(list, f)
    numChildren = popfirst!(list)
    numMeta = popfirst!(list)
    return f(ntimes(numChildren,list,f), splice!(list,1:numMeta))
end

function ntimes(n,list, f)
    result = Int16[]
    for i in 1:n
        push!(result,process(list,f))
    end
    return result
end

function aggregate(child,metadata)
    return sum(child) + sum(metadata)
end

function childValue((child,value),index)
    if length(child)>index
        value+=child[index+1]
    else
        value+= 0
    end
    return (child,value)
end

function findRootValue(child,metadata)
    pushfirst!(child,0)
    if length(child) == 1
        return sum(metadata)
    else
        return reduce(childValue,metadata,init=(child,0))[2]
    end
end

solve("input.txt")


text:
module DiffieHellman

using Random
#using Serialization
using SecureIO

using CryptoGroups

const _default_rng = Ref{RandomDevice}()
function __init__()
    _default_rng[] = RandomDevice()
end

default_rng() = _default_rng[]

function rngint(rng::AbstractRNG, len::Integer)
    max_n = ( BigInt(1) << len ) - 1
    if len > 2
        min_n = BigInt(1) << (len - 1)
        return rand(rng, min_n:max_n)
    end
    return rand(rng, 1:max_n)
end

# Probably before user descides to contact the server he asks for the certificate of the public key and checks if that is valid.
function diffie(s,sign::Function,verify::Function,G::AbstractGroup,rng::AbstractRNG)
    #@show "Diffie", typeof(s), hasmethod(serialize,(typeof(s),Any))
    serialize(s,G)

    B,Bsign = deserialize(s)

    #@show verify(B,Bsign)

    if verify(B,Bsign)

        t = security(G)
        a = rngint(rng,t)
        A = binary(G^a)

        serialize(s,(A,sign(A)))

        Bb = typeof(G)(B,G)
        key = value(Bb^a)
        return key
    else
        return Error("Key exchange failed.")
    end
end

diffie(io,sign::Function,verify::Function,G::AbstractGroup) = diffie(io,sign,verify,G,default_rng())

"""
This one returns a secret connection between two fixed parties. The signature function sign returns signature and the group with respect to which the signature was signed.
"""
function hellman(s,sign::Function,verify::Function,rng::AbstractRNG)
    #@show "Hellman"
    G = deserialize(s)

    t = security(G)
    b = rngint(rng,t)

    B = binary(G^b)
    serialize(s,(B,sign(B)))

    A,Asign = deserialize(s)

    if verify(A,Asign)
        Aa = typeof(G)(A,G)
        key = value(Aa^b)
        return key
    else
        return Error("Key exchange failed.")
    end
end

hellman(io,sign::Function,verify::Function) = hellman(io,sign,verify,default_rng())

export diffie, hellman

end # module


text:
struct GhostDerivativeOperator{T, E<:AbstractDiffEqLinearOperator{T}, F<:AbstractBC{T}} <: AbstractDiffEqLinearOperator{T}
    L :: E
    Q :: F
end

function *(L::AbstractDiffEqLinearOperator{T}, Q::AbstractBC{T}) where{T}
    return GhostDerivativeOperator{T, typeof(L), typeof(Q)}(L,Q)
end

function *(L::AbstractDiffEqCompositeOperator{T}, Q::AbstractBC{T}) where{T}
    return sum(map(op -> op * Q, L.ops))
end

function *(A::GhostDerivativeOperator{T1}, u::AbstractArray{T2}) where {T1,T2}
    #TODO Implement a function domaincheck(L::AbstractDiffEqLinearOperator, u) to see if components of L along each dimension match the size of u
    x = zeros(promote_type(T1,T2), unpadded_size(u))
    LinearAlgebra.mul!(x, A.L, A.Q*u)
    return x
end

function \(A::GhostDerivativeOperator, u::AbstractArray) # FIXME should have T1,T2 and promote result
    #TODO implement check that A has compatible size with u
    s = size(u)
    (A_l,A_b) = sparse(A, s)
    x = A_l \ Vector(reshape(u, length(u)) .- A_b) #Has to be converted to vector to work, A_b being sparse was causing a conversion to sparse.
    return reshape(x, s)
end

function \(A::GhostDerivativeOperator, u::AbstractVector) # FIXME as above, should promote_type(T1,T2)
    # TODO: is this specialization to u::AbstractVector really any faster?
    @assert length(u) == size(A.L, 1)
    (A_l,A_b) = sparse(A,length(u))
    A_l \ Vector(u .- A_b)
end

# update coefficients
function DiffEqBase.update_coefficients!(A::GhostDerivativeOperator,u,p,t)
    DiffEqBase.update_coefficients!(A.L,u,p,t)
end

# Implement multiplication for coefficients
function *(c::Number, A::GhostDerivativeOperator)
    (c * A.L) * A.Q
end

function *(c::Vector{<:Number}, A::GhostDerivativeOperator)
    (c * A.L) * A.Q
end

function *(coeff_func::Function, A::GhostDerivativeOperator)
    (coeff_func*A.L)*A.Q
end

# length and sizes
Base.ndims(A::GhostDerivativeOperator) = 2
Base.size(A::GhostDerivativeOperator) = (A.L.len, A.L.len)
Base.size(A::GhostDerivativeOperator,i::Integer) = size(A)[i]
Base.length(A::GhostDerivativeOperator) = reduce(*, size(A))

@inline ==(A1::GhostDerivativeOperator, A2::GhostDerivativeOperator) = A1.L == A2.L && A1.Q == A2.Q


text:
type Trade
    tradeID::Int
    kind::Int
    Pes::Float64
    Peb::Float64
    As::Int #Agent - Seller
    Ab::Int #Agent - Buyer
    costNw::Float64
    function Trade(tradeID, kind, Pes, Peb, As, Ab, costNw)
        w = new(tradeID, kind, Pes, Peb, As, Ab, costNw)
        return w
    end
end

type Agent
    agentID::Int
    location::Int
    kind::Int #1-Prosumer 2-Seller 3-Buyer
    pgMin::Float64
    pgMax::Float64
    costFn::Array{Float64,1}
    Pd::Float64
    Qd::Float64
    function Agent(agentID, location, kind, pgMin, pgMax, costFn, Pd, Qd)
        a = new(agentID, location, kind, pgMin, pgMax, costFn, Pd, Qd)
        return a
    end
end


text:
@testset "Rounding" begin
    f0,f1,f2 = 1.4,1.5,2.5

    @testset for T in (Posit8, Posit16, Posit32)
        @test round(f0) == Float64(round(T(f0)))
        @test round(f1) == Float64(round(T(f1)))
        @test round(f2) == Float64(round(T(f2)))
    end
end

@testset "Floor" begin

    fs = [-1.9,-1.5,-1.1,-1.0,1.0,1.1,1.5,1.9]

    @testset for T in (Posit8, Posit16, Posit32)
        for f in fs
            @test floor(f) == Float64(floor(T(f)))
            @test ceil(f) == Float64(ceil(T(f)))
        end
    end
end


text:
# testing a linear operator

function test_linearmap(F)
	for FF in [F, transpose(F)]
		# testing if (F) is linear
		x1=randn(size(FF,2))
		x2=randn(size(FF,2))
		x12=x1.+x2
		d12=FF*x12
		d1=FF*x1
		d2=FF*x2
		d12new=d1.+d2
		@test Misfits.error_squared_euclidean!(nothing, d12, d12new, nothing, norm_flag=true)<1e-25
	end

	adjtest(F)

	return nothing
end

function adjtest(F)
	x=randn(size(F,2))
	y=randn(size(F,1))
	a=LinearAlgebra.dot(y,F*x)
	b=LinearAlgebra.dot(x,adjoint(F)*y)
	c=LinearAlgebra.dot(x, transpose(F)*F*x)
	println("must be positive: ", c)
	println("adjoint test: ", a, "\t", b)
	@test isapprox(a,b,rtol=1e-5)
	@test c>0.0
	return nothing
end


text:
"""
Compute the partial density at the indicated ``k``-Point and return it (in Fourier space).
"""
function compute_partial_density!(ρ, basis, kpt, ψk, occupation)
    @assert length(occupation) == size(ψk, 2)

    # Build the partial density ρk_real for this k-Point
    ρk_real = [zeros(eltype(basis), basis.fft_size) for it = 1:Threads.nthreads()]
    ψnk_real = [zeros(complex(eltype(basis)), basis.fft_size) for it = 1:Threads.nthreads()]
    Threads.@threads for n = 1:size(ψk, 2)
        ψnk = @views ψk[:, n]
        tid = Threads.threadid()
        G_to_r!(ψnk_real[tid], basis, kpt, ψnk)
        ρk_real[tid] .+= occupation[n] .* abs2.(ψnk_real[tid])
    end
    for it = 2:Threads.nthreads()
        ρk_real[1] .+= ρk_real[it]
    end
    ρk_real = ρk_real[1]

    # Check sanity of the density (real, positive and normalized)
    T = real(eltype(ρk_real))
    real_checked(ρk_real)
    if all(occupation .> 0)
        minimum(real(ρk_real)) < 0 && @warn("Negative ρ detected",
                                            min_ρ=minimum(real(ρk_real)))
    end
    n_electrons = sum(ρk_real) * basis.model.unit_cell_volume / prod(basis.fft_size)
    if abs(n_electrons - sum(occupation)) > sqrt(eps(T))
        @warn("Mismatch in number of electrons", sum_ρ=n_electrons,
              sum_occupation=sum(occupation))
    end

    # FFT and return
    r_to_G!(ρ, basis, complex(ρk_real))
end

"""
    compute_density(basis::PlaneWaveBasis, ψ::AbstractVector, occupation::AbstractVector)

Compute the density and spin density for a wave function `ψ` discretized on the plane-wave
grid `basis`, where the individual k-Points are occupied according to `occupation`.
`ψ` should be one coefficient matrix per k-Point. If the `Model` underlying the basis
is not collinear the spin density is `nothing`.
"""
@timing function compute_density(basis::PlaneWaveBasis, ψ, occupation)
    n_k = length(basis.kpoints)
    n_spin = basis.model.n_spin_components

    # Sanity checks
    @assert n_k == length(ψ)
    @assert n_k == length(occupation)
    for ik in 1:n_k
        @assert length(G_vectors(basis.kpoints[ik])) == size(ψ[ik], 1)
        @assert length(occupation[ik]) == size(ψ[ik], 2)
    end
    @assert n_k > 0

    # Allocate an accumulator for ρ in each thread for each spin component
    ρaccus = [[similar(view(ψ[1], :, 1), basis.fft_size) for iσ in 1:n_spin]
              for ithread in 1:Threads.nthreads()]

    # TODO Better load balancing ... the workload per kpoint depends also on
    #      the number of symmetry operations. We know heuristically that the Gamma
    #      point (first k-Point) has least symmetry operations, so we will put
    #      some extra workload there if things do not break even
    kpt_per_thread = [ifelse(i <= n_k, [i], Vector{Int}()) for i in 1:Threads.nthreads()]
    if n_k >= Threads.nthreads()
        kblock = floor(Int, length(basis.kpoints) / Threads.nthreads())
        kpt_per_thread = [collect(1:length(basis.kpoints) - (Threads.nthreads() - 1) * kblock)]
        for ithread in 2:Threads.nthreads()
            push!(kpt_per_thread, kpt_per_thread[end][end] .+ collect(1:kblock))
        end
        @assert kpt_per_thread[end][end] == length(basis.kpoints)
    end

    Threads.@threads for (ikpts, ρaccu) in collect(zip(kpt_per_thread, ρaccus))
        for iσ in 1:n_spin
            ρaccu[iσ] .= 0
        end
        ρ_k = similar(ψ[1][:, 1], basis.fft_size)
        for ik in ikpts
            kpt = basis.kpoints[ik]
            compute_partial_density!(ρ_k, basis, kpt, ψ[ik], occupation[ik])
            lowpass_for_symmetry!(ρ_k, basis)
            # accumulates all the symops of ρ_k into ρaccu
            accumulate_over_symmetries!(ρaccu[kpt.spin], ρ_k, basis, basis.ksymops[ik])
        end
    end

    # Count the number of k-points modulo spin
    count = sum(length(basis.ksymops[ik]) for ik in 1:length(basis.kpoints)) ÷ n_spin
    count = mpi_sum(count, basis.comm_kpts)
    ρs = [sum(getindex.(ρaccus, iσ)) / count for iσ in 1:n_spin]
    mpi_sum!.(ρs, Ref(basis.comm_kpts))

    @assert basis.model.spin_polarization in (:none, :spinless, :collinear)
    if basis.model.spin_polarization == :collinear
        ρtot  = from_fourier(basis, ρs[1] + ρs[2])  # up + down
        ρspin = from_fourier(basis, ρs[1] - ρs[2])  # up - down
    else
        ρtot  = from_fourier(basis, ρs[1])
        ρspin = nothing
    end

    ρtot, ρspin
end


text:
module RegistryActivity

# Note about the statistics of the General registry: about 15k versions and 600
# packages were removed in October 2019 with
# https://github.com/JuliaRegistries/General/pull/4169.  It's complicated to
# deal with this because some packages had only some versions removed, others
# got completely removed and then they came back into the registry later.

using Pkg, Dates, LibGit2, TOML

export registry_activity

"""
    general_registry() -> String

Guess the path of the General registry.
"""
general_registry() =
    first(joinpath(d, "registries", "General") for d in Pkg.depots() if isfile(joinpath(d, "registries", "General", "Registry.toml")))

"""
    clone_registry(; origin="https://github.com/JuliaRegistries/General.git", dest=mktempdir())

Do a git clone of `origin` to `dest`.
"""
function clone_registry(;
                        origin="https://github.com/JuliaRegistries/General.git",
                        dest=mktempdir(),
                        )
    LibGit2.clone(origin, dest)
    return dest
end

"""
    get_date(c)

Get the UTC date of the Git commit `c`.
"""
function get_date(c)
    committer = LibGit2.committer(c)
    return Date(unix2datetime(committer.time + 60 * committer.time_offset))
end

"""
    extrema_dates(path)

Get the first and last date of commits in the git repository at `path`.
"""
function extrema_dates(path)
    repo = GitRepo(path)
    # Find the dates of all commits in the repo
    dates = LibGit2.with(LibGit2.GitRevWalker(repo)) do walker
        LibGit2.map((oid, repo) -> get_date(LibGit2.GitCommit(repo, oid)), walker, by=LibGit2.Consts.SORT_TIME)
    end
    return extrema(dates)
end

"""
    months_list(path, start_month, end_month)

Get the range of months when there have been commit to the Git repository at
`path`, between `start_month` and `end_month`.
"""
function months_list(path, start_month, end_month)
    min_month, max_month = extrema_dates(path)
    min_month = max(Date(year(min_month), month(min_month)), start_month)
    max_month = min(Date(year(max_month), month(max_month)), end_month)
    return min_month:Month(1):max_month
end

"""
    checkout_date(path, date::Date; branch="master")

Check out `branch` of Git repository at `path` on `date`.
"""
function checkout_date(path, date::Date;
                       branch="master",
                       )
    revision = readchomp(`git -C $(path) rev-list -n 1 --first-parent --before="$(date) 00:00:00 +0000" $(branch)`)
    LibGit2.checkout!(GitRepo(path), revision)
end

"""
    count_packages(path, packages=TOML.parsefile(joinpath(path, "Registry.toml"))["packages"])

Count the number of packages in the registry at `path`, using the list of
`packages`.
"""
function count_packages(path,
                        packages=TOML.parsefile(joinpath(path, "Registry.toml"))["packages"])
    return length(packages)
end

"""
    _non_yanked_versions(file)

Get the list of non-yanked versions of the package for which `file` is the path
to its `Versions.toml`.
"""
function _non_yanked_versions(file)
    toml = TOML.parsefile(file)
    collect(k for (k, v) in toml if !(haskey(v, "yanked") && v["yanked"]))
end

"""
    count_versions(path, packages=TOML.parsefile(joinpath(path, "Registry.toml"))["packages"])

Count the number of versions in the registry at `path`, using the list of
`packages`.
"""
function count_versions(path,
                        packages=TOML.parsefile(joinpath(path, "Registry.toml"))["packages"])
    paths = joinpath.(path, getindex.(values(packages), "path"))
    isempty(paths) && return 0
    return sum(isfile(joinpath(p, "Versions.toml")) ? length(_non_yanked_versions(joinpath(p, "Versions.toml"))) : 0 for p in paths)
end

"""
    registry_activity(path=clone_registry();
                      branch="master",
                      start_month=Date(2018, 2),
                      end_month=Date(year(now(UTC)), month(now(UTC))),
                      filter=p->!(occursin(r"_jll\$", last(p)["name"]) || last(p)["name"]=="julia"),
                      )

Get the 3-tuple `months, packages, versions` of vector of months, number of
packages, number of versions for the registry at `path` (by default a Git clone
of the General registry with [`clone_registry()`](@ref)), checking out the Git
`branch`.  The search is restricted between `start_month` and `end_month`.  You
can filter the packages matching the `filter` condition, which is a function
taking as argument an item in the dictionary of a parsed `Registry.toml`, that
is a pair whose key is the UUID of the package and the value is a dictionary
with information about the package, usually including the name and the path to
the directory where the package is recorded.
"""
function registry_activity(path=clone_registry();
                           branch="master",
                           # By default start from the first month there has
                           # been some actual activity in General
                           start_month=Date(2018, 2),
                           end_month=Date(year(now(UTC)), month(now(UTC))),
                           # By default exclude JLL packages and the virtual "julia" package
                           filter=p->!(occursin(r"_jll$", last(p)["name"]) || last(p)["name"]=="julia"),
                           )
    repo = LibGit2.GitRepo(path)
    commit = LibGit2.GitHash(LibGit2.lookup_branch(repo, branch))
    LibGit2.checkout!(repo, string(commit))
    months = months_list(path, start_month, end_month)
    packages = Vector{Int}(undef, length(months))
    versions = Vector{Int}(undef, length(months))
    for idx in eachindex(months)
        # Check out the registry at midnight of the first day of the month after
        checkout_date(path, months[idx] + Month(1))
        packages_list = if isfile(joinpath(path, "registry.toml"))
            # The registry file used to be call "registry.toml" at some point.
            TOML.parsefile(joinpath(path, "registry.toml"))["packages"]
        else
            TOML.parsefile(joinpath(path, "Registry.toml"))["packages"]
        end
        filter!(filter, packages_list)
        packages[idx] = count_packages(path, packages_list)
        versions[idx] = count_versions(path, packages_list)
    end
    # Restore target commit.
    LibGit2.checkout!(repo, string(commit))
    return months, packages, versions
end

end # module


text:
#holds the main parameter object
mutable struct Parameters
  σ²G::Float64
  ζ²G::Float64
  ζ²P::Float64
  SG::Float64
  SGP::Float64
end

#named constructor for testing purposes
function Parameters(;  σ²G::Float64 = 0.0,  ζ²G::Float64 = 0.0, ζ²P::Float64 = 0.0,
  SG::Float64 = 0.0,  SGP::Float64 = 0.0)

  return Parameters(σ²G, ζ²G, ζ²P, SG, SGP)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA