# using Revise 
using Flower
using YAML
using Printf

# ---------------------------------------------------------------------------------------------
# FUNCTION BARRIER: Wrapper pour exécuter la simulation
# Cela permet d'isoler le code compilé des variables dynamiques (BC) générées par eval()
# ---------------------------------------------------------------------------------------------

function execute_simulation_step(num, gp, gu, gv, op, phS, phL, sim, phys, time_scheme, bc_dict)
    
    # Préparation des arguments optionnels (Keyword Arguments)
    # On construit un NamedTuple uniquement avec les BC qui existent vraiment
    kwargs_bc = Dict{Symbol, Any}()
    
    # Liste des clés possibles
    keys_to_check = [:BC_uL, :BC_uS, :BC_vL, :BC_vS, :BC_pL, :BC_pS, :BC_u, :BC_int, :BC_trans_scal, :BC_phi_ele]
    
    for k in keys_to_check
        if haskey(bc_dict, k)
            kwargs_bc[k] = bc_dict[k]
        end
    end

    # Appel avec le splatting operator (...) qui déballe le dictionnaire en arguments
    run_forward!(
        num, gp, gu, gv, op, phS, phL;
        periodic_x = (sim.periodic_x == 1),
        periodic_y = (sim.periodic_y == 1),
        auto_reinit = sim.auto_reinit,
        time_scheme = time_scheme,
        electrolysis = true,
        navier_stokes = true,
        ns_advection = (sim.ns_advection == 1),
        ns_liquid_phase = (sim.solve_Navier_Stokes_liquid_phase == 1),
        verbose = true,
        show_every = sim.show_every,
        electrolysis_convection = (sim.electrolysis_convection == 1),
        electrolysis_liquid_phase = true,
        electrolysis_phase_change_case = sim.electrolysis_phase_change_case,
        imposed_velocity = sim.imposed_velocity,
        adapt_timestep_mode = sim.adapt_timestep_mode,
        non_dimensionalize = sim.non_dimensionalize,
        mode_2d = sim.mode_2d,
        breakup = sim.breakup,
        # On injecte ici les BC dynamiques présentes
        kwargs_bc...
    )
end

# Lecture des arguments et du fichier YAML
localARGS = ARGS
if length(localARGS) > 0
    yamlfile = localARGS[1]
    printstyled(color=:magenta, @sprintf "\n YAML file ")
    print(yamlfile)
    print("\n")
    yamlpath = yamlfile
    yamlnamelist = split(yamlfile, "/")
    yamlname = yamlnamelist[end]
end

data = YAML.load_file(yamlpath)

# Dictionnaires de configuration
prop_dict = PropertyDict(data)
study = PropertyDict(prop_dict.study)

# Paramètres de convergence
nb_grid_points = study.meshes
n_cases = length(nb_grid_points)
print("\n number of points ", nb_grid_points, "\n")
timesteps = study.timesteps

io = PropertyDict(prop_dict.plot) 
flower = PropertyDict(prop_dict.flower)
mesh = PropertyDict(flower.mesh)
sim = PropertyDict(flower.simulation)
phys = PropertyDict(flower.physics)
macros = PropertyDict(flower.macros) 

# Initialisation des variables globales via macro
eval(Meta.parseall(macros.print_parameters))

# ----------------------------------------------------------------
# Initialisation PDI (IO)
# ----------------------------------------------------------------
if io.pdi > 0
    @debug "Before PDI init"
    yml_file = yamlfile
    conf = @ccall "libparaconf".PC_parse_path(yml_file::Cstring)::PC_tree_t
    getsubyml = @ccall "libparaconf".PC_get(conf::PC_tree_t, ".pdi"::Cstring)::PC_tree_t  
    local pdi_status = @ccall "libpdi".PDI_init(getsubyml::PC_tree_t)::Cint
    
    # Metadonnées MPI
    mpi_coords_x = 1; mpi_coords_y = 1
    mpi_max_coords_x = 1; mpi_max_coords_y = 1
    local nx = 32
    local ny = 32
    nstep = 0
    phys_time = 0.0

    local PDI_status = @ccall "libpdi".PDI_multi_expose("init_PDI"::Cstring, 
            "mpi_coords_x"::Cstring, mpi_coords_x::Ref{Clonglong}, PDI_OUT::Cint,
            "mpi_coords_y"::Cstring, mpi_coords_x::Ref{Clonglong}, PDI_OUT::Cint,
            "mpi_max_coords_x"::Cstring, mpi_max_coords_x::Ref{Clonglong}, PDI_OUT::Cint,
            "mpi_max_coords_y"::Cstring, mpi_max_coords_y::Ref{Clonglong}, PDI_OUT::Cint,
            "nx"::Cstring, nx::Ref{Clonglong}, PDI_OUT::Cint,
            "ny"::Cstring, ny::Ref{Clonglong}, PDI_OUT::Cint,
            "nb_transported_scalars"::Cstring, phys.nb_transported_scalars::Ref{Clonglong}, PDI_OUT::Cint,
            "nb_levelsets"::Cstring, phys.nb_levelsets::Ref{Clonglong}, PDI_OUT::Cint,
            "nstep"::Cstring, nstep::Ref{Clonglong}, PDI_OUT::Cint,
            C_NULL::Ptr{Cvoid})::Cint

    print("\n PDI INIT status ", PDI_status)
