using Pkg
Pkg.activate("/app/environment")

# Handle dynamic environment fallback
if !haskey(Pkg.project().dependencies, "Plots")
    Pkg.add("Plots")
end

ENV["GKSwstype"] = "100" # Headless rendering for Docker container
using Plots

# Load framework core dependencies smoothly
juliamsi_path = "/app/JuliaMSI"
push!(LOAD_PATH, joinpath(juliamsi_path, "src"))
include(joinpath(juliamsi_path, "src", "MSI_src.jl"))
using .MSI_src

using CSV, DataFrames
using Statistics
using Plots.Measures

# 1. PATH CONFIGURATIONS
data_path = "/app/data/Leaf.imzML"
mask_path = "/app/data/region_of_interest.png"
clusters_path = "/app/msi_clusters_results.csv"

println("Loading data arrays...")
data = OpenMSIData(data_path)
clusters_df = CSV.read(clusters_path, DataFrame)

# 2. RUN FRAMEWORK ANALYTICS TO DISCOVER PHYSICAL BOUNDARIES
println("Syncing Dataset Precomputed Boundaries...")
MSI_src.precompute_analytics(data)

# Extract binned spectrum arrays
mzs_avg, ints_avg = MSI_src.get_average_spectrum(data, mask_path=mask_path)

# --- DYNAMIC AXIS SAFEGUARD ---
# Query the atomic fields populated by precompute_analytics
# fallback to standard values if the fields return uninitialized values
min_vis_mz = data.global_min_mz[] > 0.0 ? data.global_min_mz[] : 50.0
max_vis_mz = data.global_max_mz[] > 0.0 ? data.global_max_mz[] : 700.0

println("  Detected Dynamic Min m/z Limit: ", round(min_vis_mz, digits=2))
println("  Detected Dynamic Max m/z Limit: ", round(max_vis_mz, digits=2))

# 3. Configure Base Plot Layer with Dynamic X-limits
p = plot(mzs_avg, ints_avg, 
    linecolor = :lightgray,  # Light base so cluster sticks pop out visually
    linewidth = 0.7,
    label = "Average Spectrum (Reference)",
    xlabel = "m/z",
    ylabel = "Relative Intensity",
    title = "Cluster Mapping on Isotopic Pattern",
    size = (1200, 600),
    legend = :outerright,
    grid = false,       
    framestyle = :box,  
    xlims = (min_vis_mz, max_vis_mz), # Dynamically locks view strictly to valid mass range
    left_margin = 15mm,
    bottom_margin = 15mm
)

# 4. Use the exact qualitative 10-color palette seen in your image
cluster_colors = palette(:tab10)

# 5. Map Cluster Vectors to Binned Positions
println("Mapping chemical clusters onto binned coordinates...")
for cluster_id in sort(unique(clusters_df.cluster))
    sub_df = filter(row -> row.cluster == cluster_id, clusters_df)
    label_ready = false
    
    for target_mz in sub_df.mz
        # Find the closest matching bin channel index in the 2000-bin array
        idx = argmin(abs.(mzs_avg .- target_mz))
        
        m = mzs_avg[idx]
        i = ints_avg[idx]
        
        # Draw the vertical peak highlighter over the base spectrum
        plot!(p, [m, m], [0, i],
            linecolor = cluster_colors[cluster_id + 1],
            linewidth = 1.2,
            alpha = 0.9,
            label = label_ready ? "" : "Cluster $cluster_id"
        )
        label_ready = true
    end
end

# 6. Save assets to shared folder mount
savefig("/app/espectro_msi_paper_style.png")
println("Saved publication-ready plot: /app/espectro_msi_paper_style.png")
savefig("/app/espectro_msi_clusters_julia.png")
println("Saved backup copy: /app/espectro_msi_clusters_julia.png")