end

# Tableaux pour stocker les erreurs
error_list_l1 = zeros(n_cases)
error_list_l2 = zeros(n_cases)
error_list_linfty = zeros(n_cases)
error_list_l1_mixed = zeros(n_cases)
error_list_l2_mixed = zeros(n_cases)
error_list_linfty_mixed = zeros(n_cases)
error_list_l1_full = zeros(n_cases)
error_list_l2_full = zeros(n_cases)
error_list_linfty_full = zeros(n_cases)

cell_volume_list = zeros(n_cases)
base_directory = pwd()

# ----------------------------------------------------------------
# Boucle de Convergence : Timestep
# ----------------------------------------------------------------
for timestep in timesteps
    
    print("\n Timestep ", timestep)
    timestep_to_string = @sprintf "timestep_%.4e" timestep
    timestep_to_string = replace(timestep_to_string, "." => "_")

    mkpath(timestep_to_string)
    cd(timestep_to_string)

    # ------------------------------------------------------------
    # Boucle de Convergence : Mesh
    # ------------------------------------------------------------
    for (i, study_nb_grid_points) in enumerate(nb_grid_points)
        
        mesh_to_string = @sprintf "mesh_%.5i" study_nb_grid_points
        mesh_ratio = Int((mesh.ymax - mesh.ymin)/ (mesh.xmax - mesh.xmin))
        print("\n mesh ratio ", mesh_ratio) 

        study_nb_grid_points_x = study_nb_grid_points
        study_nb_grid_points_y = study_nb_grid_points * mesh_ratio
        print("\n nx ", study_nb_grid_points_x, " ny ", study_nb_grid_points_y)

        mkpath(mesh_to_string)
        cd(mesh_to_string)

        # Init regular grid
        scalar_mesh_x = collect(LinRange(mesh.xmin, mesh.xmax, study_nb_grid_points_x + 1))    
        scalar_mesh_y = collect(LinRange(mesh.ymin, mesh.ymax, study_nb_grid_points_y + 1))

        @debug "Before Numerical"

        # Aliases from yml to Flower.jl
        aliases = Dict(
            :x => :scalar_mesh_x,
            :y => :scalar_mesh_y,
            :xcoord => :intfc_x,
            :ycoord => :intfc_y,
            :R => :radius,
            :max_iterations => :max_iter,
            :save_every => :max_iter,
            :ϵ => :epsilon,
            :ϵwall => :epsilon_wall,
            :nLS => :nb_levelsets,
            :θd => :temperature0,
            :β => :beta,
            :σ => :sigma,
            :δreinit => :delta_reinit,
            :n_ext_cl => :n_ext,
            :timestep_0 => :timestep,
            :timestep_n => :timestep,
            :io_pdi => :pdi,
        )

        extra = Dict(
            :scalar_mesh_x => scalar_mesh_x,
            :scalar_mesh_y => scalar_mesh_y,
        )

        default = Numerical{Float64,Int}(
            x = scalar_mesh_x,
            y = scalar_mesh_y,
            timestep_n = timestep,
            timestep_0 = timestep,
            electrolysis_reaction_symb = Symbol(phys.electrolysis_reaction),
            bulk_velocity_symb = Symbol(phys.bulk_velocity),
        )

        print("\n ns_advection ", (sim.ns_advection == 1))

        global num_new = safefill_with_aliases_and_extra_already_init(Numerical{Float64,Int}, default,
                                 sim, phys, io, aliases, extra)

        global num = num_new

        if num.verbosity > 0
            print("\n Parameters num ", num)
        end

        Broadcast.broadcastable(num::Numerical) = Ref(num) 
        
        global gp, gu, gv = init_meshes(num)
        global op, phS, phL = init_fields(num, gp, gu, gv)

        gp.LS[1].u .= 1.0 # deactivate interface

		
		# ------------------------------------------------------------
        # EVALUATION DYNAMIQUE (Macros & BC)
        # ------------------------------------------------------------
        # 1. On exécute TOUTES les macros pour créer les variables dans Main
        eval(Meta.parseall(macros.init_fields))
        eval(Meta.parseall(macros.boundaries))
        eval(Meta.parseall(macros.interface)) # <--- DÉPLACÉ ICI (IMPORTANT)

        # 2. Capture des BC dans un dictionnaire pour éviter le "World Age Problem"
        bc_dict = Dict{Symbol, Any}()
        bc_symbols = [:BC_uL, :BC_uS, :BC_vL, :BC_vS, :BC_pL, :BC_pS, :BC_u, :BC_int, :BC_trans_scal, :BC_phi_ele]
        
        for sym in bc_symbols
            if isdefined(Main, sym)
                bc_dict[sym] = getfield(Main, sym)
            else
                # Si la variable n'est pas définie, on ne l'ajoute pas au dict.
                # Cela laissera run_forward! utiliser sa valeur par défaut (si elle existe)
                # ou plantera avec une erreur plus claire si elle est obligatoire.
            end
        end
        # ------------------------------------------------------------

        if num.io_pdi > 0
            # Send meta-data to PDI
            nx = gp.nx
            ny = gp.ny
            try
                local PDI_status = @ccall "libpdi".PDI_multi_expose("init_PDI"::Cstring, 
                        "mpi_coords_x"::Cstring, mpi_coords_x::Ref{Clonglong}, PDI_OUT::Cint,
                        "mpi_coords_y"::Cstring, mpi_coords_x::Ref{Clonglong}, PDI_OUT::Cint,
                        "mpi_max_coords_x"::Cstring, mpi_max_coords_x::Ref{Clonglong}, PDI_OUT::Cint,
                        "mpi_max_coords_y"::Cstring, mpi_max_coords_y::Ref{Clonglong}, PDI_OUT::Cint,
                        "nx"::Cstring, nx::Ref{Clonglong}, PDI_OUT::Cint,
                        "ny"::Cstring, ny::Ref{Clonglong}, PDI_OUT::Cint,
                        "nb_transported_scalars"::Cstring, phys.nb_transported_scalars::Ref{Clonglong}, PDI_OUT::Cint,
                        "nb_levelsets"::Cstring, phys.nb_levelsets::Ref{Clonglong}, PDI_OUT::Cint,
                        "nstep"::Cstring, num.current_iter::Ref{Clonglong}, PDI_OUT::Cint,
                        "nb_Navier_slip_BC"::Cstring, num.nNavier::Ref{Clonglong}, PDI_OUT::Cint,
                        "timestep"::Cstring, timestep::Ref{Cdouble}, PDI_OUT::Cint,
                        C_NULL::Ptr{Cvoid})::Cint
            catch error
                print("\n Bug init_PDI \n")
                print("\n error", error)
            end
        end

        eval(Meta.parseall(macros.interface))

        if sim.time_scheme == "FE"
            time_scheme = FE
        else
            time_scheme = CN
        end

        # --- Sorties PDI Initiale ---
        if num.io_pdi > 0
            try
                iLSpdi = 1 
                LStable = zeros(gp)
                if phys.nb_levelsets > 1
                    LStable = gp.LS[2].u
                end
                
                us = zeros(gp)
                vs = zeros(gp)
                interpolate_grid_liquid!(gp, gu, gv, phL.u, phL.v, us, vs)
                current_radius = phys.radius

                PDI_status = @ccall "libpdi".PDI_multi_expose("write_initialization"::Cstring,
                    "nstep"::Cstring, num.current_iter::Ref{Clonglong}, PDI_OUT::Cint,
                    "time"::Cstring, phys_time::Ref{Cdouble}, PDI_OUT::Cint,
                    "u_1D"::Cstring, phL.uD::Ptr{Cdouble}, PDI_OUT::Cint,
                    "v_1D"::Cstring, phL.vD::Ptr{Cdouble}, PDI_OUT::Cint,
                    "p_1D"::Cstring, phL.pD::Ptr{Cdouble}, PDI_OUT::Cint,
                    "levelset_p"::Cstring, gp.LS[iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
                    "levelset_u"::Cstring, gu.LS[iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
                    "levelset_v"::Cstring, gv.LS[iLSpdi].u::Ptr{Cdouble}, PDI_OUT::Cint,
                    "levelset_p_wall"::Cstring, LStable::Ptr{Cdouble}, PDI_OUT::Cint,
                    "phi_ele_1D"::Cstring, phL.phi_eleD::Ptr{Cdouble}, PDI_OUT::Cint,   
                    "velocity_x"::Cstring, us::Ptr{Cdouble}, PDI_OUT::Cint,   
                    "velocity_y"::Cstring, vs::Ptr{Cdouble}, PDI_OUT::Cint,      
                    "radius"::Cstring, current_radius::Ref{Cdouble}, PDI_OUT::Cint,  
                    C_NULL::Ptr{Cvoid})::Cint
            catch error
                printstyled(color=:red, @sprintf "\n PDI error \n")
                print(error)
            end
        end 

        # --- Sorties PDI Post-process bubble ---
        velocity_y_bubble = 0.0 .* gp.LS[end].geoS.dcap[:,:,5]
        volume_fraction = gp.LS[end].geoS.dcap[:,:,5]
        iLSpdi = 1
        velocity_y = velocity_y_bubble

        PDI_status = @ccall "libpdi".PDI_multi_expose("post_processing_rising_bubble"::Cstring,
            "nstep"::Cstring, num.current_iter ::Ref{Clonglong}, PDI_OUT::Cint,
            "velocity_y"::Cstring, velocity_y::Ptr{Cdouble}, PDI_OUT::Cint,      
            "volume_fraction"::Cstring, volume_fraction::Ptr{Cdouble}, PDI_OUT::Cint,
            "volume_liq_cell"::Cstring, gp.LS[iLSpdi].geoL.dcap[:,:,5]::Ptr{Cdouble}, PDI_OUT::Cint, 
            "volume_cell"::Cstring, gp.LS[iLSpdi].geoS.dcap[:,:,5]::Ptr{Cdouble}, PDI_OUT::Cint, 
            "mesh_p_x"::Cstring, gp.x::Ptr{Cdouble}, PDI_OUT::Cint,
            "mesh_p_y"::Cstring, gp.y::Ptr{Cdouble}, PDI_OUT::Cint,
            "dcap_1"::Cstring, gp.LS[iLSpdi].geoS.dcap[:,:,1]::Ptr{Cdouble}, PDI_OUT::Cint, 
            "dcap_2"::Cstring, gp.LS[iLSpdi].geoS.dcap[:,:,2]::Ptr{Cdouble}, PDI_OUT::Cint, 
            "dcap_3"::Cstring, gp.LS[iLSpdi].geoS.dcap[:,:,3]::Ptr{Cdouble}, PDI_OUT::Cint, 
            "dcap_4"::Cstring, gp.LS[iLSpdi].geoS.dcap[:,:,4]::Ptr{Cdouble}, PDI_OUT::Cint, 
            C_NULL::Ptr{Cvoid})::Cint

        # -----------------------------------------------------------------------------------------
        # APPEL DU SOLVEUR via FUNCTION BARRIER + INVOKELATEST
        # C'est ici que la correction opère : on passe le dictionnaire 'bc_dict' via invokelatest
        # -----------------------------------------------------------------------------------------
        @debug "Before run_forward wrapper"
        Base.invokelatest(
            execute_simulation_step, 
            num, gp, gu, gv, op, phS, phL, sim, phys, time_scheme, bc_dict
        )
        @debug "After run"
        # -----------------------------------------------------------------------------------------

        # Calcul d'erreurs (si demandé)
        if study.compute_errors != "None"

            LIQUID = gp.ind.all_indices[gp.LS[1].geoL.cap[:,:,5] .> (1-1e-16)]
            MIXED = gp.ind.all_indices[gp.LS[1].geoL.cap[:,:,5] .<= (1-1e-16) .&& gp.LS[1].geoL.cap[:,:,5] .> 1e-16]

            if study.compute_errors == "Poiseuille"
                l1, l2, linfty = relative_errors(phL.v, vPoiseuille, vcat(LIQUID, MIXED), gp.LS[1].geoL.cap[:,:,5], num.Δ)
                l1_mixed, l2_mixed, linfty_mixed = relative_errors(phL.v, vPoiseuille, MIXED, gp.LS[1].geoL.cap[:,:,5], num.Δ)
                l1_full, l2_full, linfty_full = relative_errors(phL.v, vPoiseuille, LIQUID, gp.LS[1].geoL.cap[:,:,5], num.Δ)
            
            elseif study.compute_errors == "diffusion_full_cell"
                L = mesh.xmax - mesh.xmin
                D = num.diffusion_coeff[2]
                analytical_nx = 1000
                diffusion_time_scale = L^2 / D
                ele_current_i = -1.59e4
                F1 = ele_current_i / (2 * num.Faraday * D)
                analytical_dx = L / analytical_nx
                c0 = num.concentration0[2] 
                max_nb_Fourier_series = 1000

                function full_cell_u1(x, t)
                    return F1 * x
                end
                function full_cell_A_n(n, L, F1, dx, c0)
                    if n == 0
                        return -F1*L + 2*c0 
                    else
                        return (-2*F1*L)/(n*π)^2*((-1)^n-1) 
                    end
                end
                function full_cell_u2(x, t, L, D, max_nb_Fourier_series, c0, analytical_dx)
                    result = sum(full_cell_A_n(n, L, F1, analytical_dx, c0) * cos(n * π * x / L) * exp(-D * (n * π / L)^2 * t) for n in 1:max_nb_Fourier_series)
                    result += full_cell_A_n(0, L, F1, analytical_dx, c0) * cos(0 * π * x / L) * exp(-D * (0 * π / L)^2 * t) / 2
                    return result
                end
                function full_cell_u(x, t, L, D, max_nb_Fourier_series, c0, analytical_dx)
                    return full_cell_u1(x, t) + full_cell_u2(x, t, L, D, max_nb_Fourier_series, c0, analytical_dx)
                end
                
                concentration_profile = full_cell_u.(gp.x[2,:], num.time, L, num.diffusion_coeff[2], max_nb_Fourier_series, c0, analytical_dx)
                
                l1, l2, linfty = relative_errors(phL.trans_scal[:,:,2], concentration_profile, vcat(LIQUID, MIXED), gp.LS[1].geoL.cap[:,:,5], num.Δ)
                l1_mixed, l2_mixed, linfty_mixed = relative_errors(phL.trans_scal[:,:,2], concentration_profile, MIXED, gp.LS[1].geoL.cap[:,:,5], num.Δ)
                l1_full, l2_full, linfty_full = relative_errors(phL.trans_scal[:,:,2], concentration_profile, LIQUID, gp.LS[1].geoL.cap[:,:,5], num.Δ)
            end

            error_list_l1[i] = l1
            error_list_l2[i] = l2
            error_list_linfty[i] = linfty
            error_list_l1_mixed[i] = l1_mixed
            error_list_l2_mixed[i] = l2_mixed
            error_list_linfty_mixed[i] = linfty_mixed
            error_list_l1_full[i] = l1_full
            error_list_l2_full[i] = l2_full
            error_list_linfty_full[i] = linfty_full
            cell_volume_list[i] = minimum(gp.LS[1].geoL.dcap[:,:,5])
            min_cell_volume = minimum(gp.LS[1].geoL.dcap[:,:,5])

            local PDI_status = @ccall "libpdi".PDI_multi_expose("convergence_study_iter"::Cstring, 
                "study_nb_grid_points"::Cstring, study_nb_grid_points::Ref{Clong}, PDI_OUT::Cint,
                "study_timestep"::Cstring, timestep::Ref{Cdouble}, PDI_OUT::Cint,
                "study_l1_rel_error"::Cstring, l1::Ref{Cdouble}, PDI_OUT::Cint,
                "study_l2_rel_error"::Cstring, l2::Ref{Cdouble}, PDI_OUT::Cint,
                "study_linfty_rel_error"::Cstring, linfty::Ref{Cdouble}, PDI_OUT::Cint,
                "study_l1_rel_error_full_cells"::Cstring, l1_full::Ref{Cdouble}, PDI_OUT::Cint,
                "study_l2_rel_error_full_cells"::Cstring, l2_full::Ref{Cdouble}, PDI_OUT::Cint,
                "study_linfty_rel_error_full_cells"::Cstring, linfty_full::Ref{Cdouble}, PDI_OUT::Cint,
                "study_l1_rel_error_partial_cells"::Cstring, l1_mixed::Ref{Cdouble}, PDI_OUT::Cint,
                "study_l2_rel_error_partial_cells"::Cstring, l2_mixed::Ref{Cdouble}, PDI_OUT::Cint,
                "study_linfty_rel_error_partial_cells"::Cstring, linfty_mixed::Ref{Cdouble}, PDI_OUT::Cint,
                "domain_length"::Cstring, L0::Ref{Cdouble}, PDI_OUT::Cint,
                "min_cell_volume"::Cstring, min_cell_volume::Ref{Cdouble}, PDI_OUT::Cint,
                C_NULL::Ptr{Cvoid})::Cint
        end

        current_directory = pwd()
        if current_directory == base_directory*"/"*timestep_to_string*"/"*mesh_to_string
            cd("..")
        else 
            print("\n Error in mesh subdirectory", current_directory, " not ", base_directory*"/"*timestep_to_string*"/"*mesh_to_string)
            exit()
        end

    end # end mesh convergence

    current_directory = pwd()
    if current_directory == base_directory*"/"*timestep_to_string
        cd("..")
    else
        print("\n Error in timestep subdirectory", current_directory, " not ", base_directory*"/"*timestep_to_string)
        exit()
    end

end # end timestep convergence

# Finalisation PDI et Tests
min_cell_volume = minimum(gp.LS[1].geoL.cap[:,:,5])
print("\n min_cell_volume ", min_cell_volume, " type ", typeof(min_cell_volume))

if study.compute_errors != "None"
    local PDI_status = @ccall "libpdi".PDI_multi_expose("convergence_study"::Cstring, 
    "n_tests"::Cstring, n_cases::Ref{Clonglong}, PDI_OUT::Cint,
    "nx_list"::Cstring, nb_grid_points::Ptr{Clonglong}, PDI_OUT::Cint,
    "cell_volume_list"::Cstring, cell_volume_list::Ptr{Cdouble}, PDI_OUT::Cint,
    "l1_rel_error"::Cstring, error_list_l1::Ptr{Cdouble}, PDI_OUT::Cint,
    "l2_rel_error"::Cstring, error_list_l2::Ptr{Cdouble}, PDI_OUT::Cint,
    "linfty_rel_error"::Cstring, error_list_linfty::Ptr{Cdouble}, PDI_OUT::Cint,
    "l1_rel_error_full_cells"::Cstring, error_list_l1_full::Ptr{Cdouble}, PDI_OUT::Cint,
    "l2_rel_error_full_cells"::Cstring, error_list_l2_full::Ptr{Cdouble}, PDI_OUT::Cint,
    "linfty_rel_error_full_cells"::Cstring, error_list_linfty_full::Ptr{Cdouble}, PDI_OUT::Cint,
    "l1_rel_error_partial_cells"::Cstring, error_list_l1_mixed::Ptr{Cdouble}, PDI_OUT::Cint,
    "l2_rel_error_partial_cells"::Cstring, error_list_l2_mixed::Ptr{Cdouble}, PDI_OUT::Cint,
    "linfty_rel_error_partial_cells"::Cstring, error_list_linfty_mixed::Ptr{Cdouble}, PDI_OUT::Cint,
    "domain_length"::Cstring, L0::Ref{Cdouble}, PDI_OUT::Cint,
    "min_cell_volume"::Cstring, min_cell_volume::Ref{Cdouble}, PDI_OUT::Cint,
    C_NULL::Ptr{Cvoid})::Cint
end
 
if io.pdi > 0
    try
        local PDI_status = @ccall "libpdi".PDI_finalize()::Cint
    catch error
        printstyled(color=:red, @sprintf "\n PDI error \n")
        print(error)
    end
end

printstyled(color=:red, @sprintf "\n After PDI \n")

if haskey(macros, "test_end")
    eval(Meta.parseall(macros.test_end))
end