From 8b6bed0ff52882095adfe823b396160a4686ba71 Mon Sep 17 00:00:00 2001 From: Pixelguy14 Date: Thu, 12 Mar 2026 12:51:16 -0600 Subject: [PATCH 1/6] Hot fix to the precalculation atomic variable accountability, and other base threads safety --- app.jl | 12 ++++++-- src/MSIData.jl | 66 +++++++++++++++++++++--------------------- src/Precalculations.jl | 16 +++++----- 3 files changed, 51 insertions(+), 43 deletions(-) diff --git a/app.jl b/app.jl index 2966aba..11058a5 100644 --- a/app.jl +++ b/app.jl @@ -1434,7 +1434,11 @@ end # --- 2. Parameter Assembly with Validation --- current_pipeline_step = "Configuring parameters..." println("DEBUG: Configuring parameters and validating enabled steps...") - ref_peaks = Dict{Float64, String}(p["mz"] => p["label"] for p in reference_peaks_list) + ref_peaks = Dict{Float64, String}( + parse(Float64, string(p["mz"])) => string(get(p, "label", "")) + for p in reference_peaks_list + if tryparse(Float64, string(p["mz"])) !== nothing + ) final_params = Dict{Symbol, Dict{Symbol, Any}}() validation_errors = String[] @@ -1944,7 +1948,11 @@ end try msg = "Recalculating suggestions..." - ref_peaks = Dict(p["mz"] => p["label"] for p in reference_peaks_list) + ref_peaks = Dict{Float64, String}( + parse(Float64, string(p["mz"])) => string(get(p, "label", "")) + for p in reference_peaks_list + if tryparse(Float64, string(p["mz"])) !== nothing + ) recommended_params = main_precalculation(msi_data, reference_peaks=ref_peaks) diff --git a/src/MSIData.jl b/src/MSIData.jl index 995370b..96a01ca 100644 --- a/src/MSIData.jl +++ b/src/MSIData.jl @@ -277,9 +277,9 @@ mutable struct MSIData # Buffer Pool for binary data operations buffer_pool::SimpleBufferPool - # Pre-computed analytics/metadata - use Threads.Atomic for compatibility - global_min_mz::Threads.Atomic{Float64} - global_max_mz::Threads.Atomic{Float64} + # Pre-computed analytics/metadata - use Base.Threads.Atomic for compatibility + global_min_mz::Base.Threads.Atomic{Float64} + global_max_mz::Base.Threads.Atomic{Float64} spectrum_stats_df::Union{DataFrame, Nothing} bloom_filters::Union{Vector{<:BloomFilter}, Nothing} analytics_ready::AtomicFlag @@ -289,7 +289,7 @@ mutable struct MSIData obj = new(source, metadata, instrument_meta, dims, coordinate_map, Dict(), [], cache_size, ReentrantLock(), SimpleBufferPool(), - Threads.Atomic{Float64}(Inf), Threads.Atomic{Float64}(-Inf), + Base.Threads.Atomic{Float64}(Inf), Base.Threads.Atomic{Float64}(-Inf), nothing, nothing, AtomicFlag(), nothing) # Ensure file handles are closed when the object is garbage collected @@ -319,8 +319,8 @@ Sets the global m/z range for the MSIData object. - `nothing` """ function set_global_mz_range!(data::MSIData, min_mz::Float64, max_mz::Float64) - Threads.atomic_xchg!(data.global_min_mz, min_mz) - Threads.atomic_xchg!(data.global_max_mz, max_mz) + Base.Threads.atomic_xchg!(data.global_min_mz, min_mz) + Base.Threads.atomic_xchg!(data.global_max_mz, max_mz) end """ @@ -337,8 +337,8 @@ Gets the global m/z range for the MSIData object. """ function get_global_mz_range(data::MSIData) # Atomic read using atomic_add! with zero - min_mz = Threads.atomic_add!(data.global_min_mz, 0.0) - max_mz = Threads.atomic_add!(data.global_max_mz, 0.0) + min_mz = Base.Threads.atomic_add!(data.global_min_mz, 0.0) + max_mz = Base.Threads.atomic_add!(data.global_max_mz, 0.0) return (min_mz, max_mz) end @@ -896,8 +896,8 @@ function precompute_analytics(msi_data::MSIData) num_spectra = length(msi_data.spectra_metadata) # Initialize thread-local variables for global stats - thread_local_min_mz = [Inf for _ in 1:Threads.nthreads()] - thread_local_max_mz = [-Inf for _ in 1:Threads.nthreads()] + thread_local_min_mz = [Inf for _ in 1:Base.Threads.nthreads()] + thread_local_max_mz = [-Inf for _ in 1:Base.Threads.nthreads()] # Initialize vectors for per-spectrum stats tics = Vector{Float64}(undef, num_spectra) @@ -937,7 +937,7 @@ function precompute_analytics(msi_data::MSIData) # Update thread-local global m/z range local_min, local_max = extrema(mz) - thread_id = Threads.threadid() + thread_id = Base.Threads.threadid() thread_local_min_mz[thread_id] = min(thread_local_min_mz[thread_id], local_min) thread_local_max_mz[thread_id] = max(thread_local_max_mz[thread_id], local_max) min_mzs[idx] = local_min @@ -963,16 +963,16 @@ function precompute_analytics(msi_data::MSIData) # Only update if we found valid ranges if isfinite(g_min_mz) && isfinite(g_max_mz) && g_min_mz < g_max_mz - Threads.atomic_xchg!(msi_data.global_min_mz, g_min_mz) - Threads.atomic_xchg!(msi_data.global_max_mz, g_max_mz) + Base.Threads.atomic_xchg!(msi_data.global_min_mz, g_min_mz) + Base.Threads.atomic_xchg!(msi_data.global_max_mz, g_max_mz) else # Fallback: calculate from first spectrum try if num_spectra > 0 mz, _ = GetSpectrum(msi_data, 1) if !isempty(mz) - Threads.atomic_xchg!(msi_data.global_min_mz, minimum(mz)) - Threads.atomic_xchg!(msi_data.global_max_mz, maximum(mz)) + Base.Threads.atomic_xchg!(msi_data.global_min_mz, minimum(mz)) + Base.Threads.atomic_xchg!(msi_data.global_max_mz, maximum(mz)) end end catch e @@ -1040,10 +1040,10 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_ local global_min_mz, global_max_mz - if Threads.atomic_add!(msi_data.global_min_mz, 0.0) !== Inf + if Base.Threads.atomic_add!(msi_data.global_min_mz, 0.0) !== Inf println(" Using pre-computed m/z range.") - global_min_mz = Threads.atomic_add!(msi_data.global_min_mz, 0.0) - global_max_mz = Threads.atomic_add!(msi_data.global_max_mz, 0.0) + global_min_mz = Base.Threads.atomic_add!(msi_data.global_min_mz, 0.0) + global_max_mz = Base.Threads.atomic_add!(msi_data.global_max_mz, 0.0) else # 1. First Pass: Find the global m/z range by reading the fast .ibd file pass1_start_time = time_ns() @@ -1077,14 +1077,14 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_ min_mz = global_min_mz # Initialize thread-local intensity sums - thread_local_intensity_sums = [zeros(Float64, num_bins) for _ in 1:Threads.nthreads()] - thread_local_spectra_processed = [0 for _ in 1:Threads.nthreads()] + thread_local_intensity_sums = [zeros(Float64, num_bins) for _ in 1:Base.Threads.nthreads()] + thread_local_spectra_processed = [0 for _ in 1:Base.Threads.nthreads()] # 3. Second Pass: Optimized binning pass2_start_time = time_ns() println(" Pass 2: Summing intensities into $num_bins bins...") _iterate_spectra_fast(msi_data, masked_indices === nothing ? nothing : collect(masked_indices)) do idx, mz, intensity - thread_id = Threads.threadid() + thread_id = Base.Threads.threadid() thread_local_spectra_processed[thread_id] += 1 if isempty(mz) return @@ -1109,7 +1109,7 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_ # Combine thread-local results intensity_sum = zeros(Float64, num_bins) num_spectra_processed = 0 - for i in 1:Threads.nthreads() + for i in 1:Base.Threads.nthreads() intensity_sum .+= thread_local_intensity_sums[i] num_spectra_processed += thread_local_spectra_processed[i] end @@ -1157,10 +1157,10 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000, masked_i local global_min_mz, global_max_mz - if Threads.atomic_add!(msi_data.global_min_mz, 0.0) !== Inf + if Base.Threads.atomic_add!(msi_data.global_min_mz, 0.0) !== Inf println(" Using pre-computed m/z range.") - global_min_mz = Threads.atomic_add!(msi_data.global_min_mz, 0.0) - global_max_mz = Threads.atomic_add!(msi_data.global_max_mz, 0.0) + global_min_mz = Base.Threads.atomic_add!(msi_data.global_min_mz, 0.0) + global_max_mz = Base.Threads.atomic_add!(msi_data.global_max_mz, 0.0) else # --- Pass 1: Find m/z range and cache it --- pass1_start_time = time_ns() @@ -1200,11 +1200,11 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000, masked_i min_mz = global_min_mz # Initialize thread-local intensity sums - thread_local_intensity_sums = [zeros(Float64, num_bins) for _ in 1:Threads.nthreads()] - thread_local_spectra_processed = [0 for _ in 1:Threads.nthreads()] + thread_local_intensity_sums = [zeros(Float64, num_bins) for _ in 1:Base.Threads.nthreads()] + thread_local_spectra_processed = [0 for _ in 1:Base.Threads.nthreads()] _iterate_spectra_fast(msi_data, masked_indices === nothing ? nothing : collect(masked_indices)) do idx, mz, intensity - thread_id = Threads.threadid() + thread_id = Base.Threads.threadid() thread_local_spectra_processed[thread_id] += 1 if isempty(mz) return @@ -1224,7 +1224,7 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000, masked_i # Combine thread-local results intensity_sum = zeros(Float64, num_bins) num_spectra_processed = 0 - for i in 1:Threads.nthreads() + for i in 1:Base.Threads.nthreads() intensity_sum .+= thread_local_intensity_sums[i] num_spectra_processed += thread_local_spectra_processed[i] end @@ -1722,7 +1722,7 @@ end A smart dispatcher for internal, high-performance iteration over spectra. It automatically chooses between serial and parallel execution based on the -number of available threads (`Threads.nthreads()`). +number of available threads (`Base.Threads.nthreads()`). # Arguments - `f`: A function to call for each spectrum, with the signature `f(index, mz, intensity)`. @@ -1745,7 +1745,7 @@ function _iterate_spectra_fast(f::Function, data::MSIData, indices_to_iterate::U # Default to all indices if not specified all_indices = (indices_to_iterate === nothing) ? (1:length(data.spectra_metadata)) : indices_to_iterate - if Threads.nthreads() > 1 && length(all_indices) > 100 # Heuristic: only parallelize for enough work + if Base.Threads.nthreads() > 1 && length(all_indices) > 100 # Heuristic: only parallelize for enough work _iterate_spectra_fast_parallel(f, data, all_indices) else _iterate_spectra_fast_serial(f, data, all_indices) @@ -1794,11 +1794,11 @@ Each thread gets its own file handle, eliminating contention. """ function _iterate_spectra_fast_parallel(f::Function, data::MSIData, indices::AbstractVector) # Split indices into chunks for each thread - n_chunks = Threads.nthreads() + n_chunks = Base.Threads.nthreads() chunk_size = ceil(Int, length(indices) / n_chunks) chunks = collect(Iterators.partition(indices, chunk_size)) - Threads.@threads for chunk in chunks + Base.Threads.@threads for chunk in chunks # Each thread gets its own file handle based on source type if data.source isa ImzMLSource local_handle = open(data.source.ibd_handle.path, "r") diff --git a/src/Precalculations.jl b/src/Precalculations.jl index a99b74f..786c4af 100644 --- a/src/Precalculations.jl +++ b/src/Precalculations.jl @@ -81,12 +81,12 @@ function analyze_mass_accuracy( println("Analyzing mass accuracy for $(length(spectrum_indices)) spectra...") all_ppm_errors = Float64[] - total_matched_peaks = Atomic{Int}(0) - total_spectra_processed = Atomic{Int}(0) + total_matched_peaks = Base.Threads.Atomic{Int}(0) + total_spectra_processed = Base.Threads.Atomic{Int}(0) results_lock = ReentrantLock() _iterate_spectra_fast(msi_data, spectrum_indices) do idx, mz, intensity - atomic_add!(total_spectra_processed, 1) + Base.Threads.atomic_add!(total_spectra_processed, 1) if !validate_spectrum(mz, intensity) @warn "Spectrum $idx is invalid, skipping mass accuracy analysis for it." return @@ -111,7 +111,7 @@ function analyze_mass_accuracy( lock(results_lock) do push!(all_ppm_errors, min_ppm_error) end - atomic_add!(total_matched_peaks, 1) + Base.Threads.atomic_add!(total_matched_peaks, 1) end end end @@ -671,8 +671,8 @@ function analyze_peak_characteristics(msi_data::MSIData, instrument_analysis::Di peak_counts = Int[] fwhm_values = Float64[] - spectra_analyzed = Atomic{Int}(0) - peaks_analyzed = Atomic{Int}(0) + spectra_analyzed = Base.Threads.Atomic{Int}(0) + peaks_analyzed = Base.Threads.Atomic{Int}(0) results_lock = ReentrantLock() _iterate_spectra_fast(msi_data, spectrum_indices) do idx, mz, intensity @@ -680,7 +680,7 @@ function analyze_peak_characteristics(msi_data::MSIData, instrument_analysis::Di return end - atomic_add!(spectra_analyzed, 1) + Base.Threads.atomic_add!(spectra_analyzed, 1) meta = msi_data.spectra_metadata[idx] # Detect peaks with lower SNR threshold to find more peaks @@ -709,7 +709,7 @@ function analyze_peak_characteristics(msi_data::MSIData, instrument_analysis::Di r2 = _fit_gaussian_and_r2(mz, intensity, peak_idx, 5) push!(r_squared_values, r2) end - atomic_add!(peaks_analyzed, 1) + Base.Threads.atomic_add!(peaks_analyzed, 1) if peaks_analyzed[] <= 3 println("DEBUG: Peak at m/z $(peak.mz), FWHM = $(fwhm_ppm) ppm, R² = $(r_squared_values[end])") -- 2.43.0 From e0ef601a9b61a89eedcc2755d0572e50676ce354 Mon Sep 17 00:00:00 2001 From: Pixelguy14 Date: Mon, 6 Apr 2026 15:41:38 -0600 Subject: [PATCH 2/6] hopefully fixed the JSON Lock for Windows --- julia_imzML_visual.jl | 162 +++++++++++++++++++----------------------- mask.jl | 4 +- 2 files changed, 75 insertions(+), 91 deletions(-) diff --git a/julia_imzML_visual.jl b/julia_imzML_visual.jl index 502e17d..33acb19 100644 --- a/julia_imzML_visual.jl +++ b/julia_imzML_visual.jl @@ -903,6 +903,18 @@ function warmup_init() end end +""" + clean_registry_path(registry_path::String) + +Sanitizes the file path to prevent OS/file system issues, particularly on Windows, +by standardizing separators and removing stripping whitespace. +""" +function clean_registry_path(registry_path::String) + # Standardize to forward slashes internally for uniformity, then normpath formats for OS + clean_path = normpath(strip(replace(registry_path, "\\" => "/"))) + return abspath(clean_path) +end + """ load_registry(registry_path) @@ -914,9 +926,9 @@ Loads the dataset registry from a JSON file, with retries for robustness on Wind # Returns - A dictionary containing the registry data. Returns an empty dictionary if the file doesn't exist or fails to parse. """ -function load_registry(registry_path) +function load_registry(registry_path::String) return lock(REGISTRY_LOCK) do - clean_path = abspath(registry_path) + clean_path = clean_registry_path(registry_path) if isfile(clean_path) for attempt in 1:5 try @@ -937,93 +949,82 @@ function load_registry(registry_path) end """ - extract_metadata(msi_data::MSIData, source_path::String) + extract_metadata(loaded_data::MSIData, full_route::String) -Extracts key metadata from an `MSIData` object for display. - -# Arguments -- `msi_data`: The `MSIData` object. -- `source_path`: The path to the source data file. - -# Returns -- A dictionary containing summary statistics and metadata. +Extracts key metadata from an MSIData object and file path for UI display. +Returns a dictionary with a "summary" key containing parameter-value pairs. """ -function extract_metadata(msi_data::MSIData, source_path::String) - df = msi_data.spectrum_stats_df - if df === nothing - # This can happen if precompute_analytics hasn't been run - # We can still return basic info - return Dict( - "summary" => [ - Dict("parameter" => "File Name", "value" => basename(source_path)), - Dict("parameter" => "Number of Spectra", "value" => length(msi_data.spectra_metadata)), - Dict("parameter" => "Image Dimensions", "value" => "$(msi_data.image_dims[1]) x $(msi_data.image_dims[2])"), - ], - "global_min_mz" => nothing, - "global_max_mz" => nothing - ) +function extract_metadata(loaded_data::MSIData, full_route::String) + stats = [] + + # Basic File Info + push!(stats, Dict("parameter" => "Filename", "value" => basename(full_route))) + push!(stats, Dict("parameter" => "Source Path", "value" => full_route)) + + # Image Dimensions + w, h = loaded_data.image_dims + if w > 0 && h > 0 + push!(stats, Dict("parameter" => "Image Dimensions", "value" => "$w x $h")) + else + push!(stats, Dict("parameter" => "Image Dimensions", "value" => "N/A (Non-imaging)")) end - - summary_stats = [ - Dict("parameter" => "File Name", "value" => basename(source_path)), - Dict("parameter" => "Number of Spectra", "value" => length(msi_data.spectra_metadata)), - Dict("parameter" => "Image Dimensions", "value" => "$(msi_data.image_dims[1]) x $(msi_data.image_dims[2])"), - Dict("parameter" => "Global Min m/z", "value" => @sprintf("%.4f", Threads.atomic_add!(msi_data.global_min_mz, 0.0))), - Dict("parameter" => "Global Max m/z", "value" => @sprintf("%.4f", Threads.atomic_add!(msi_data.global_max_mz, 0.0))), - Dict("parameter" => "Mean TIC", "value" => @sprintf("%.2e", mean(df.TIC))), - Dict("parameter" => "Mean BPI", "value" => @sprintf("%.2e", mean(df.BPI))), - Dict("parameter" => "Mean # Points", "value" => @sprintf("%.1f", mean(df.NumPoints))), - ] - - if hasproperty(df, :Mode) - centroid_count = count(==(MSI_src.CENTROID), df.Mode) - profile_count = count(==(MSI_src.PROFILE), df.Mode) - unknown_count = count(==(MSI_src.UNKNOWN), df.Mode) - - push!(summary_stats, Dict("parameter" => "Centroid Spectra", "value" => string(centroid_count))) - push!(summary_stats, Dict("parameter" => "Profile Spectra", "value" => string(profile_count))) - if unknown_count > 0 - push!(summary_stats, Dict("parameter" => "Unknown Mode Spectra", "value" => string(unknown_count))) + + # Spectrum Counts + num_spectra = length(loaded_data.spectra_metadata) + push!(stats, Dict("parameter" => "Total Spectra", "value" => string(num_spectra))) + + # m/z Range + min_mz, max_mz = get_global_mz_range(loaded_data) + if isfinite(min_mz) && isfinite(max_mz) + push!(stats, Dict("parameter" => "Global m/z Range", "value" => "$(round(min_mz, digits=4)) - $(round(max_mz, digits=4))")) + else + push!(stats, Dict("parameter" => "Global m/z Range", "value" => "Unknown")) + end + + # Instrument Metadata + if loaded_data.instrument_metadata !== nothing + inst = loaded_data.instrument_metadata + push!(stats, Dict("parameter" => "Instrument Model", "value" => inst.instrument_model)) + push!(stats, Dict("parameter" => "Polarity", "value" => titlecase(string(inst.polarity)))) + push!(stats, Dict("parameter" => "Acquisition Mode", "value" => titlecase(string(inst.acquisition_mode)))) + + if inst.resolution !== nothing + push!(stats, Dict("parameter" => "Instrument Resolution", "value" => string(inst.resolution))) end end - - return Dict( - "summary" => summary_stats, - "global_min_mz" => msi_data.global_min_mz, - "global_max_mz" => msi_data.global_max_mz - ) + + return Dict("summary" => stats) end """ - save_registry(registry_path, registry_data) + save_registry(registry_path::String, registry_data::Dict) -Saves the dataset registry to a JSON file, ensuring thread-safe access and robustness on Windows. - -# Arguments -- `registry_path`: Path to the `registry.json` file. -- `registry_data`: The dictionary containing the registry data to save. +Saves the dataset registry to a JSON file atomically to prevent data loss +and sharing violations on Windows. Writes to a `.tmp` file first. """ -function save_registry(registry_path, registry_data) +function save_registry(registry_path::String, registry_data) lock(REGISTRY_LOCK) do - # Ensure the path is clean and valid for the OS - clean_path = abspath(registry_path) + clean_path = clean_registry_path(registry_path) dir = dirname(clean_path) if !isdir(dir) mkpath(dir) end - # On Windows, opening files for writing can occasionally fail if another process - # (like an anti-virus or file indexer) holds a transient lock. - # Retry with a small delay to improve robustness. - for attempt in 1:5 + temp_path = clean_path * ".tmp" + + for attempt in 1:6 try - open(clean_path, "w") do f + open(temp_path, "w") do f JSON.print(f, registry_data, 4) end - return true # Success + + # Replace final file with temp atomically. Force=true maps to ReplaceFile/MoveFileEx + mv(temp_path, clean_path, force=true) + return true catch e - if attempt == 5 + if attempt == 6 @error "Final attempt failed to write to registry.json: $(sprint(showerror, e))" + isfile(temp_path) && rm(temp_path, force=true) return false else @warn "Attempt $attempt to write to registry.json failed (Error: $(sprint(showerror, e))), retrying in 0.2s..." @@ -1047,29 +1048,21 @@ Adds or updates an entry in the dataset registry JSON file. - `metadata`: Optional dictionary of metadata to store. - `is_imzML`: Boolean indicating if the source is an imzML file. """ -function update_registry(registry_path, dataset_name, source_path, metadata=nothing, is_imzML=false) +function update_registry(registry_path::String, dataset_name::String, source_path::String, metadata=nothing, is_imzML=false) lock(REGISTRY_LOCK) do registry = load_registry(registry_path) - # Get existing entry if it exists, otherwise create new one existing_entry = get(registry, dataset_name, Dict{String,Any}()) - - # Start with existing data and update only the basic fields entry = copy(existing_entry) entry["source_path"] = source_path entry["processed_date"] = string(now()) entry["is_imzML"] = is_imzML - # Only update metadata if provided if metadata !== nothing entry["metadata"] = metadata end - # Note: has_mask and mask_path are preserved from existing_entry if they exist - registry[dataset_name] = entry - - # Use the robust save function save_registry(registry_path, registry) end end @@ -1194,24 +1187,14 @@ function process_file_safely(file_path, masses, params, progress_message_ref, ov end end -function update_registry_mask_fields(registry_path, dataset_name, has_mask, mask_path) +function update_registry_mask_fields(registry_path::String, dataset_name::String, has_mask::Bool, mask_path::String) lock(REGISTRY_LOCK) do - registry = if isfile(registry_path) - try - JSON.parsefile(registry_path, dicttype=Dict{String,Any}) - catch e - @error "Failed to parse registry.json while updating mask fields: $e" - Dict{String,Any}() - end - else - Dict{String,Any}() - end + registry = load_registry(registry_path) if haskey(registry, dataset_name) registry[dataset_name]["has_mask"] = has_mask registry[dataset_name]["mask_path"] = mask_path else - # Create new entry if it doesn't exist registry[dataset_name] = Dict{String,Any}( "has_mask" => has_mask, "mask_path" => mask_path, @@ -1221,7 +1204,6 @@ function update_registry_mask_fields(registry_path, dataset_name, has_mask, mask ) end - # Use the robust save function save_registry(registry_path, registry) end end diff --git a/mask.jl b/mask.jl index 45f95e2..4486654 100644 --- a/mask.jl +++ b/mask.jl @@ -12,7 +12,9 @@ using MSI_src using .MSI_src: MSIData, process_image_pipeline, REGISTRY_LOCK # Plot Handling -include("./julia_imzML_visual.jl") +if !@isdefined(increment_image) + include("./julia_imzML_visual.jl") +end # Image Processing Pipeline using ImageBinarization -- 2.43.0 From 70fff16170ca5784651f408bc9e311477300fca4 Mon Sep 17 00:00:00 2001 From: Pixelguy14 Date: Tue, 7 Apr 2026 14:06:02 -0600 Subject: [PATCH 3/6] more secure locks and a script to create a system image of juliaMSI for faster boot times, still experimental but everything else remains pristine --- .gitignore | 4 +- Project.toml | 5 +- README.md | 34 ++++++++- app.jl | 14 ++++ build_sysimage.jl | 35 ++++++++++ julia_imzML_visual.jl | 25 ++++--- precompile_script.jl | 158 ++++++++++++++++++++++++++++++++++++++++++ 7 files changed, 262 insertions(+), 13 deletions(-) create mode 100644 build_sysimage.jl create mode 100644 precompile_script.jl diff --git a/.gitignore b/.gitignore index 0bec6dd..b219dd4 100644 --- a/.gitignore +++ b/.gitignore @@ -7,4 +7,6 @@ log/* R_original_scripts/ test/results/ test/binarization_results/ -Manifest.toml \ No newline at end of file +Manifest.toml +MSI_sysimage.dll +MSI_sysimage.so \ No newline at end of file diff --git a/Project.toml b/Project.toml index 72450c3..48a8ef4 100644 --- a/Project.toml +++ b/Project.toml @@ -35,6 +35,7 @@ Loess = "4345ca2d-374a-55d4-8d30-97f9976e7612" Mmap = "a63ad114-7e13-5084-954f-fe012c677804" NativeFileDialog = "e1fe445b-aa65-4df4-81c1-2041507f0fd4" NaturalSort = "c020b1a1-e9b0-503a-9c33-f039bfc54a85" +PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" @@ -52,9 +53,9 @@ Accessors = "0.1" Base64 = "1.11" CSV = "0.10" CairoMakie = "0.13" +CodecBase = "0.3" ColorSchemes = "3.30" Colors = "0.12" -CodecBase = "0.3" ContinuousWavelets = "1.1" DataFrames = "1.7" Dates = "1.11" @@ -63,7 +64,6 @@ GLMakie = "0.11" Genie = "5.31" GenieFramework = "2.8" HistogramThresholding = "0.3" -Interpolations = "0.15" ImageBinarization = "0.3" ImageComponentAnalysis = "0.2" ImageContrastAdjustment = "0.3" @@ -72,6 +72,7 @@ ImageFiltering = "0.7" ImageMorphology = "0.4" ImageSegmentation = "1.9" Images = "0.26" +Interpolations = "0.15" JLD2 = "0.6" JSON = "0.21" Libz = "1.0" diff --git a/README.md b/README.md index f63ed57..e83f4b9 100644 --- a/README.md +++ b/README.md @@ -22,7 +22,7 @@ https://codeberg.org/LabABI/JuliaMSI Example of a correct route:
~/Downloads/JuliaMSI-main/juliamsi 2. Without entering the Julia environment, launch the project in your terminal with the following command (which works for all operating systems): - ``` + ```bash julia --threads auto --project=. start_MSI_GUI.jl ``` 3. After the script has finished loading, you can open a [page](http://127.0.0.1:1481/) in your browser with the web app running. @@ -34,6 +34,38 @@ Minimum system requirements: 4 core processor, 8 GB RAM
JuliaMSI is a Graphical User Interface for a library of MSI tools in Julia: https://github.com/CINVESTAV-LABI/julia_mzML_imzML +## Build the System Image (One-time) +Alternatively, you can generate the .so file by running the next build script in your directory: + +```bash +julia --project=. build_sysimage.jl +``` +This may take 5–10 minutes as it merges all dependencies and JuliaMSI functions into a single binary. The resulting .so/.dll file will be around **300MB–600MB** because it contains the pre-compiled machine code for your entire environment. +A sysimage built on Linux (.so) will not work on Windows. +You must run the build_sysimage.jl script once on each target operating system. + +Once MSI_sysimage.so is created in your directory, adapt your command like this: +```bash +julia --project=. -e 'using Pkg; Pkg.precompile()' +# This will find MSI_sysimage.so, .dll, or .dylib automatically: +julia --threads auto --project=. --sysimage MSI_sysimage* start_MSI_GUI.jl +``` +The command above loads a pre-compiled version of your environment (note, you have to use either the sysimage command or the normal start command), speeding up the booting delay. Your GUI or scripts using JuliaMSI should start almost instantly and process files significantly faster than before. +Always run the build script **in the same project root** where you intend to run the app. + +In scripts like test/run_tests.jl, you no longer need to manually include("../src/MSI_src.jl"). You can simply treat it like a globally installed package to load the precompiled binary: + +```julia +using MSI_src +# ... rest of your script +``` +### Run with the --sysimage flag +Launch your tests using the same flag to bypass all compilation overhead: + +```bash +julia --threads auto --project=. --sysimage MSI_sysimage.so test/run_tests.jl +``` + ## License JuliaMSI is published under the terms of the MIT License. diff --git a/app.jl b/app.jl index 11058a5..f4bdec4 100644 --- a/app.jl +++ b/app.jl @@ -2422,8 +2422,22 @@ end else push!(newly_created_folders, replace(basename(file_path), r"\.imzML$"i => "")) end + + # --- Periodic Memory Reclamation --- + # Force GC and return memory to the OS after processing each large file + GC.gc(true) + if Sys.islinux() + ccall(:malloc_trim, Cint, (Cint,), 0) + end + current_step += 1 end + + # Final memory cleanup for the batch + GC.gc(true) + if Sys.islinux() + ccall(:malloc_trim, Cint, (Cint,), 0) + end # --- 3. Final Report --- total_time_end = round(time() - total_time_start, digits=3) diff --git a/build_sysimage.jl b/build_sysimage.jl new file mode 100644 index 0000000..f1c8e43 --- /dev/null +++ b/build_sysimage.jl @@ -0,0 +1,35 @@ +# build_sysimage.jl +using PackageCompiler + +println("--- Starting Custom System Image Build ---") +println("This process can take 5-10 minutes.") + +# Define the image path with OS-specific extension +ext = Sys.iswindows() ? ".dll" : Sys.isapple() ? ".dylib" : ".so" +sysimage_path = "MSI_sysimage" * ext + +# List dependencies to include in the sysimage for maximum stability +# We have removed graphical heavy-lifters (CairoMakie, PlotlyBase) to ensure +# the build completes successfully on systems with standard RAM. +packages_to_include = [ + :Genie, + :GenieFramework, + :PlotlyBase, + :Libz, + :DataFrames, + :Serialization, + :Printf, + :JSON +] + +create_sysimage( + packages_to_include, + sysimage_path = sysimage_path, + precompile_execution_file = "precompile_script.jl", + incremental = true, + filter_stdlibs = false +) + +println("--- System Image Build Successful! ---") +println("To start Julia with this image, use:") +println("julia --threads auto --project=. --sysimage $(sysimage_path) start_MSI_GUI.jl") diff --git a/julia_imzML_visual.jl b/julia_imzML_visual.jl index 33acb19..95ae054 100644 --- a/julia_imzML_visual.jl +++ b/julia_imzML_visual.jl @@ -1131,23 +1131,24 @@ function process_file_safely(file_path, masses, params, progress_message_ref, ov # --- Extract metadata --- metadata = extract_metadata(local_msi_data, file_path) - # --- Save Slices --- + # --- Save Slices (Parallelized) --- mkpath(output_dir) - - for (mass_idx, mass) in enumerate(masses) - progress_message_ref = "File $(params.fileIdx)/$(params.nFiles): Saving slice for m/z=$mass" - + + progress_lock = ReentrantLock() + + # Parallelize over the requested masses to fully utilize multi-core CPUs + Threads.@threads for mass in masses slice = slice_dict[mass] text_nmass = replace(string(mass), "." => "_") bitmap_filename = params.triqE ? "TrIQ_$(text_nmass).bmp" : "MSI_$(text_nmass).bmp" colorbar_filename = params.triqE ? "colorbar_TrIQ_$(text_nmass).png" : "colorbar_MSI_$(text_nmass).png" + local sliceQuant, bounds if all(iszero, slice) sliceQuant = zeros(UInt8, size(slice)) bounds = (0.0, 1.0) # Default bounds for empty slice - @warn "No intensity data for m/z = $mass in $(dataset_name)" else - # Get both quantized data AND bounds in one call + # TrIQ and quantization are computationally intensive and now run in parallel if params.triqE sliceQuant, bounds = TrIQ(slice, params.colorL, params.triqP, mask_matrix=mask_matrix_for_triq) else @@ -1156,18 +1157,24 @@ function process_file_safely(file_path, masses, params, progress_message_ref, ov if params.medianF sliceQuant = round.(UInt8, median_filter(sliceQuant)) - # Note: bounds remain the same after median filter end end + # Disk I/O (Saving BMP) save_bitmap(joinpath(output_dir, bitmap_filename), sliceQuant, ViridisPalette) + if !all(iszero, slice) - # Now pass the precomputed bounds to colorbar generation + # CairoMakie colorbar generation is now also parallelized generate_colorbar_image(slice, params.colorL, joinpath(output_dir, colorbar_filename), bounds; use_triq=params.triqE, triq_prob=params.triqP, mask_path=mask_path) end + + lock(progress_lock) do + progress_message_ref = "File $(params.fileIdx)/$(params.nFiles): Saved slice for m/z=$mass" + end end + is_imzML = local_msi_data.source isa ImzMLSource update_registry(params.registry, dataset_name, file_path, metadata, is_imzML) return (true, "") diff --git a/precompile_script.jl b/precompile_script.jl new file mode 100644 index 0000000..050df9a --- /dev/null +++ b/precompile_script.jl @@ -0,0 +1,158 @@ +# precompile_script.jl +# This script exercises the core MSI processing kernels (imzML and mzML) +# to ensure they are precompiled into the custom system image. + +using MSI_src +using Mmap +using Base64 +using PlotlyBase + +# --- DYNAMIC WARM-UP WORKLOAD --- # + +function create_warmup_datasets(dir, T_mz::Type, T_int::Type) + # 1. Create minimal imzML/ibd + imzml_path = joinpath(dir, "warmup_$(T_mz)_$(T_int).imzML") + ibd_path = joinpath(dir, "warmup_$(T_mz)_$(T_int).ibd") + + n_points = 10 + n_pixels = 1 + mz_bytes = n_points * sizeof(T_mz) + int_bytes = n_points * sizeof(T_int) + + # Write some dummy binary data + open(ibd_path, "w") do f + write(f, rand(T_mz, n_points)) + write(f, rand(T_int, n_points)) + end + + # Minimal but VALID imzML structure that passes axes_config_img + imzml_content = """ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +""" + write(imzml_path, imzml_content) + + # 2. Create minimal mzML (Base64) + mzml_path = joinpath(dir, "warmup_$(T_mz)_$(T_int).mzML") + b64_mz = base64encode(rand(T_mz, n_points)) + b64_int = base64encode(rand(T_int, n_points)) + + mzml_content = """ + + + + + + + + + $b64_mz + + + + + $b64_int + + + + + + PLACEHOLDER + ID_OFFSET + +""" + spec_start = findfirst(" string(spec_start)) + mzml_content = replace(mzml_content, "ID_OFFSET" => string(index_start)) + + write(mzml_path, mzml_content) + + return imzml_path, mzml_path +end + +println("Starting robust precompilation warmup (imzML + mzML)...") + +try + mktempdir() do tmp_dir + for (T_mz, T_int) in [(Float32, Float32), (Float64, Float32)] + println("Exercising pathways for $T_mz and $T_int...") + imzml_path, mzml_path = create_warmup_datasets(tmp_dir, T_mz, T_int) + + try + # Exercise imzML pathway + imzml_data = MSI_src.load_imzml_lazy(imzml_path; use_mmap=true) + MSI_src.get_multiple_mz_slices(imzml_data, [10.0, 20.0], 0.1) + + # Exercise mzML pathway + mzml_data = MSI_src.load_mzml_lazy(mzml_path) + MSI_src.GetSpectrum(mzml_data, 1) + catch e + @warn "Warmup failed for $T_mz/$T_int: $e" + end + end + + # --- COMMON KERNELS --- + println("Precompiling image processing kernels...") + dummy_slice = rand(Float32, 4, 4) + MSI_src.TrIQ(dummy_slice, 256, 1.0) + MSI_src.save_bitmap(joinpath(tmp_dir, "d.bmp"), zeros(UInt8, 4, 4), MSI_src.ViridisPalette) + + # Plotly Warmup + p = PlotlyBase.Plot(PlotlyBase.scatter(x=1:2, y=[1,2])) + end +catch e + @warn "Global Warmup failed: $e" +catch + # Silent fail for interrupt +end +println("Precompilation workload completed.") -- 2.43.0 From 8c37ea64aa872787152686076f1f7c789be6e2a3 Mon Sep 17 00:00:00 2001 From: Pixelguy14 Date: Thu, 9 Apr 2026 12:30:43 -0600 Subject: [PATCH 4/6] fixed disk and memory bloat when loading several files, created zero-copy mmaps for msidata type objects for easier cache storage and less memory overhead, improved load times for files in general across the whole JuliaMSI environment --- Project.toml | 1 + app.jl | 197 ++++++----- src/Common.jl | 31 +- src/FusedPipeline.jl | 93 +++++ src/MSIData.jl | 532 ++++++++++++++++++----------- src/MSI_src.jl | 16 + src/Precalculations.jl | 42 ++- src/Preprocessing.jl | 69 ++-- src/ResourcePool.jl | 95 ++++++ src/StreamingKernels.jl | 317 +++++++++++++++++ src/StreamingPipeline.jl | 402 ++++++++++++++++++++++ src/imzML.jl | 404 ++++++++++++---------- src/mzML.jl | 80 +++-- start_MSI_GUI.jl | 46 +++ test/new_benchmark_mmap.jl | 48 +++ test/run_precalculation_example.jl | 3 +- test/run_preprocessing.jl | 8 +- test/run_tests.jl | 12 +- test/test_streaming_pipeline.jl | 124 +++++++ 19 files changed, 1986 insertions(+), 534 deletions(-) create mode 100644 src/FusedPipeline.jl create mode 100644 src/ResourcePool.jl create mode 100644 src/StreamingKernels.jl create mode 100644 src/StreamingPipeline.jl create mode 100644 test/new_benchmark_mmap.jl create mode 100644 test/test_streaming_pipeline.jl diff --git a/Project.toml b/Project.toml index 48a8ef4..e171608 100644 --- a/Project.toml +++ b/Project.toml @@ -43,6 +43,7 @@ ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" SavitzkyGolay = "c4bf5708-b6a6-4fbe-bcd0-6850ed671584" Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StipplePlotly = "ec984513-233d-481d-95b0-a3b58b97af2b" diff --git a/app.jl b/app.jl index f4bdec4..534748c 100644 --- a/app.jl +++ b/app.jl @@ -29,6 +29,9 @@ if !@isdefined(increment_image) include("./julia_imzML_visual.jl") end +const global_msi_data = Ref{Union{MSIData, Nothing}}(nothing) + + # --- Memory Validation Logging --- if get(ENV, "GENIE_ENV", "dev") != "prod" function get_rss_mb() @@ -60,7 +63,7 @@ if get(ENV, "GENIE_ENV", "dev") != "prod" println("--- MEMORY LOG [$(context)] ---") println(" Timestamp: $(now())") println(" Process RSS: $(rss_mb) MB") - println(" msi_data size: $(msi_data_size_mb) MB") + println(" global_msi_data[] size: $(msi_data_size_mb) MB") println(" Cumulative GC time: $(gc_time_s) s") println("--------------------------") end @@ -193,6 +196,15 @@ macro ui_log(message, level="INFO", log_entries) end =# +# --- CRITICAL: Disable Stipple's session-to-disk persistence --- +# Stipple's ModelStorage registers on(field) handlers that serialize the ENTIRE +# ReactiveModel to disk via GenieSessionFileSession on every UI state change. +# With 253 reactive variables including multi-MB Plotly traces, this generates +# gigabytes of orphaned session files in /tmp/jl_XXXXXX, exhausting disk space. +# For a single-user desktop application, session persistence is unnecessary. +Stipple.enable_model_storage(false) +Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage is disabled") during layout render + # Reactive code to make the UI interactive @app begin # == Notification & Logs == @@ -476,7 +488,7 @@ end # == DATA MANAGEMENT VARIABLES == # Centralized MSIData object - @out msi_data::Union{MSIData, Nothing} = nothing + # global_msi_data[] is now global to avoid Genie Session memory leaks # Image file management @out text_nmass="" # For specific mass charge image creation @@ -661,7 +673,7 @@ end try # 1. Clear large data objects explicitly - msi_data = nothing + global_msi_data[] = nothing feature_matrix_result = nothing bin_info_result = nothing @@ -717,30 +729,42 @@ end sure the file can be processed by later steps like mainProcess =# @onbutton btnSearch begin - is_processing = true - push!(__model__) + # 0. Robustness Guard: Prevent double-trigger during processing + if is_processing + println("DEBUG: btnSearch ignored because another process is already running.") + return + end + btnSearch = false # Manual reset of the trigger + + # 1. Grab the file path from the main task picked_route = pick_file(; filterlist="imzML,imzml,mzML,mzml") if isnothing(picked_route) || isempty(picked_route) - is_processing = false return end - # --- Close previous dataset if one is open --- - if msi_data !== nothing - println("DEBUG: Closing previously loaded dataset before opening new one: $(basename(full_route))") - close(msi_data) - msi_data = nothing - GC.gc() - if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) - end - end - + # 2. Update reactive state synchronously + is_processing = true msg = "Opening file: $(basename(picked_route))..." + SpectraEnabled = false + btnMetadataDisable = true + push!(__model__) - try - dataset_name = replace(basename(picked_route), r"(\.(imzML|imzml|mzML|mzml))$"i => "") + # 3. Spawn background computational thread + Threads.@spawn begin + try + # --- Close previous dataset if one is open --- + if global_msi_data[] !== nothing + println("DEBUG: Closing previously loaded dataset before opening new one...") + close(global_msi_data[]) + global_msi_data[] = nothing + GC.gc() + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) + end + end + + dataset_name = replace(basename(picked_route), r"(\.(imzML|imzml|mzML|mzml))$"i => "") registry = load_registry(registry_path) existing_entry = get(registry, dataset_name, nothing) @@ -757,8 +781,8 @@ end dims = parse.(Int, split(dims_str, " x ")) imgWidth, imgHeight = dims[1], dims[2] - msi_data = nothing # Ensure data is not held in memory - log_memory_usage("Fast Load (msi_data cleared)", msi_data) + global_msi_data[] = nothing # Ensure data is not held in memory + log_memory_usage("Fast Load (global_msi_data[] cleared)", global_msi_data[]) btnMetadataDisable = false SpectraEnabled = true selected_folder_main = dataset_name @@ -1000,10 +1024,10 @@ end image_available_folders = deepcopy(img_folders) selected_folder_main = dataset_name - msi_data = loaded_data + global_msi_data[] = loaded_data # Determine plot mode from loaded data - df = msi_data.spectrum_stats_df + df = global_msi_data[].spectrum_stats_df if df !== nothing && "Mode" in names(df) profile_count = count(==(MSI_src.PROFILE), df.Mode) total_count = length(df.Mode) @@ -1013,26 +1037,29 @@ end last_plot_mode = "lines" # Default end - log_memory_usage("Full Load", msi_data) + log_memory_usage("Full Load", global_msi_data[]) eTime = round(time() - sTime, digits=3) msg = "Active file loaded in $(eTime) seconds. Dataset '$(dataset_name)' is ready for analysis." SpectraEnabled = true - catch e - msi_data = nothing - msg = "Error loading active file: $e" - warning_msg = true - SpectraEnabled = false - btnMetadataDisable = true - @error "File loading failed" exception=(e, catch_backtrace()) - finally - GC.gc() # Trigger garbage collection - if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS + catch e + global_msi_data[] = nothing + msg = "Error loading active file: $e" + warning_msg = true + SpectraEnabled = false + btnMetadataDisable = true + @error "File loading failed" exception=(e, catch_backtrace()) + push!(__model__) # Force sending error back to UI immediately + finally + GC.gc() # Trigger garbage collection + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS + end + is_processing = false + push!(__model__) # Clear spinner loop end - is_processing = false end end @@ -1232,6 +1259,11 @@ end This reactive handler job is to run the full preprocessing pipeline on the selected dataset. =# @onbutton run_full_pipeline begin + if is_processing + println("DEBUG: run_full_pipeline ignored because another process is already running.") + return + end + run_full_pipeline = false # Manual reset is_processing = true push!(__model__) overall_progress = 0.0 @@ -1259,9 +1291,9 @@ end end target_path = entry["source_path"] - # Ensure msi_data is for the currently selected file and load if needed + # Ensure global_msi_data[] is for the currently selected file and load if needed # NOTE: For the pipeline, we will open a DEDICATED instance to avoid race conditions - # with the global msi_data used for plotting/interactive exploration. + # with the global global_msi_data[] used for plotting/interactive exploration. println("DEBUG: Opening isolated MSIData instance for pipeline stability...") pipeline_msi_data = OpenMSIData(target_path) @@ -1686,9 +1718,9 @@ end # Determine plot mode for this specific spectrum spectrum_mode_for_plot = "lines" # Default to lines - if msi_data.spectrum_stats_df !== nothing && "Mode" in names(msi_data.spectrum_stats_df) - if selected_spectrum_id_for_plot > 0 && selected_spectrum_id_for_plot <= length(msi_data.spectrum_stats_df.Mode) - mode = msi_data.spectrum_stats_df.Mode[selected_spectrum_id_for_plot] + if global_msi_data[].spectrum_stats_df !== nothing && "Mode" in names(global_msi_data[].spectrum_stats_df) + if selected_spectrum_id_for_plot > 0 && selected_spectrum_id_for_plot <= length(global_msi_data[].spectrum_stats_df.Mode) + mode = global_msi_data[].spectrum_stats_df.Mode[selected_spectrum_id_for_plot] if mode == MSI_src.CENTROID spectrum_mode_for_plot = "stem" end @@ -1939,7 +1971,7 @@ end @onbutton recalculate_suggestions_btn begin is_processing = true push!(__model__) - if msi_data === nothing + if global_msi_data[] === nothing msg = "Please load a file first." warning_msg = true return @@ -1953,7 +1985,7 @@ end for p in reference_peaks_list if tryparse(Float64, string(p["mz"])) !== nothing ) - recommended_params = main_precalculation(msi_data, reference_peaks=ref_peaks) + recommended_params = main_precalculation(global_msi_data[], reference_peaks=ref_peaks) for (step_name, params) in recommended_params @@ -2348,6 +2380,11 @@ end end @onbutton mainProcess @time begin + if is_processing + println("DEBUG: mainProcess ignored because another process is already running.") + return + end + mainProcess = false # Manual reset # --- UI State Update --- overall_progress = 0.0 progress_message = "Preparing batch process..." @@ -2572,21 +2609,21 @@ end return end - if msi_data === nothing || full_route != target_path - if msi_data !== nothing - close(msi_data) + if global_msi_data[] === nothing || full_route != target_path + if global_msi_data[] !== nothing + close(global_msi_data[]) end msg = "Reloading $(basename(target_path)) for analysis..." full_route = target_path - msi_data = OpenMSIData(target_path) + global_msi_data[] = OpenMSIData(target_path) if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing raw_min = entry["metadata"]["global_min_mz"] raw_max = entry["metadata"]["global_max_mz"] min_val = isa(raw_min, Dict) ? get(raw_min, "value", raw_min) : raw_min max_val = isa(raw_max, Dict) ? get(raw_max, "value", raw_max) : raw_max - set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val)) + set_global_mz_range!(global_msi_data[], convert(Float64, min_val), convert(Float64, max_val)) else - precompute_analytics(msi_data) + precompute_analytics(global_msi_data[]) end end @@ -2599,7 +2636,7 @@ end end end - plotdata, plotlayout, xSpectraMz, ySpectraMz = meanSpectrumPlot(msi_data, selected_folder_main, mask_path=mask_path_for_plot) + plotdata, plotlayout, xSpectraMz, ySpectraMz = meanSpectrumPlot(global_msi_data[], selected_folder_main, mask_path=mask_path_for_plot) plotdata_before = plotdata plotlayout_before = plotlayout last_plot_type = "mean" @@ -2607,7 +2644,7 @@ end fTime = time() eTime = round(fTime - sTime, digits=3) msg = "Plot loaded in $(eTime) seconds" - log_memory_usage("Mean Plot Generated", msi_data) + log_memory_usage("Mean Plot Generated", global_msi_data[]) catch e msg = "Could not generate mean spectrum plot: $e" warning_msg = true @@ -2661,21 +2698,21 @@ end return end - if msi_data === nothing || full_route != target_path - if msi_data !== nothing - close(msi_data) + if global_msi_data[] === nothing || full_route != target_path + if global_msi_data[] !== nothing + close(global_msi_data[]) end msg = "Reloading $(basename(target_path)) for analysis..." full_route = target_path - msi_data = OpenMSIData(target_path) + global_msi_data[] = OpenMSIData(target_path) if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing raw_min = entry["metadata"]["global_min_mz"] raw_max = entry["metadata"]["global_max_mz"] min_val = isa(raw_min, Dict) ? get(raw_min, "value", raw_min) : raw_min max_val = isa(raw_max, Dict) ? get(raw_max, "value", raw_max) : raw_max - set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val)) + set_global_mz_range!(global_msi_data[], convert(Float64, min_val), convert(Float64, max_val)) else - precompute_analytics(msi_data) + precompute_analytics(global_msi_data[]) end end local mask_path_for_plot::Union{String, Nothing} = nothing @@ -2687,7 +2724,7 @@ end end end - plotdata, plotlayout, xSpectraMz, ySpectraMz = sumSpectrumPlot(msi_data, selected_folder_main, mask_path=mask_path_for_plot) + plotdata, plotlayout, xSpectraMz, ySpectraMz = sumSpectrumPlot(global_msi_data[], selected_folder_main, mask_path=mask_path_for_plot) plotdata_before = plotdata plotlayout_before = plotlayout last_plot_type = "sum" @@ -2695,7 +2732,7 @@ end fTime = time() eTime = round(fTime - sTime, digits=3) msg = "Total plot loaded in $(eTime) seconds" - log_memory_usage("Sum Plot Generated", msi_data) + log_memory_usage("Sum Plot Generated", global_msi_data[]) catch e msg = "Could not generate total spectrum plot: $e" warning_msg = true @@ -2752,21 +2789,21 @@ end return end - if msi_data === nothing || full_route != target_path - if msi_data !== nothing - close(msi_data) + if global_msi_data[] === nothing || full_route != target_path + if global_msi_data[] !== nothing + close(global_msi_data[]) end msg = "Reloading $(basename(target_path)) for analysis..." full_route = target_path - msi_data = OpenMSIData(target_path) + global_msi_data[] = OpenMSIData(target_path) if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing raw_min = entry["metadata"]["global_min_mz"] raw_max = entry["metadata"]["global_max_mz"] min_val = isa(raw_min, Dict) ? get(raw_min, "value", raw_min) : raw_min max_val = isa(raw_max, Dict) ? get(raw_max, "value", raw_max) : raw_max - set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val)) + set_global_mz_range!(global_msi_data[], convert(Float64, min_val), convert(Float64, max_val)) else - precompute_analytics(msi_data) + precompute_analytics(global_msi_data[]) end end @@ -2781,7 +2818,7 @@ end # Convert to positive coordinates for processing y_positive = yCoord < 0 ? abs(yCoord) : yCoord - plotdata, plotlayout, xSpectraMz, ySpectraMz, spectrum_id = xySpectrumPlot(msi_data, xCoord, y_positive, imgWidth, imgHeight, selected_folder_main, mask_path=mask_path_for_plot) + plotdata, plotlayout, xSpectraMz, ySpectraMz, spectrum_id = xySpectrumPlot(global_msi_data[], xCoord, y_positive, imgWidth, imgHeight, selected_folder_main, mask_path=mask_path_for_plot) plotdata_before = plotdata plotlayout_before = plotlayout last_plot_type = "single" @@ -2822,7 +2859,7 @@ end fTime = time() eTime = round(fTime - sTime, digits=3) msg = "Plot loaded in $(eTime) seconds" - log_memory_usage("XY Plot Generated", msi_data) + log_memory_usage("XY Plot Generated", global_msi_data[]) catch e msg = "Could not retrieve spectrum: $e" warning_msg = true @@ -2867,21 +2904,21 @@ end return end - if msi_data === nothing || full_route != target_path - if msi_data !== nothing - close(msi_data) + if global_msi_data[] === nothing || full_route != target_path + if global_msi_data[] !== nothing + close(global_msi_data[]) end msg = "Reloading $(basename(target_path)) for analysis..." full_route = target_path - msi_data = OpenMSIData(target_path) + global_msi_data[] = OpenMSIData(target_path) if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing raw_min = entry["metadata"]["global_min_mz"] raw_max = entry["metadata"]["global_max_mz"] min_val = isa(raw_min, Dict) ? get(raw_min, "value", raw_min) : raw_min max_val = isa(raw_max, Dict) ? get(raw_max, "value", raw_max) : raw_max - set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val)) + set_global_mz_range!(global_msi_data[], convert(Float64, min_val), convert(Float64, max_val)) else - precompute_analytics(msi_data) + precompute_analytics(global_msi_data[]) end end @@ -2895,7 +2932,7 @@ end end # Call the new nSpectrumPlot function - plotdata, plotlayout, xSpectraMz, ySpectraMz, spectrum_id = nSpectrumPlot(msi_data, idSpectrum, selected_folder_main, mask_path=mask_path_for_plot) + plotdata, plotlayout, xSpectraMz, ySpectraMz, spectrum_id = nSpectrumPlot(global_msi_data[], idSpectrum, selected_folder_main, mask_path=mask_path_for_plot) plotdata_before = plotdata plotlayout_before = plotlayout last_plot_type = "single" @@ -2905,7 +2942,7 @@ end fTime = time() eTime = round(fTime - sTime, digits=3) msg = "Plot loaded in $(eTime) seconds" - log_memory_usage("nSpectrum Plot Generated", msi_data) + log_memory_usage("nSpectrum Plot Generated", global_msi_data[]) catch e msg = "Could not retrieve spectrum: $e" warning_msg = true @@ -3186,7 +3223,7 @@ end # This handler will now correctly load the first image from the newly selected folder. @onchange selected_folder_main begin - # The msi_data object lifecycle is managed by the btnSearch handler. + # The global_msi_data[] object lifecycle is managed by the btnSearch handler. # This handler is now only for updating the UI images when the folder changes. if !isempty(selected_folder_main) @@ -3418,7 +3455,7 @@ end fTime = time() eTime = round(fTime - sTime, digits=3) msg = "Plot loaded in $(eTime) seconds" - log_memory_usage("Mean Plot Generated", msi_data) + log_memory_usage("Mean Plot Generated", global_msi_data[]) catch e msg = "Failed to load and process image: $e" warning_msg = true @@ -3476,7 +3513,7 @@ end fTime = time() eTime = round(fTime - sTime, digits=3) msg = "Plot loaded in $(eTime) seconds" - log_memory_usage("Mean Plot Generated", msi_data) + log_memory_usage("Mean Plot Generated", global_msi_data[]) catch e msg = "Failed to load and process image: $e" warning_msg = true @@ -3577,7 +3614,7 @@ end # To include a visualization in the spectrum plot indicating where is the selected mass @onchange Nmass begin if !isempty(xSpectraMz) - df = msi_data.spectrum_stats_df + df = global_msi_data[].spectrum_stats_df plot_as_lines = false # Default to stem if df !== nothing && hasproperty(df, :Mode) && !isempty(df.Mode) profile_count = count(==(MSI_src.PROFILE), df.Mode) @@ -3836,7 +3873,7 @@ end eTime=round(fTime-sTime,digits=3) is_initializing = false # Hide loading screen when initialization is complete (current code is hidden due to incompatibility) msg = "The app took $(eTime) seconds to get ready." - log_memory_usage("App Ready", msi_data) + log_memory_usage("App Ready", global_msi_data[]) end # is_processing = false GC.gc() # Trigger garbage collection diff --git a/src/Common.jl b/src/Common.jl index 668e816..de1139c 100644 --- a/src/Common.jl +++ b/src/Common.jl @@ -1,5 +1,34 @@ -# src/Common.jl - Updated with BloomFilter using Base.Threads +using Mmap + +# POSIX madvise constants +const MADV_NORMAL = 0 +const MADV_RANDOM = 1 +const MADV_SEQUENTIAL = 2 +const MADV_WILLNEED = 3 +const MADV_DONTNEED = 4 + +""" + posix_madvise(buffer::AbstractArray, advice::Integer) + +A safe wrapper for the OS `madvise` system call. Signals the kernel about the +access pattern for a memory-mapped region. Currently supports Linux/Unix systems. +""" +function posix_madvise(buffer::AbstractArray, advice::Integer) + @static if Sys.isunix() + try + ptr = pointer(buffer) + len = sizeof(buffer) + # ccall(:madvise, return_type, (arg_types...), args...) + ret = ccall(:madvise, Int32, (Ptr{Cvoid}, Csize_t, Int32), ptr, len, Int32(advice)) + return ret == 0 + catch + return false + end + else + return false # Not supported on this OS + end +end # --- Buffer Pooling --- """ diff --git a/src/FusedPipeline.jl b/src/FusedPipeline.jl new file mode 100644 index 0000000..47f835a --- /dev/null +++ b/src/FusedPipeline.jl @@ -0,0 +1,93 @@ +# src/FusedPipeline.jl +# This file defines a high-performance in-place preprocessing pipeline +# that minimizes allocations by using reused buffers from ResourcePool. + +using .MSI_src # Ensure it can see the module's contents if needed + +""" + SpectralPipeline + +Holds a sequence of preprocessing steps and the necessary buffers to execute them in-place. +""" +struct SpectralPipeline + steps::Vector{AbstractPreprocessingStep} +end + +""" + apply_pipeline!(mz::Vector{Float64}, intensity::Vector{Float64}, pipeline::SpectralPipeline; data::MSIData) + +Applies all steps in the pipeline to the spectrum arrays in-place. +Uses internal buffers from data.resource_pool where needed. +""" +function apply_pipeline!(mz::Vector{Float64}, intensity::Vector{Float64}, pipeline::SpectralPipeline, data::MSIData) + # Process each step in sequence + for step in pipeline.steps + apply_step!(mz, intensity, step, data) + end + return mz, intensity +end + +# --- Basic in-place implementations of core preprocessing steps --- + +function apply_step!(mz, int, step::Normalization, data) + if step.method === :tic + s = sum(int) + if s > 0 + int ./= s + end + elseif step.method === :median + m = median(int) + if m > 0 + int ./= m + end + end + return int +end + +function apply_step!(mz, int, step::Smoothing, data) + if step.method === :savitzky_golay + # SavitzkyGolay.savitzky_golay currently allocates, but we can't easily fix that here + # without refactoring the library. However, we can use a pooled vector for its output + # then copy back to intensity. + # [Wait: For now we'll call smoothed_y = smooth_spectrum_core(int, ...)] + # We'll use a resource from the pool to avoid fresh allocation + temp_buf = acquire(data.resource_pool) + resize!(temp_buf, length(int)) + + # Call existing core which returns a new vector, unfortunately + # But we'll copy it back to 'int' to maintain in-place pipeline + smoothed = smooth_spectrum_core(int; method=step.method, window=step.window, order=step.order) + copyto!(int, smoothed) + + release!(data.resource_pool, temp_buf) + end + return int +end + +function apply_step!(mz, int, step::BaselineCorrection, data) + if step.method === :snip + # SNIP is easy to make in-place! + iterations = (step.iterations === nothing) ? 100 : step.iterations + snip_baseline_inplace!(int, iterations) + end + return int +end + +""" + snip_baseline_inplace!(y, iterations) + +In-place implementation of the Sensitive Nonlinear Iterative Peak clipping algorithm. +""" +function snip_baseline_inplace!(y::Vector{Float64}, iterations::Int) + n = length(y) + n < 3 && return y + + # We still need one temporary buffer for the SNIP iteration to read from the previous state + # Actually, we can just return the baseline and subtract it, but to BE in-place, + # we need a temporary to hold the baseline during calculation. + + # For now, we'll use the existing _snip_baseline_impl and subtract + baseline = _snip_baseline_impl(y, iterations=iterations) + y .-= baseline + return y +end diff --git a/src/MSIData.jl b/src/MSIData.jl index 96a01ca..b40a7ef 100644 --- a/src/MSIData.jl +++ b/src/MSIData.jl @@ -6,7 +6,7 @@ including caching and iteration logic, for handling large mzML and imzML dataset efficiently. """ -using Base64, Libz, Serialization, Printf, DataFrames, Base.Threads, StatsBase +using Base64, Libz, Serialization, Printf, DataFrames, Base.Threads, StatsBase, Mmap const FILE_HANDLE_LOCK = ReentrantLock() @@ -117,9 +117,11 @@ A data source for `.imzML` files, holding a handle to the binary `.ibd` file and the expected format for m/z and intensity arrays. """ struct ImzMLSource <: MSDataSource - ibd_handle::Union{IO, ThreadSafeFileHandle} + ibd_handles::Vector{IO} # HandlePool: One handle per thread mz_format::Type intensity_format::Type + mmap_data::Union{Vector{UInt8}, Nothing} + is_any_compressed::Bool # Cached for zero-allocation dispatch end """ @@ -129,9 +131,35 @@ A data source for `.mzML` files, holding a handle to the `.mzML` file itself (which contains the binary data encoded in Base64) and the expected data formats. """ struct MzMLSource <: MSDataSource - file_handle::Union{IO, ThreadSafeFileHandle} + file_handles::Vector{IO} # HandlePool: One handle per thread mz_format::Type intensity_format::Type + mmap_data::Union{Vector{UInt8}, Nothing} +end + +# --- HandlePool Helpers --- # +""" + get_handle(source::ImzMLSource) -> IO + get_handle(source::MzMLSource) -> IO + +Retrieves a thread-local file handle from the source's pool. +""" +function get_handle(source::ImzMLSource) + tid = Threads.threadid() + if tid <= length(source.ibd_handles) + return source.ibd_handles[tid] + else + return source.ibd_handles[1] + end +end + +function get_handle(source::MzMLSource) + tid = Threads.threadid() + if tid <= length(source.file_handles) + return source.file_handles[tid] + else + return source.file_handles[1] + end end """ @@ -157,6 +185,10 @@ struct SpectrumAsset # For mzML, axis_type is needed to distinguish mz from intensity. # For imzML, this can be ignored as the order is fixed. axis_type::Symbol + + # Pre-computed analytics + min_val::Float64 + max_val::Float64 end """ @@ -241,6 +273,24 @@ struct SpectrumMetadata int_asset::SpectrumAsset end +""" + SpectrumMetadataBinary + +A fixed-size version of SpectrumMetadata for high-speed binary serialization. +Used for the metadata cache (.cache files). +""" +struct SpectrumMetadataBinary + x::Int32 + y::Int32 + mode::Int8 + mz_offset::Int64 + mz_encoded_len::Int32 + int_offset::Int64 + int_encoded_len::Int32 + min_mz::Float32 # Persistent analytics + max_mz::Float32 # Persistent analytics +end + """ MSIData @@ -275,7 +325,7 @@ mutable struct MSIData cache_lock::ReentrantLock # Buffer Pool for binary data operations - buffer_pool::SimpleBufferPool + buffer_pool::SimpleBufferPool # Existing UInt8 pool (mostly for mzML base64) # Pre-computed analytics/metadata - use Base.Threads.Atomic for compatibility global_min_mz::Base.Threads.Atomic{Float64} @@ -287,17 +337,26 @@ mutable struct MSIData function MSIData(source, metadata, instrument_meta, dims, coordinate_map, cache_size) obj = new(source, metadata, instrument_meta, dims, coordinate_map, - Dict(), [], cache_size, ReentrantLock(), + Dict(), [], min(10, cache_size), ReentrantLock(), SimpleBufferPool(), - Base.Threads.Atomic{Float64}(Inf), Base.Threads.Atomic{Float64}(-Inf), + Base.Threads.Atomic{Float64}(0.0), Base.Threads.Atomic{Float64}(0.0), nothing, nothing, AtomicFlag(), nothing) - # Ensure file handles are closed when the object is garbage collected + # Initialize mz bounds cleanly instead of Inf + Base.Threads.atomic_xchg!(obj.global_min_mz, 1e9) + Base.Threads.atomic_xchg!(obj.global_max_mz, -1e9) + + + # Ensure all file handles in the pool are closed when the object is garbage collected finalizer(obj) do o - if o.source isa ImzMLSource && isopen(o.source.ibd_handle) - close(o.source.ibd_handle) - elseif o.source isa MzMLSource && isopen(o.source.file_handle) - close(o.source.file_handle) + if o.source isa ImzMLSource + for h in o.source.ibd_handles + isopen(h) && close(h) + end + elseif o.source isa MzMLSource + for h in o.source.file_handles + isopen(h) && close(h) + end end end return obj @@ -394,9 +453,7 @@ Gets the spectrum statistics for the MSIData object. - `stats_df::DataFrame`: The statistics DataFrame. """ function get_spectrum_stats(data::MSIData) - lock(data.cache_lock) do - return data.spectrum_stats_df - end + return data.spectrum_stats_df end """ @@ -413,10 +470,14 @@ It is good practice to call this method when you are finished with an `MSIData` - `nothing` """ function Base.close(data::MSIData) - if data.source isa ImzMLSource && isopen(data.source.ibd_handle) - close(data.source.ibd_handle) - elseif data.source isa MzMLSource && isopen(data.source.file_handle) - close(data.source.file_handle) + if data.source isa ImzMLSource + for handle in data.source.ibd_handles + isopen(handle) && close(handle) + end + elseif data.source isa MzMLSource + for handle in data.source.file_handles + isopen(handle) && close(handle) + end end # Clear cache @@ -552,38 +613,36 @@ by this function and is assumed to be handled by the caller if necessary. - A `Vector` of the appropriate type containing the decoded data. """ function read_binary_vector(data::MSIData, io::IO, asset::SpectrumAsset) - if asset.offset < 0 || asset.offset >= filesize(io) - throw(FileFormatError("Invalid asset offset: $(asset.offset) for file size $(filesize(io))")) + if asset.offset < 0 + throw(FileFormatError("Invalid asset offset: $(asset.offset)")) end - seek(io, asset.offset) - raw_b64 = read(io, asset.encoded_length) - - # Use String directly to avoid intermediate allocations - b64_string = String(raw_b64) + # Use mmap view if available for Base64 (mzML) + b64_string = (data.source isa MzMLSource && data.source.mmap_data !== nothing) ? + String(view(data.source.mmap_data, (asset.offset + 1):(asset.offset + asset.encoded_length))) : + String(read(seek(io, asset.offset), asset.encoded_length)) local decoded_bytes::Vector{UInt8} if asset.is_compressed - # Direct Base64 decode to temporary, then decompress temp_decoded = Base64.base64decode(b64_string) decoded_bytes = Libz.inflate(temp_decoded) else - # Direct Base64 decode decoded_bytes = Base64.base64decode(b64_string) end - # Calculate number of elements - n_elements = length(decoded_bytes) ÷ sizeof(asset.format) + alignment = sizeof(asset.format) + n_elements = length(decoded_bytes) ÷ alignment - if n_elements * sizeof(asset.format) != length(decoded_bytes) + if n_elements * alignment != length(decoded_bytes) throw(FileFormatError("Size of decoded byte array is not a multiple of the element size.")) end - # Reinterpret the byte array as an array of the target type. This does not copy. reinterpreted_array = reinterpret(asset.format, decoded_bytes) - # Allocate the final output array and convert byte order while copying. + # Optimization: Use a temporary array for byte order conversion. + # We could use the ResourcePool here if we wanted to return a Float64 vector, + # but currently we return the native format. out_array = [ltoh(x) for x in reinterpreted_array] return out_array @@ -625,35 +684,51 @@ and converting it from little-endian to the host's native byte order. # Returns - A tuple `(mz, intensity)` containing the two requested data arrays. """ -function read_spectrum_from_disk(source::ImzMLSource, meta::SpectrumMetadata) - # For imzML, the binary data is raw, not base64 encoded. - # The `encoded_length` field in this case holds the number of points. +@inline function read_spectrum_from_disk(source::ImzMLSource, meta::SpectrumMetadata) + # 1. Use Mmap logic if available (implemented in previous step) + if source.mmap_data !== nothing + mz_offset = meta.mz_asset.offset + mz_byte_len = sizeof(source.mz_format) * meta.mz_asset.encoded_length + int_offset = meta.int_asset.offset + int_byte_len = sizeof(source.intensity_format) * meta.int_asset.encoded_length + + mz_view_raw = view(source.mmap_data, (mz_offset + 1):(mz_offset + mz_byte_len)) + int_view_raw = view(source.mmap_data, (int_offset + 1):(int_offset + int_byte_len)) + + mz_reinterpreted = reinterpret(source.mz_format, mz_view_raw) + int_reinterpreted = reinterpret(source.intensity_format, int_view_raw) + + # TRUE Zero-copy logic: avoid allocations if host matches file endianness (LE for .ibd) + if Base.ENDIAN_BOM == 0x04030201 # Little Endian Host (Common for Linux/X86) + mz = mz_reinterpreted + intensity = int_reinterpreted + else + # On Big Endian hosts, we MUST allocate and byte-swap + mz = ltoh.(mz_reinterpreted) + intensity = ltoh.(int_reinterpreted) + end + + validate_spectrum_data(mz, intensity, meta.id) + return mz, intensity + end + + # 2. Use HandlePool logic if Mmap is not available + handle = get_handle(source) + mz = Array{source.mz_format}(undef, meta.mz_asset.encoded_length) intensity = Array{source.intensity_format}(undef, meta.int_asset.encoded_length) - # Validate offsets before reading - file_size = filesize(source.ibd_handle) + # Note: No lock() required here because we are using a thread-local handle! + seek(handle, meta.mz_asset.offset) + read!(handle, mz) - mz_end = meta.mz_asset.offset + sizeof(source.mz_format) * meta.mz_asset.encoded_length - if meta.mz_asset.offset < 0 || mz_end > file_size - throw(FileFormatError("Invalid m/z data offset/length for spectrum $(meta.id): offset=$(meta.mz_asset.offset), end=$mz_end, file_size=$file_size")) - end + seek(handle, meta.int_asset.offset) + read!(handle, intensity) - int_end = meta.int_asset.offset + sizeof(source.intensity_format) * meta.int_asset.encoded_length - if meta.int_asset.offset < 0 || int_end > file_size - throw(FileFormatError("Invalid intensity data offset/length for spectrum $(meta.id): offset=$(meta.int_asset.offset), end=$int_end, file_size=$file_size")) - end - - # Use the new atomic read_at! method for thread-safety - read_at!(source.ibd_handle, mz, meta.mz_asset.offset) - read_at!(source.ibd_handle, intensity, meta.int_asset.offset) - - # imzML data is little-endian. Convert to host byte order. mz .= ltoh.(mz) intensity .= ltoh.(intensity) validate_spectrum_data(mz, intensity, meta.id) - return mz, intensity end @@ -683,6 +758,87 @@ function read_spectrum_from_disk(data::MSIData, source::MzMLSource, meta::Spectr return mz, intensity end +# --- Metadata Caching (Sprint 1: Milestone 4) --- # + +""" + save_metadata_cache(data::MSIData, cache_path::String) + +Serializes the spectrum metadata to a custom binary format for near-instant loading. +""" +function save_metadata_cache(data::MSIData, cache_path::String) + open(cache_path, "w") do io + # Write magic number and version (v2 adds min_mz/max_mz to SpectrumMetadataBinary) + write(io, "JMSI") + write(io, Int32(2)) + + # Write number of spectra + num_spectra = length(data.spectra_metadata) + write(io, Int32(num_spectra)) + + # Write global formats (assuming uniform for now) + # We'll write the names of the types as strings for safety + write(io, string(data.source.mz_format)) + write(io, "\n") + write(io, string(data.source.intensity_format)) + write(io, "\n") + + # Convert to binary structs and write in one block + binary_metadata = Vector{SpectrumMetadataBinary}(undef, num_spectra) + for i in 1:num_spectra + m = data.spectra_metadata[i] + binary_metadata[i] = SpectrumMetadataBinary( + m.x, m.y, Int8(m.mode), + m.mz_asset.offset, m.mz_asset.encoded_length, + m.int_asset.offset, m.int_asset.encoded_length, + Float32(m.mz_asset.min_val), Float32(m.mz_asset.max_val) + ) + end + write(io, binary_metadata) + end + @debug "Metadata cache saved to $cache_path" +end + +""" + load_metadata_cache(cache_path::String, mz_format::Type, int_format::Type) -> Vector{SpectrumMetadata} + +Loads spectrum metadata from a custom binary cache file. +""" +function load_metadata_cache(cache_path::String, mz_format::Type, int_format::Type) + open(cache_path, "r") do io + magic = read(io, 4) + if String(magic) != "JMSI" + error("Invalid cache file format.") + end + version = read(io, Int32) + if version != 2 + error("Unsupported cache version $version (expected 2). Delete the .cache file to regenerate.") + end + num_spectra = read(io, Int32) + + # Skip format strings (we already have them from the header or caller) + readline(io) + readline(io) + + # Read all binary metadata in one swoop + binary_metadata = Vector{SpectrumMetadataBinary}(undef, num_spectra) + read!(io, binary_metadata) + + # Convert back to SpectrumMetadata + spectra_metadata = Vector{SpectrumMetadata}(undef, num_spectra) + for i in 1:num_spectra + b = binary_metadata[i] + mz_asset = SpectrumAsset(mz_format, false, b.mz_offset, b.mz_encoded_len, :mz, Float64(b.min_mz), Float64(b.max_mz)) + int_asset = SpectrumAsset(int_format, false, b.int_offset, b.int_encoded_len, :intensity, 0.0, 0.0) + + spectra_metadata[i] = SpectrumMetadata( + b.x, b.y, "", :sample, SpectrumMode(b.mode), + mz_asset, int_asset + ) + end + return spectra_metadata + end +end + # --- Public API --- # """ @@ -921,7 +1077,7 @@ function precompute_analytics(msi_data::MSIData) println("Processing chunk $chunk_start - $chunk_end / $num_spectra") # Process current chunk - _iterate_spectra_fast(msi_data, collect(chunk_range)) do idx, mz, intensity + _iterate_spectra_fast(msi_data, chunk_range) do idx, mz, intensity # Store metadata modes[idx] = msi_data.spectra_metadata[idx].mode @@ -1400,87 +1556,98 @@ end # --- High-performance Internal Iterator --- # """ - read_compressed_array(io::IO, asset::SpectrumAsset, format::Type) + read_compressed_array(data::MSIData, io::IO, asset::SpectrumAsset, ::Type{T}) where {T} Reads a single data array (m/z or intensity) from an `.ibd` file stream, handling both compressed and uncompressed data. -This is an internal function designed for high-performance iteration. It assumes -the file stream `io` is already positioned at the correct offset. - -- If `asset.is_compressed` is true, it reads `asset.encoded_length` bytes of - compressed data, inflates them using zlib, and reinterprets the result as a - vector of the given `format`. -- If false, it reads `asset.encoded_length` *elements* of uncompressed data - directly into a vector. +This is an internal function designed for high-performance reading and decompressing of binary arrays. +Uses type parameters and buffer pooling to minimize allocations and maximize speed. # Arguments +- `data`: The `MSIData` object. - `io`: The IO stream of the `.ibd` file. - `asset`: The `SpectrumAsset` for the array. -- `format`: The data type of the elements in the array. +- `::Type{T}`: The target format of the data. # Returns -- A `Vector` containing the data. +- A `Vector{T}` containing the data. # Throws -- An error if zlib decompression fails, which can indicate corrupt data or - an incorrect offset in the `.imzML` metadata. +- An error if zlib decompression fails. """ -function read_compressed_array(data::MSIData, io::IO, asset::SpectrumAsset, format::Type) +function read_compressed_array(data::MSIData, io::IO, asset::SpectrumAsset, ::Type{T}) where {T} # Add validation before seeking if asset.offset < 0 || asset.offset >= filesize(io) throw(FileFormatError("Invalid asset offset: $(asset.offset) for file size $(filesize(io))")) end - seek(io, asset.offset) - - if asset.is_compressed - # Get buffer for compressed bytes - compressed_bytes_buffer = get_buffer!(data.buffer_pool, asset.encoded_length) - readbytes!(io, compressed_bytes_buffer, asset.encoded_length) - - println("DEBUG: Decompressing data - offset=$(asset.offset), compressed_bytes=$(length(compressed_bytes_buffer))") - - local decompressed_bytes_buffer - try - # Estimate decompressed size (can be larger than compressed) - # A common heuristic is 4x compressed size, but zlib can be more efficient - # For now, let Libz.inflate handle allocation, then copy to pooled buffer - # This is a temporary allocation, will be optimized later if needed - temp_decompressed = Libz.inflate(compressed_bytes_buffer) - - decompressed_bytes_buffer = get_buffer!(data.buffer_pool, length(temp_decompressed)) - copyto!(decompressed_bytes_buffer, temp_decompressed) + # Optimization: Use Mmap if available to avoid seek and copy + mmap_data = (data.source isa ImzMLSource) ? data.source.mmap_data : nothing - println("DEBUG: Decompression successful - decompressed_bytes=$(length(decompressed_bytes_buffer))") - catch e - @error "ZLIB DECOMPRESSION FAILED. This is likely due to an incorrect offset or corrupt data in the .ibd file." - @error "Asset offset: $(asset.offset), Encoded length: $(asset.encoded_length)" - # Print first 16 bytes to stderr for diagnosis - bytes_to_print = min(16, length(compressed_bytes_buffer)) - @error "First $bytes_to_print bytes of the data chunk we tried to decompress:" - println(stderr, view(compressed_bytes_buffer, 1:bytes_to_print)) - rethrow(e) - finally - release_buffer!(data.buffer_pool, compressed_bytes_buffer) + if asset.is_compressed + local decompressed_view + if mmap_data !== nothing + # Zero-copy access to the compressed segment + # Base64 should be read using a view as well + compressed_view = view(mmap_data, (asset.offset + 1):(asset.offset + asset.encoded_length)) + decompressed_view = Libz.inflate(compressed_view) + else + # Fallback to standard IO + seek(io, asset.offset) + compressed_bytes_buffer = get_buffer!(data.buffer_pool, Int(asset.encoded_length)) + try + readbytes!(io, compressed_bytes_buffer, asset.encoded_length) + decompressed_view = Libz.inflate(view(compressed_bytes_buffer, 1:asset.encoded_length)) + finally + release_buffer!(data.buffer_pool, compressed_bytes_buffer) + end end - # Use an IOBuffer to safely read the data - bytes_io = IOBuffer(decompressed_bytes_buffer) - n_elements = bytes_io.size ÷ sizeof(format) - array = Array{format}(undef, n_elements) - read!(bytes_io, array) + local array + try + # Pre-allocate the typed output array + n_elements = length(decompressed_view) ÷ sizeof(T) + array = Vector{T}(undef, n_elements) + + # Use unsafe_copyto! for zero-overhead copy into the typed array + unsafe_copyto!(reinterpret(Ptr{UInt8}, pointer(array)), pointer(decompressed_view), length(decompressed_view)) + + catch e + @error "ZLIB DECOMPRESSION FAILED at offset $(asset.offset)" + rethrow(e) + end - release_buffer!(data.buffer_pool, decompressed_bytes_buffer) return array else - # Read uncompressed data directly - # For uncompressed imzML, encoded_length is the number of elements - array = Vector{format}(undef, asset.encoded_length) - read!(io, array) - return array + # Read uncompressed data + if mmap_data !== nothing + # SAFETY: Check address alignment (sizeof(T) must divide asset.offset) + # Since mmap_data itself is page-aligned, we only check the offset. + alignment = sizeof(T) + if asset.offset % alignment == 0 + # ZERO-COPY Path + end_pos = asset.offset + asset.encoded_length * alignment + raw_view = view(mmap_data, (asset.offset + 1):end_pos) + return Vector{T}(reinterpret(T, raw_view)) + else + # ALIGNMENT FALLBACK: Memory-to-memory copy (safer than reinterpret) + array = Vector{T}(undef, asset.encoded_length) + # Raw copy from mmap to vector + n_bytes = asset.encoded_length * alignment + unsafe_copyto!(reinterpret(Ptr{UInt8}, pointer(array)), pointer(mmap_data, asset.offset + 1), n_bytes) + return array + end + else + # Fallback to standard IO + seek(io, asset.offset) + array = Vector{T}(undef, asset.encoded_length) + read!(io, array) + return array + end end end + """ _iterate_uncompressed_fast(f::Function, data::MSIData, source::ImzMLSource) @@ -1506,44 +1673,6 @@ _iterate_uncompressed_fast(data, 1) do mz, intensity end ``` """ -function _iterate_uncompressed_fast(f::Function, data::MSIData, source::ImzMLSource, indices_to_iterate::Union{AbstractVector{Int}, Nothing}) - # Optimized path for uncompressed data using buffer reuse - max_points = maximum(meta -> meta.mz_asset.encoded_length, data.spectra_metadata) - mz_buffer = Vector{source.mz_format}(undef, max_points) - int_buffer = Vector{source.intensity_format}(undef, max_points) - - # Determine which indices to iterate over - spectrum_indices = (indices_to_iterate === nothing) ? (1:length(data.spectra_metadata)) : indices_to_iterate - - for i in spectrum_indices - meta = data.spectra_metadata[i] - nPoints = meta.mz_asset.encoded_length - - if nPoints == 0 - f(i, view(mz_buffer, 0:-1), view(int_buffer, 0:-1)) - continue - end - - mz_view = view(mz_buffer, 1:nPoints) - int_view = view(int_buffer, 1:nPoints) - - if meta.mz_asset.offset < meta.int_asset.offset - seek(source.ibd_handle, meta.mz_asset.offset) - read!(source.ibd_handle, mz_view) - seek(source.ibd_handle, meta.int_asset.offset) # FIX: Added missing seek - read!(source.ibd_handle, int_view) - else - seek(source.ibd_handle, meta.int_asset.offset) - read!(source.ibd_handle, int_view) - seek(source.ibd_handle, meta.mz_asset.offset) # FIX: Added missing seek - read!(source.ibd_handle, mz_view) - end - - mz_view .= ltoh.(mz_view) - int_view .= ltoh.(int_view) - f(i, mz_view, int_view) - end -end """ _iterate_compressed_fast(f::Function, data::MSIData, source::ImzMLSource) @@ -1573,11 +1702,19 @@ end """ function _iterate_compressed_fast(f::Function, data::MSIData, source::ImzMLSource, indices_to_iterate::Union{AbstractVector{Int}, Nothing}) # Path for datasets containing at least one compressed spectrum. - # This path reads and decompresses each spectrum individually. + # Optimized to minimize allocations by reusing buffers. # Determine which indices to iterate over spectrum_indices = (indices_to_iterate === nothing) ? (1:length(data.spectra_metadata)) : indices_to_iterate + # Pre-allocate large enough buffers for the expected maximum number of points + # We estimate based on metadata if possible, or grow dynamically + max_encoded = maximum(meta -> max(meta.mz_asset.encoded_length, meta.int_asset.encoded_length), data.spectra_metadata) + + # Heuristic: decompressed size is usually larger. We'll start with 10x and grow if needed. + # But read_compressed_array currently returns a Vector, so we'll need to modify it + # to accept an optional target buffer. + for i in spectrum_indices meta = data.spectra_metadata[i] @@ -1586,7 +1723,9 @@ function _iterate_compressed_fast(f::Function, data::MSIData, source::ImzMLSourc continue end - # Read and decompress each array + # For compressed data, we currently allocate new arrays per spectrum. + # To truly minimize allocations, we'd need read_compressed_array! (in-place version). + # For now, let's ensure we are at least using the optimized type-stable version. mz_array = read_compressed_array(data, source.ibd_handle, meta.mz_asset, source.mz_format) intensity_array = read_compressed_array(data, source.ibd_handle, meta.int_asset, source.intensity_format) @@ -1630,11 +1769,8 @@ function _iterate_spectra_fast_impl(f::Function, data::MSIData, source::ImzMLSou return end - # Check if ANY spectra are compressed and dispatch to the appropriate implementation - any_compressed = any(meta -> meta.mz_asset.is_compressed || meta.int_asset.is_compressed, - data.spectra_metadata) - - if any_compressed + # Use cached compression status for zero-allocation dispatch + if source.is_any_compressed _iterate_compressed_fast(f, data, source, indices_to_iterate) else _iterate_uncompressed_fast(f, data, source, indices_to_iterate) @@ -1667,28 +1803,59 @@ end ``` """ function _iterate_spectra_fast_impl(f::Function, data::MSIData, source::MzMLSource, indices_to_iterate::Union{AbstractVector{Int}, Nothing}) - # This implementation is for mzML. To improve disk I/O, we can reorder the read - # operations to be as sequential as possible based on their offset in the file. - - # Determine which indices to iterate over spectrum_indices = (indices_to_iterate === nothing) ? (1:length(data.spectra_metadata)) : indices_to_iterate - - # Create a vector of (index, offset) tuples to be sorted indices_with_offsets = [(i, data.spectra_metadata[i].mz_asset.offset) for i in spectrum_indices] - - # Sort by offset to make disk access more sequential sort!(indices_with_offsets, by = x -> x[2]) + handle = get_handle(source) + for (i, _) in indices_with_offsets meta = data.spectra_metadata[i] - - mz = read_binary_vector(data, source.file_handle, meta.mz_asset) - intensity = read_binary_vector(data, source.file_handle, meta.int_asset) - + # For MzML, we still have some allocations due to Base64 decoding, + # but we use the thread-local handle. + mz = read_binary_vector(data, handle, meta.mz_asset) + intensity = read_binary_vector(data, handle, meta.int_asset) f(i, mz, intensity) end end +function _iterate_uncompressed_fast(f::Function, data::MSIData, source::ImzMLSource, indices_to_iterate::Union{AbstractVector{Int}, Nothing}) + spectrum_indices = (indices_to_iterate === nothing) ? (1:length(data.spectra_metadata)) : indices_to_iterate + + # Intialize local buffers for this thread's sequential iteration + mz_buf = Vector{Float64}() + int_buf = Vector{Float64}() + + for i in spectrum_indices + meta = data.spectra_metadata[i] + + # Use Mmap views if available (zero-copy if LE) + if source.mmap_data !== nothing + # Optimized Mmap path (same as read_spectrum_from_disk but potentially avoiding copies) + mz, intensity = read_spectrum_from_disk(source, meta) + f(i, mz, intensity) + else + # Read into our loop buffers to avoid continuous allocation + handle = get_handle(source) + + # Resize buffers if necessary (minimal reallocation) + resize!(mz_buf, meta.mz_asset.encoded_length) + resize!(int_buf, meta.int_asset.encoded_length) + + seek(handle, meta.mz_asset.offset) + read!(handle, mz_buf) + seek(handle, meta.int_asset.offset) + read!(handle, int_buf) + + # Convert in-place if possible + mz_buf .= ltoh.(mz_buf) + int_buf .= ltoh.(int_buf) + + f(i, mz_buf, int_buf) + end + end +end + """ _iterate_spectra_fast_serial(f::Function, data::MSIData, indices_to_iterate=nothing) @@ -1793,32 +1960,17 @@ Each thread gets its own file handle, eliminating contention. - Best for bulk processing operations """ function _iterate_spectra_fast_parallel(f::Function, data::MSIData, indices::AbstractVector) - # Split indices into chunks for each thread - n_chunks = Base.Threads.nthreads() - chunk_size = ceil(Int, length(indices) / n_chunks) - chunks = collect(Iterators.partition(indices, chunk_size)) + n_total = length(indices) + n_threads = Base.Threads.nthreads() - Base.Threads.@threads for chunk in chunks - # Each thread gets its own file handle based on source type - if data.source isa ImzMLSource - local_handle = open(data.source.ibd_handle.path, "r") - local_source = ImzMLSource(local_handle, data.source.mz_format, data.source.intensity_format) - - try - # Use the appropriate implementation with thread-local source - _iterate_spectra_fast_impl(f, data, local_source, chunk) - finally - close(local_handle) - end - elseif data.source isa MzMLSource - local_handle = open(data.source.file_handle.path, "r") - local_source = MzMLSource(local_handle, data.source.mz_format, data.source.intensity_format) - - try - _iterate_spectra_fast_impl(f, data, local_source, chunk) - finally - close(local_handle) - end + # Manual chunking to avoid allocations of Iterators.partition and collect + Base.Threads.@threads for t in 1:n_threads + start_idx = ((t - 1) * n_total ÷ n_threads) + 1 + end_idx = (t * n_total) ÷ n_threads + + if start_idx <= end_idx + chunk = view(indices, start_idx:end_idx) + _iterate_spectra_fast_impl(f, data, data.source, chunk) end end end diff --git a/src/MSI_src.jl b/src/MSI_src.jl index 779d5a2..aab9265 100644 --- a/src/MSI_src.jl +++ b/src/MSI_src.jl @@ -20,6 +20,7 @@ export OpenMSIData, MSIData, _iterate_spectra_fast, validate_spectrum, + get_mz_slice, REGISTRY_LOCK # Define shared registry lock @@ -60,18 +61,33 @@ export apply_baseline_correction, apply_intensity_transformation, save_feature_matrix +# Sprint 2: Streaming Pipeline API +export process_dataset!, + PipelineConfig, + StreamingStep, + normalize_inplace!, + transform_inplace!, + smooth_inplace!, + baseline_subtract_inplace!, + detect_peaks_streaming, + calibrate_inplace! + # Include all source files directly into the main module include("BloomFilters.jl") include("Common.jl") +include("ResourcePool.jl") include("MSIData.jl") include("ParserHelpers.jl") include("mzML.jl") include("imzML.jl") include("MzmlConverter.jl") include("Preprocessing.jl") +include("FusedPipeline.jl") include("ImageProcessing.jl") include("Precalculations.jl") include("PreprocessingPipeline.jl") +include("StreamingKernels.jl") +include("StreamingPipeline.jl") using Setfield # For immutable struct updates diff --git a/src/Precalculations.jl b/src/Precalculations.jl index 786c4af..f648aae 100644 --- a/src/Precalculations.jl +++ b/src/Precalculations.jl @@ -1314,18 +1314,16 @@ function _fit_gaussian_and_r2(mz::AbstractVector{<:Real}, intensity::AbstractVec end_idx = min(n, peak_idx + half_window) # Ensure there's enough data to fit - if (end_idx - start_idx + 1) < 3 + count = end_idx - start_idx + 1 + if count < 3 return 0.0 end - x_data = mz[start_idx:end_idx] - y_data = intensity[start_idx:end_idx] - # Estimate Gaussian parameters # Amplitude (A): peak intensity - A_est = intensity[peak_idx] + A_est = float(intensity[peak_idx]) # Mean (μ): m/z at peak intensity - mu_est = mz[peak_idx] + mu_est = float(mz[peak_idx]) # Standard deviation (σ): related to FWHM. FWHM = 2 * sqrt(2 * ln(2)) * σ ≈ 2.355 * σ # So, σ ≈ FWHM / 2.355 fwhm_delta_m = _calculate_fwhm_delta_m(mz, intensity, peak_idx) @@ -1339,19 +1337,27 @@ function _fit_gaussian_and_r2(mz::AbstractVector{<:Real}, intensity::AbstractVec return 0.0 end - # Gaussian function - gaussian(x, A, mu, sigma) = A * exp.(-(x .- mu).^2 ./ (2 * sigma^2)) + # Calculate SS_res, mean_y in a single pass to avoid allocations + SS_res = 0.0 + sum_y = 0.0 + + @inbounds for i in start_idx:end_idx + x_val = float(mz[i]) + y_val = float(intensity[i]) + + # Gaussian function estimate + y_est = A_est * exp(-((x_val - mu_est)^2) / (2 * sigma_est^2)) + + SS_res += (y_val - y_est)^2 + sum_y += y_val + end - # Generate estimated Gaussian curve - y_est = gaussian(x_data, A_est, mu_est, sigma_est) - - # Calculate pseudo R-squared - # R^2 = 1 - (SS_res / SS_tot) - # SS_res = sum((y_data - y_est).^2) - # SS_tot = sum((y_data - mean(y_data)).^2) - - SS_res = sum((y_data .- y_est).^2) - SS_tot = sum((y_data .- mean(y_data)).^2) + mean_y = sum_y / count + SS_tot = 0.0 + + @inbounds for i in start_idx:end_idx + SS_tot += (float(intensity[i]) - mean_y)^2 + end if SS_tot == 0 return 1.0 # Perfect fit if all y_data are the same diff --git a/src/Preprocessing.jl b/src/Preprocessing.jl index 95cb15d..ecf977f 100644 --- a/src/Preprocessing.jl +++ b/src/Preprocessing.jl @@ -12,6 +12,7 @@ generation. # ============================================================================= using Statistics # For mean, median +using SparseArrays using StatsBase # For mad (Median Absolute Deviation) using SavitzkyGolay # For SavitzkyGolay filtering using Dates # For now() @@ -40,7 +41,7 @@ A struct to hold the final feature matrix generated from the preprocessing pipel - `sample_ids::Vector{Int}`: A vector of identifiers for each sample (row) in the `matrix`. """ struct FeatureMatrix - matrix::Array{Float64,2} + matrix::AbstractMatrix{Float64} mz_bins::Vector{Tuple{Float64,Float64}} sample_ids::Vector{Int} end @@ -487,32 +488,23 @@ Estimates the baseline of a spectrum using the SNIP algorithm (internal implemen function _snip_baseline_impl(y::AbstractVector{<:Real}; iterations::Int=100) n = length(y) - # Initialize two buffers. b1 holds the current baseline estimate, b2 for the next. - # Always convert to Float64 to ensure type stability and avoid copying if already correct type + # Initialize the baseline estimate array once b1 = collect(float.(y)) - b2 = similar(b1) - current_b = b1 - next_b = b2 - for k in 1:iterations - # Calculate next baseline estimate into `next_b` based on `current_b` - # Boundary conditions - if n > 1 - next_b[1] = min(current_b[1], current_b[2]) - next_b[n] = min(current_b[n], current_b[n-1]) - end - - @inbounds for i in 2:n-1 - next_b[i] = min(current_b[i], 0.5 * (current_b[i-1] + current_b[i+1])) - end + prev_val = b1[1] + b1[1] = min(b1[1], b1[2]) - # Swap references for the next iteration (no data copy here) - current_b, next_b = next_b, current_b + @inbounds for i in 2:n-1 + curr_val = b1[i] + b1[i] = min(curr_val, 0.5 * (prev_val + b1[i+1])) + prev_val = curr_val + end + b1[n] = min(b1[n], prev_val) end - # Return the final baseline estimate (which is in current_b after the last swap) - return current_b + # Return the final baseline estimate + return b1 end """ @@ -720,18 +712,32 @@ function detect_peaks_profile_core(mz::AbstractVector{<:Real}, y::AbstractVector n = length(y) n < 3 && return NamedTuple{(:mz, :intensity, :fwhm, :shape_r2, :snr, :prominence), Tuple{Float64, Float64, Float64, Float64, Float64, Float64}}[] - noise_level = mad(y, normalize=true) + eps(Float64) - ys = smooth_spectrum_core(y; method=:savitzky_golay, window=max(5, 2*half_window+1), order=2) # Use smoothed data for detection + # Fast, non-allocating noise estimation + mean_y = sum(y) / n + noise_level = (sum(abs.(y .- mean_y)) / n) * 1.5 + eps(Float64) + + ys = smooth_spectrum_core(y; method=:savitzky_golay, window=max(5, 2*half_window+1), order=2) candidate_peak_indices = Int[] - for i in 2:n-1 + sizehint!(candidate_peak_indices, div(n, 10)) # Pre-allocate memory capacity + + @inbounds for i in 2:n-1 left = max(1, i - half_window) right = min(n, i + half_window) - # Prominence check - prominence = ys[i] - max(minimum(@view ys[left:i]), minimum(@view ys[i:right])) + # Avoid @view allocation in tight loop by manually computing minimums and maximums + min_left = ys[left] + for j in left:i; min_left = min(min_left, ys[j]); end - if ys[i] >= maximum(@view ys[left:right]) && + min_right = ys[i] + for j in i:right; min_right = min(min_right, ys[j]); end + + prominence = ys[i] - max(min_left, min_right) + + max_local = ys[left] + for j in left:right; max_local = max(max_local, ys[j]); end + + if ys[i] >= max_local && (ys[i] > snr_threshold * noise_level) && (prominence > min_peak_prominence * ys[i]) push!(candidate_peak_indices, i) @@ -768,7 +774,14 @@ function detect_peaks_profile_core(mz::AbstractVector{<:Real}, y::AbstractVector left = max(1, p_idx - half_window) right = min(n, p_idx + half_window) - prominence = ys[p_idx] - max(minimum(@view ys[left:p_idx]), minimum(@view ys[p_idx:right])) + + min_left = ys[left] + for j in left:p_idx; min_left = min(min_left, ys[j]); end + + min_right = ys[p_idx] + for j in p_idx:right; min_right = min(min_right, ys[j]); end + + prominence = ys[p_idx] - max(min_left, min_right) push!(detected_peaks, (mz=peak_mz, intensity=peak_int, fwhm=fwhm_ppm, shape_r2=shape_r2, snr=peak_snr, prominence=prominence)) end diff --git a/src/ResourcePool.jl b/src/ResourcePool.jl new file mode 100644 index 0000000..c0e0413 --- /dev/null +++ b/src/ResourcePool.jl @@ -0,0 +1,95 @@ +# src/ResourcePool.jl +using Base.Threads + +""" + ResourcePool{T} + +A thread-safe pool for reusing objects of type `T` to minimize allocations and GC pressure. +Specifically designed for high-performance computing tasks where large buffers are needed +repeatedly across multiple threads. + +# Fields: +- `pool::Vector{T}`: The underlying storage for idle resources. +- `lock::ReentrantLock`: Ensures thread-safe access to the pool. +- `max_size::Int`: Maximum number of resources to hold in the pool. +- `constructor::Function`: A function to create a new resource if the pool is empty. +""" +mutable struct ResourcePool{T} + pool::Vector{T} + lock::ReentrantLock + max_size::Int + constructor::Function +end + +""" + aligned_vector(::Type{T}, n::Int; alignment::Int=64) where T + +Creates a `Vector{T}` that is aligned to `alignment` bytes. +Note: In modern Julia, standard vectors are often 16 or 64 byte aligned, but for +HPC we ensure this by allocating slightly more and using a view, or using +specific pointers. For simplicity and performance, we use a small hack: +allocating a larger array and taking a 64-byte aligned view. +""" +function aligned_vector(::Type{T}, n::Int; alignment::Int=64) where T + # Allocate enough space to find an aligned starting point + raw = Vector{UInt8}(undef, n * sizeof(T) + alignment) + ptr = Int(pointer(raw)) + off = (alignment - (ptr % alignment)) % alignment + # Return a reinterpret view of the aligned segment + return reinterpret(T, view(raw, (off + 1):(off + n * sizeof(T)))) +end + +""" + ResourcePool{T}(constructor::Function; max_size::Int=2 * nthreads()) + +Creates a new `ResourcePool` for resources of type `T`. +""" +function ResourcePool{T}(constructor::Function; max_size::Int=2 * nthreads()) where T + return ResourcePool{T}(T[], ReentrantLock(), max_size, constructor) +end + +""" + acquire(pool::ResourcePool{T}) -> T + +Retrieves a resource from the pool. If the pool is empty, a new resource is created +using the constructor. +""" +function acquire(pool::ResourcePool{T}) where T + lock(pool.lock) do + if !isempty(pool.pool) + return pop!(pool.pool) + end + end + # Create new resource outside of lock to minimize contention + return pool.constructor() +end + +""" + release!(pool::ResourcePool{T}, resource::T) + +Returns a resource to the pool for later reuse. If the pool is already at `max_size`, +the resource is allowed to be garbage collected. +""" +function release!(pool::ResourcePool{T}, resource::T) where T + lock(pool.lock) do + if length(pool.pool) < pool.max_size + push!(pool.pool, resource) + end + end + return nothing +end + +""" + with_resource(f::Function, pool::ResourcePool{T}) + +Acquires a resource from the pool, executes the function `f(resource)`, and +automatically releases the resource back to the pool when finished. +""" +function with_resource(f::Function, pool::ResourcePool{T}) where T + resource = acquire(pool) + try + return f(resource) + finally + release!(pool, resource) + end +end diff --git a/src/StreamingKernels.jl b/src/StreamingKernels.jl new file mode 100644 index 0000000..8f4b1a0 --- /dev/null +++ b/src/StreamingKernels.jl @@ -0,0 +1,317 @@ +# src/StreamingKernels.jl +# ============================================================================ +# In-Place Spectral Kernels for the Streaming Pipeline +# +# These functions operate on raw (mz, intensity) views from the Sprint 1 +# Mmap engine. They write results back to the input buffers using .= to +# achieve zero-allocation processing per spectrum. +# +# Design contract: +# - All !-suffixed functions modify their arguments in-place +# - If a kernel needs temporary storage, it borrows from data.resource_pool +# - No function creates MutableSpectrum objects +# ============================================================================ + +using Statistics: mean, median + +# ============================================================================= +# Category A: Purely Streamable Kernels +# ============================================================================= + +""" + normalize_inplace!(intensity::AbstractVector{<:Real}, method::Symbol) + +Normalizes intensity values in-place. Supports :tic, :median, :rms. +Zero-allocation for the normalization itself. +""" +@inline function normalize_inplace!(intensity::AbstractVector{<:Real}, method::Symbol) + if method === :tic + s = sum(intensity) + if s > 0 + intensity ./= s + end + elseif method === :median + m = median(intensity) + if m > 0 + intensity ./= m + end + elseif method === :rms + s = sqrt(sum(abs2, intensity) / length(intensity)) + if s > 0 + intensity ./= s + end + end + return intensity +end + +""" + transform_inplace!(intensity::AbstractVector{Float64}, method::Symbol) + +Applies intensity transformation in-place. Supports :sqrt, :log1p, :log, :log2, :log10. +""" +@inline function transform_inplace!(intensity::AbstractVector{Float64}, method::Symbol) + if method === :sqrt + @inbounds @simd for i in eachindex(intensity) + intensity[i] = sqrt(max(0.0, intensity[i])) + end + elseif method === :log1p + @inbounds @simd for i in eachindex(intensity) + intensity[i] = log1p(max(0.0, intensity[i])) + end + elseif method === :log + @inbounds @simd for i in eachindex(intensity) + intensity[i] = log(max(eps(Float64), intensity[i])) + end + elseif method === :log2 + @inbounds @simd for i in eachindex(intensity) + intensity[i] = log2(max(eps(Float64), intensity[i])) + end + elseif method === :log10 + @inbounds @simd for i in eachindex(intensity) + intensity[i] = log10(max(eps(Float64), intensity[i])) + end + end + return intensity +end + +""" + smooth_inplace!(intensity::AbstractVector{Float64}, data::MSIData; + method::Symbol=:savitzky_golay, window::Int=9, order::Int=2) + +Smooths intensity in-place using a temporary buffer from the resource pool. +The SavitzkyGolay library allocates internally, but we copy the result back +to the original buffer and return the pool buffer. +""" +function smooth_inplace!(intensity::AbstractVector{Float64}, scratch::AbstractVector{Float64}, data::MSIData; + method::Symbol=:savitzky_golay, window::Int=9, order::Int=2) + n = length(intensity) + if n < 3 + return intensity + end + + if method === :savitzky_golay + win = isodd(window) ? window : window + 1 + if n < win + return intensity + end + # SavitzkyGolay handles its own math but causes mild allocation. + res = SavitzkyGolay.savitzky_golay(collect(intensity), win, order) + @inbounds for i in eachindex(intensity) + intensity[i] = max(0.0, res.y[i]) + end + elseif method === :moving_average + copyto!(scratch, intensity) + + half_w = div(window, 2) + @inbounds for i in 1:n + s_idx = max(1, i - half_w) + e_idx = min(n, i + half_w) + s = 0.0 + @simd for j in s_idx:e_idx + s += scratch[j] + end + intensity[i] = max(0.0, s / (e_idx - s_idx + 1)) + end + end + return intensity +end + +""" + baseline_subtract_inplace!(intensity::AbstractVector{Float64}, data::MSIData; + method::Symbol=:snip, iterations::Int=100, window::Int=20) + +Subtracts baseline from intensity in-place. Uses two pool buffers for the +SNIP ping-pong iteration to avoid any heap allocation in the hot loop. +""" +function baseline_subtract_inplace!(intensity::AbstractVector{Float64}, scratch::AbstractVector{Float64}, data::MSIData; + method::Symbol=:snip, iterations::Int=100, window::Int=20) + n = length(intensity) + if n < 3 + return intensity + end + + if method === :snip + copyto!(scratch, intensity) + + for k in 1:iterations + prev_val = scratch[1] + scratch[1] = min(scratch[1], scratch[2]) + + @inbounds for i in 2:n-1 + curr_val = scratch[i] + scratch[i] = min(curr_val, 0.5 * (prev_val + scratch[i+1])) + prev_val = curr_val + end + scratch[n] = min(scratch[n], prev_val) + end + + @inbounds @simd for i in 1:n + intensity[i] = max(0.0, intensity[i] - scratch[i]) + end + elseif method === :convex_hull + baseline = convex_hull_baseline(intensity) + @inbounds @simd for i in eachindex(intensity) + intensity[i] = max(0.0, intensity[i] - baseline[i]) + end + elseif method === :median + baseline = median_baseline(intensity; window=window) + @inbounds @simd for i in eachindex(intensity) + intensity[i] = max(0.0, intensity[i] - baseline[i]) + end + end + return intensity +end + +""" + detect_peaks_streaming(mz::AbstractVector, intensity::AbstractVector; + method::Symbol=:profile, snr_threshold::Float64=3.0, + half_window::Int=10, min_peak_prominence::Float64=0.1, + merge_peaks_tolerance::Float64=0.002) + +Detects peaks and returns a vector of (mz, intensity) tuples. +This delegates to existing _core functions but returns a lightweight format +suitable for sparse accumulation (no NamedTuple overhead in the hot path). +""" +function detect_peaks_streaming(callback::Function, mz::AbstractVector{Float64}, intensity::AbstractVector{Float64}, scratch::AbstractVector{Float64}; + method::Symbol=:profile, snr_threshold::Float64=3.0, + half_window::Int=10, min_peak_prominence::Float64=0.1, + merge_peaks_tolerance::Float64=0.002) + n = length(intensity) + if n < 3 + return + end + + if method === :profile || method === :wavelet + # Zero-allocation noisy estimation (using mean of bottom half) + sum_i = 0.0 + @simd for i in 1:n + sum_i += intensity[i] + end + mean_i = sum_i / n + + sum_noise = 0.0 + count_noise = 0 + @inbounds for i in 1:n + if intensity[i] < mean_i + sum_noise += intensity[i] + count_noise += 1 + end + end + # Use * 1.5 as an approximation to MAD + noise_level = count_noise > 0 ? (sum_noise / count_noise) * 1.5 + eps(Float64) : mean_i + eps(Float64) + + # We will use the scratch buffer to store candidate indices to avoid allocating `Int[]` + # Because scratch is Float64, we can safely store integer indices up to 2^53 exactly. + num_candidates = 0 + + @inbounds for i in 2:n-1 + if intensity[i] > snr_threshold * noise_level + left = max(1, i - half_window) + right = min(n, i + half_window) + + is_max = true + for j in left:right + if intensity[j] > intensity[i] + is_max = false + break + end + end + + if is_max + # Compute prominence + min_left = intensity[i] + for j in left:i + if intensity[j] < min_left + min_left = intensity[j] + end + end + min_right = intensity[i] + for j in i:right + if intensity[j] < min_right + min_right = intensity[j] + end + end + + prominence = intensity[i] - max(min_left, min_right) + + if prominence > min_peak_prominence * intensity[i] + num_candidates += 1 + scratch[num_candidates] = i + end + end + end + end + + # Merge close peaks + if num_candidates > 0 + if merge_peaks_tolerance > 0 + last_idx = trunc(Int, scratch[1]) + # We emit the first peak lazily down below, so let's compact them in place + num_merged = 1 + + for i in 2:num_candidates + idx = trunc(Int, scratch[i]) + if (mz[idx] - mz[last_idx]) > merge_peaks_tolerance + num_merged += 1 + scratch[num_merged] = idx + last_idx = idx + elseif intensity[idx] > intensity[last_idx] + scratch[num_merged] = idx + last_idx = idx + end + end + num_candidates = num_merged + end + + # Emit merged peaks + for i in 1:num_candidates + idx = trunc(Int, scratch[i]) + callback(mz[idx], intensity[idx]) + end + end + + elseif method === :centroid + # Just use any value over snr_threshold * mean_noise + sum_i = sum(intensity) + mean_i = sum_i / n + noise_level = mean_i + eps(Float64) + + @inbounds for i in 1:n + if intensity[i] > snr_threshold * noise_level + callback(mz[i], intensity[i]) + end + end + end +end + +# ============================================================================= +# Category B: Conditionally Streamable Kernels (Fixed-Reference) +# ============================================================================= + +""" + calibrate_inplace!(mz::Vector{Float64}, intensity::AbstractVector, + reference_masses::Vector{Float64}; ppm_tolerance::Float64=20.0) + +Calibrates the m/z axis in-place using a fixed dictionary of internal standard +reference masses. This is streamable because the reference is constant. + +Returns `true` if calibration was applied, `false` if insufficient peaks were found. +""" +function calibrate_inplace!(mz::Vector{Float64}, intensity::AbstractVector, + reference_masses::Vector{Float64}; ppm_tolerance::Float64=20.0) + matched_peaks = find_calibration_peaks_core(mz, intensity, reference_masses; + ppm_tolerance=ppm_tolerance) + if length(matched_peaks) < 2 + return false # Insufficient reference peaks + end + + measured = sort(collect(values(matched_peaks))) + theoretical = sort(collect(keys(matched_peaks))) + itp = linear_interpolation(measured, theoretical, extrapolation_bc=Line()) + + # Apply calibration in-place + @inbounds for i in eachindex(mz) + mz[i] = itp(mz[i]) + end + return true +end diff --git a/src/StreamingPipeline.jl b/src/StreamingPipeline.jl new file mode 100644 index 0000000..ee97ac4 --- /dev/null +++ b/src/StreamingPipeline.jl @@ -0,0 +1,402 @@ +# src/StreamingPipeline.jl +# ============================================================================ +# The Streaming Pipeline Executor +# +# This module provides `process_dataset!`, the Sprint 2 master function that +# streams spectral data through an in-place kernel chain and accumulates +# results into a SparseMatrixCSC without ever holding more than 1 spectrum +# per thread in RAM. +# +# Architecture: +# 1. _iterate_spectra_fast → Mmap zero-copy views +# 2. copyto!(writable_buf, view) → makes mutable copy for kernels +# 3. Kernel chain: smooth! → baseline! → peaks → bin +# 4. Thread-local (I, J, V) sparse accumulators +# 5. Final sparse(I, J, V, num_bins, num_spectra) assembly +# +# This works alongside the existing execute_full_preprocessing in +# PreprocessingPipeline.jl — it does NOT replace the app.jl integration. +# ============================================================================ + +using SparseArrays +using Printf + +# ============================================================================= +# Configuration Structs +# ============================================================================= + +""" + StreamingStep + +Represents a single step in the streaming pipeline. +""" +struct StreamingStep + name::Symbol + params::Dict{Symbol, Any} +end + +""" + PipelineConfig + +Holds the complete configuration for a streaming pipeline execution. + +# Fields +- `steps::Vector{StreamingStep}` — ordered sequence of processing steps +- `reference_peaks::Vector{Float64}` — fixed m/z values for calibration (Category B) +- `num_bins::Int` — number of bins for the output feature matrix +- `min_peaks_per_bin::Int` — minimum peak count to keep a bin +- `frequency_threshold::Float64` — minimum fraction of spectra a bin must appear in (0.0-1.0) + +# Example +```julia +config = PipelineConfig( + steps = [ + StreamingStep(:smoothing, Dict(:method => :savitzky_golay, :window => 9, :order => 2)), + StreamingStep(:baseline_correction, Dict(:method => :snip, :iterations => 100)), + StreamingStep(:normalization, Dict(:method => :tic)), + StreamingStep(:peak_picking, Dict(:method => :profile, :snr_threshold => 3.0)), + ], + num_bins = 2000 +) +``` +""" +struct PipelineConfig + steps::Vector{StreamingStep} + reference_peaks::Vector{Float64} + num_bins::Int + min_peaks_per_bin::Int + frequency_threshold::Float64 +end + +# Convenience constructor with defaults +function PipelineConfig(; steps::Vector{StreamingStep}=StreamingStep[], + reference_peaks::Vector{Float64}=Float64[], + num_bins::Int=2000, + min_peaks_per_bin::Int=3, + frequency_threshold::Float64=0.0) + return PipelineConfig(steps, reference_peaks, num_bins, min_peaks_per_bin, frequency_threshold) +end + +# ============================================================================= +# Sparse Accumulator (Thread-Local) +# ============================================================================= + +""" + SparseAccumulator + +Thread-local accumulator for sparse matrix construction. +Collects (row, col, val) triplets that will be assembled into +a SparseMatrixCSC at the end of the pipeline. +""" +mutable struct SparseAccumulator + I::Vector{Int} # Row indices (bin indices) + J::Vector{Int} # Column indices (spectrum indices) + V::Vector{Float64} # Values (intensities) + lck::Base.Threads.SpinLock + + function SparseAccumulator(capacity_hint::Int=10000) + acc = new( + Vector{Int}(undef, 0), + Vector{Int}(undef, 0), + Vector{Float64}(undef, 0), + Base.Threads.SpinLock() + ) + sizehint!(acc.I, capacity_hint) + sizehint!(acc.J, capacity_hint) + sizehint!(acc.V, capacity_hint) + return acc + end +end + +""" + accumulate!(acc::SparseAccumulator, spectrum_idx::Int, bin_indices::AbstractVector{Int}, + intensities::AbstractVector{Float64}) + +Appends peak data for one spectrum into the sparse accumulator. +""" +@inline function accumulate!(acc::SparseAccumulator, spectrum_idx::Int, + bin_indices::AbstractVector{Int}, + intensities::AbstractVector{Float64}) + n = length(bin_indices) + for k in 1:n + @inbounds begin + push!(acc.I, bin_indices[k]) + push!(acc.J, spectrum_idx) + push!(acc.V, intensities[k]) + end + end +end + +# ============================================================================= +# The Pipeline Executor +# ============================================================================= + +""" + process_dataset!(data::MSIData, config::PipelineConfig; + progress_callback::Union{Function, Nothing}=nothing, + masked_indices::Union{AbstractVector{Int}, Nothing}=nothing) + +The Sprint 2 master streaming function. Processes an entire MSI dataset through +a kernel chain without holding more than 1 spectrum per thread in RAM. + +# Returns +- `SparseMatrixCSC{Float64, Int}`: The feature matrix (bins × spectra) +- `Vector{Float64}`: The m/z bin centers + +# Architecture +1. Ensures analytics are computed (for global m/z range) +2. Creates thread-local SparseAccumulators +3. Streams spectra via `_iterate_spectra_fast` +4. Per spectrum: copy view → kernel chain → peak detect → bin → accumulate +5. Merges accumulators → `sparse(I, J, V)` +""" +function process_dataset!(data::MSIData, config::PipelineConfig; + progress_callback::Union{Function, Nothing}=nothing, + masked_indices::Union{AbstractVector{Int}, Nothing}=nothing) + + # --- Step 1: Ensure analytics are computed (provides global m/z range) --- + if !is_set(data.analytics_ready) + println("Pre-computing analytics for streaming pipeline...") + precompute_analytics(data) + end + + # Determine global m/z range for binning + global_min_mz = Base.Threads.atomic_add!(data.global_min_mz, 0.0) + global_max_mz = Base.Threads.atomic_add!(data.global_max_mz, 0.0) + + if !isfinite(global_min_mz) || !isfinite(global_max_mz) || global_min_mz >= global_max_mz + @warn "Invalid global m/z range: [$global_min_mz, $global_max_mz]. Cannot bin peaks." + return spzeros(0, 0), Float64[] + end + + num_bins = config.num_bins + bin_edges = range(global_min_mz, stop=global_max_mz, length=num_bins + 1) + bin_centers = [(bin_edges[i] + bin_edges[i+1]) / 2 for i in 1:num_bins] + inv_bin_width = 1.0 / step(bin_edges) + + num_spectra = length(data.spectra_metadata) + indices_to_process = masked_indices === nothing ? nothing : masked_indices + + # --- Step 2: Create thread-local accumulators --- + n_threads = Base.Threads.nthreads() + accumulators = [SparseAccumulator(num_spectra * 10) for _ in 1:n_threads] + spectra_processed = Base.Threads.Atomic{Int}(0) + + # NEW: Create dedicated workspace buffers for each thread. + # This completely eliminates the need for acquire/release and prevents deadlocks. + workspaces_mz = [Vector{Float64}(undef, 0) for _ in 1:n_threads] + workspaces_int = [Vector{Float64}(undef, 0) for _ in 1:n_threads] + workspaces_scratch = [Vector{Float64}(undef, 0) for _ in 1:n_threads] + + # Pre-parse step configuration for fast dispatch in the hot loop + has_smoothing = false + has_baseline = false + has_normalization = false + has_transform = false + has_peak_picking = false + has_calibration = false + + smooth_params = Dict{Symbol, Any}() + baseline_params = Dict{Symbol, Any}() + norm_params = Dict{Symbol, Any}() + transform_params = Dict{Symbol, Any}() + peak_params = Dict{Symbol, Any}() + + for s in config.steps + if s.name === :smoothing + has_smoothing = true + smooth_params = s.params + elseif s.name === :baseline_correction + has_baseline = true + baseline_params = s.params + elseif s.name === :normalization + has_normalization = true + norm_params = s.params + elseif s.name === :stabilization || s.name === :intensity_transformation + has_transform = true + transform_params = s.params + elseif s.name === :peak_picking + has_peak_picking = true + peak_params = s.params + elseif s.name === :calibration + has_calibration = true + end + end + + reference_masses = config.reference_peaks + + # --- Step 3: Stream and process --- + start_time = time_ns() + + # Use let block to capture all variables cleanly for the closure + let data=data, accumulators=accumulators, spectra_processed=spectra_processed, + bin_edges=bin_edges, num_bins=num_bins, inv_bin_width=inv_bin_width, + global_min_mz=global_min_mz, + workspaces_mz=workspaces_mz, workspaces_int=workspaces_int, workspaces_scratch=workspaces_scratch, + has_smoothing=has_smoothing, has_baseline=has_baseline, + has_normalization=has_normalization, has_transform=has_transform, + has_peak_picking=has_peak_picking, has_calibration=has_calibration, + smooth_params=smooth_params, baseline_params=baseline_params, + norm_params=norm_params, transform_params=transform_params, + peak_params=peak_params, reference_masses=reference_masses + + _iterate_spectra_fast(data, indices_to_process) do idx, mz_view, int_view + thread_id = Base.Threads.threadid() + acc = accumulators[thread_id] + + # --- Grab Thread-Local Workspaces --- + # No locking, no blocking, guaranteed to be available + mz_buf = workspaces_mz[thread_id] + int_buf = workspaces_int[thread_id] + scratch_buf = workspaces_scratch[thread_id] + + resize!(mz_buf, length(mz_view)) + resize!(int_buf, length(int_view)) + resize!(scratch_buf, length(int_view)) + copyto!(mz_buf, mz_view) + copyto!(int_buf, int_view) + + # --- Kernel Chain (in pipeline order) --- + + # Category B: Fixed-reference calibration + if has_calibration && !isempty(reference_masses) + calibrate_inplace!(mz_buf, int_buf, reference_masses) + end + + # Category A: Intensity transformation + if has_transform + transform_inplace!(int_buf, get(transform_params, :method, :sqrt)) + end + + # Category A: Smoothing + if has_smoothing + smooth_inplace!(int_buf, scratch_buf, data; + method=get(smooth_params, :method, :savitzky_golay), + window=get(smooth_params, :window, 9), + order=get(smooth_params, :order, 2)) + end + + # Category A: Baseline correction + if has_baseline + baseline_subtract_inplace!(int_buf, scratch_buf, data; + method=get(baseline_params, :method, :snip), + iterations=get(baseline_params, :iterations, 100), + window=get(baseline_params, :window, 20)) + end + + # Category A: Normalization + if has_normalization + normalize_inplace!(int_buf, get(norm_params, :method, :tic)) + end + + # --- Peak Detection & Binning --- + if has_peak_picking + detect_peaks_streaming(mz_buf, int_buf, scratch_buf; + method=get(peak_params, :method, :profile), + snr_threshold=Float64(get(peak_params, :snr_threshold, 3.0)), + half_window=Int(get(peak_params, :half_window, 10)), + min_peak_prominence=Float64(get(peak_params, :min_peak_prominence, 0.1)), + merge_peaks_tolerance=Float64(get(peak_params, :merge_peaks_tolerance, 0.002))) do peak_mz, peak_int + + # Bin each discovered peak directly + bin_idx = trunc(Int, (peak_mz - global_min_mz) * inv_bin_width) + 1 + bin_idx = clamp(bin_idx, 1, num_bins) + + push!(acc.I, bin_idx) + push!(acc.J, idx) + push!(acc.V, peak_int) + end + else + # No peak picking: bin raw intensity directly + @inbounds for i in eachindex(mz_buf) + bin_idx = trunc(Int, (mz_buf[i] - global_min_mz) * inv_bin_width) + 1 + bin_idx = clamp(bin_idx, 1, num_bins) + + push!(acc.I, bin_idx) + push!(acc.J, idx) + push!(acc.V, int_buf[i]) + end + end + + Base.Threads.atomic_add!(spectra_processed, 1) + end + end + + # --- Step 4: Merge thread-local accumulators --- + total_entries = sum(length(acc.I) for acc in accumulators) + merged_I = Vector{Int}(undef, total_entries) + merged_J = Vector{Int}(undef, total_entries) + merged_V = Vector{Float64}(undef, total_entries) + + offset = 0 + for acc in accumulators + n = length(acc.I) + if n > 0 + copyto!(merged_I, offset + 1, acc.I, 1, n) + copyto!(merged_J, offset + 1, acc.J, 1, n) + copyto!(merged_V, offset + 1, acc.V, 1, n) + offset += n + end + end + + # --- Step 5: Assemble sparse matrix --- + # Use max combiner: when multiple peaks map to the same bin for same spectrum, + # keep the maximum intensity + feature_matrix = sparse(merged_I, merged_J, merged_V, num_bins, num_spectra, max) + + # --- Step 6: Apply frequency threshold if configured --- + if config.frequency_threshold > 0.0 + # Count how many spectra have a non-zero value in each bin + bin_presence = vec(sum(feature_matrix .> 0, dims=2)) + min_count = ceil(Int, config.frequency_threshold * num_spectra) + keep_bins = findall(bin_presence .>= min_count) + feature_matrix = feature_matrix[keep_bins, :] + bin_centers = bin_centers[keep_bins] + end + + duration = (time_ns() - start_time) / 1e9 + n_processed = spectra_processed[] + n_nonzeros = nnz(feature_matrix) + sparsity = 1.0 - n_nonzeros / (size(feature_matrix, 1) * size(feature_matrix, 2) + 1) + + @printf "Streaming pipeline complete: %d spectra processed in %.2f seconds.\n" n_processed duration + @printf "Feature matrix: %d bins × %d spectra, %d non-zeros (%.1f%% sparse)\n" size(feature_matrix, 1) size(feature_matrix, 2) n_nonzeros sparsity * 100 + @printf "RAM: %.1f MB (vs %.1f MB dense)\n" (n_nonzeros * 16) / 1e6 (size(feature_matrix, 1) * size(feature_matrix, 2) * 8) / 1e6 + + if progress_callback !== nothing + progress_callback(1.0) + end + + return feature_matrix, collect(Float64, bin_centers) +end + +""" + save_sparse_matrix(matrix::SparseMatrixCSC, output_path::String) + +Exports a highly optimized SparseMatrixCSC array to disk using the standard +Matrix Market Coordinate format (`.mtx`), guaranteeing no bottleneck or OOM crashes +for extremely large MS dataset persistence. +""" +function save_sparse_matrix(matrix::SparseMatrixCSC{Float64, Int}, output_path::String) + m, n = size(matrix) + nnz_val = nnz(matrix) + # Use streaming I/O with a large buffer for ultra-fast persistence + open(output_path, "w") do io + # Write Matrix Market Header + write(io, "%%MatrixMarket matrix coordinate real general\n") + write(io, "$m $n $nnz_val\n") + + # Directly extract CSC properties (O(1) memory, zero allocation) + row_indices = rowvals(matrix) + values_array = nonzeros(matrix) + + @inbounds for filter_j in 1:n + # nzrange returns the index bounds for non-zero elements in column 'j' + for idx in nzrange(matrix, filter_j) + i = row_indices[idx] + v = values_array[idx] + write(io, "$i $filter_j $v\n") + end + end + end +end diff --git a/src/imzML.jl b/src/imzML.jl index d173ba6..7b319ea 100644 --- a/src/imzML.jl +++ b/src/imzML.jl @@ -1,5 +1,4 @@ -# src/imzML.jl -using Images, Statistics, CairoMakie, DataFrames, Printf, ColorSchemes, StatsBase +using Images, Statistics, CairoMakie, DataFrames, Printf, ColorSchemes, StatsBase, Mmap """ This file provides a library for parsing `.imzML` and `.ibd` files in pure Julia. @@ -487,8 +486,8 @@ function parse_imzml_spectrum_block(stream::IO, hIbd::Union{IO, ThreadSafeFileHa @warn "Expected spectrum block $k but found none or reached EOF prematurely. Stopping parsing." # Fill remaining spectra_metadata with placeholder or error. for j in k:num_spectra - mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, Int64(0), 0, :mz) - int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, Int64(0), 0, :intensity) + mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, Int64(0), 0, :mz, 0.0, 0.0) + int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, Int64(0), 0, :intensity, 0.0, 0.0) spectra_metadata[j] = SpectrumMetadata(Int32(0), Int32(0), "", :sample, global_mode, mz_asset, int_asset) end break @@ -512,8 +511,8 @@ function parse_imzml_spectrum_block(stream::IO, hIbd::Union{IO, ThreadSafeFileHa if length(mz_data) != 1 || length(int_data) != 1 println("DEBUG: Spectrum $k is empty or invalid - creating placeholder metadata") - mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, Int64(0), 0, :mz) - int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, Int64(0), 0, :intensity) + mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, Int64(0), 0, :mz, 0.0, 0.0) + int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, Int64(0), 0, :intensity, 0.0, 0.0) else mz_info = mz_data[1] int_info = int_data[1] @@ -527,9 +526,9 @@ function parse_imzml_spectrum_block(stream::IO, hIbd::Union{IO, ThreadSafeFileHa end mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, mz_info.offset, - mz_is_compressed ? mz_info.encoded_length : mz_info.array_length, :mz) + mz_is_compressed ? mz_info.encoded_length : mz_info.array_length, :mz, 0.0, 0.0) int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, int_info.offset, - int_is_compressed ? int_info.encoded_length : int_info.array_length, :intensity) + int_is_compressed ? int_info.encoded_length : int_info.array_length, :intensity, 0.0, 0.0) end spectra_metadata[k] = SpectrumMetadata(x, y, "", :sample, spectrum_mode, mz_asset, int_asset) @@ -555,7 +554,7 @@ parsed information acquired by the helper functions. - `msi_data::MSIData`: The MSI data. """ -function load_imzml_lazy(file_path::String; cache_size::Int=100) +function load_imzml_lazy(file_path::String; cache_size::Int=100, use_mmap::Bool=true) println("DEBUG: Checking for .imzML file at $file_path") if !isfile(file_path) throw(FileFormatError("Provided path is not a file: $(file_path)")) @@ -569,110 +568,119 @@ function load_imzml_lazy(file_path::String; cache_size::Int=100) println("DEBUG: Opening file streams for .imzML and .ibd") stream = open(file_path, "r") - ts_hIbd = ThreadSafeFileHandle(ibd_path) + + # --- Handle Pool Optimization --- + # We open multiple handles to the same .ibd file to avoid lock contention in parallel code. + num_handles = Threads.nthreads() + ibd_handles = [open(ibd_path, "r") for _ in 1:num_handles] + + # --- Mmap Optimization with RAM Safety --- + mmap_data = nothing + if use_mmap + try + file_size = filesize(ibd_path) + free_ram = Sys.free_memory() + + if file_size > free_ram * 0.8 + @warn "Dataset size ($(round(file_size/1e9, digits=2)) GB) exceeds 80% of free RAM. Mmap will still work via 'Streaming', but expect slight I/O overhead." + end + + @debug "Memory mapping .ibd file..." + # We use the first handle for mmapping + mmap_data = Mmap.mmap(ibd_handles[1], Vector{UInt8}, file_size) + # Use POSIX shim for sequential access optimization + posix_madvise(mmap_data, MADV_SEQUENTIAL) + @debug ".ibd file mmapped successfully." + catch e + @warn "Memory mapping failed, falling back to standard I/O: $e" + end + end try - # --- NEW: Parse all header information in a more efficient single pass --- - println("DEBUG: Parsing imzML header...") + @debug "Parsing imzML header..." (instrument_meta, param_groups, imgDim) = parse_imzml_header(stream) - # The header parser will have reset the stream for the next step (spectrum parsing) - - println("--- Extracted Instrument Metadata ---") - println("Resolution: ", instrument_meta.resolution) - println("Acquisition Mode (pre-check): ", instrument_meta.acquisition_mode) - println("Calibration Status: ", instrument_meta.calibration_status) - println("Instrument Model: ", instrument_meta.instrument_model) - println("Mass Accuracy (ppm): ", instrument_meta.mass_accuracy_ppm) - println("Laser Settings: ", instrument_meta.laser_settings) - println("Polarity: ", instrument_meta.polarity) - println("------------------------------------") width, height, num_spectra = imgDim - println("DEBUG: Image dimensions: $(width)x$(height), $num_spectra spectra.") + @debug "Image dimensions: $(width)x$(height), $num_spectra spectra." - # Extract default formats from the parsed param_groups + # ... (format extraction logic stays the same) ... + # [Simplified for brevity in replacement chunk, but keeping the logic] mz_group = nothing int_group = nothing - for group in values(param_groups) - if group.Axis == 1 - mz_group = group - elseif group.Axis == 2 - int_group = group + if group.Axis == 1; mz_group = group; elseif group.Axis == 2; int_group = group; end + end + + default_mz_format = (mz_group !== nothing) ? mz_group.Format : Float64 + default_intensity_format = (int_group !== nothing) ? int_group.Format : Float64 + mz_is_compressed = (mz_group !== nothing) ? mz_group.Packed : false + int_is_compressed = (int_group !== nothing) ? int_group.Packed : false + global_mode = (mz_group !== nothing && mz_group.Mode != UNKNOWN) ? mz_group.Mode : UNKNOWN + + # Use the first handle for metadata parsing (sequential) + # --- Metadata Caching Strategy (Sprint 1) --- + cache_path = file_path * ".cache" + use_cache = isfile(cache_path) && (mtime(cache_path) > mtime(file_path)) + + local spectra_metadata + if use_cache + @debug "Found valid metadata cache at $cache_path. Loading..." + try + spectra_metadata = load_metadata_cache(cache_path, default_mz_format, default_intensity_format) + @debug "Metadata loaded from cache in O(1) time." + catch e + @warn "Failed to load cache: $e. Falling back to full XML parsing." + use_cache = false end end - if mz_group === nothing || int_group === nothing - @warn "Could not find global definitions for m/z and intensity arrays. Using hardcoded defaults (Float64)." - default_mz_format = Float64 - default_intensity_format = Float64 - mz_is_compressed = false - int_is_compressed = false - global_mode = UNKNOWN - else - default_mz_format = mz_group.Format - default_intensity_format = int_group.Format - mz_is_compressed = mz_group.Packed - int_is_compressed = int_group.Packed - global_mode = mz_group.Mode != UNKNOWN ? mz_group.Mode : int_group.Mode + # Check for compression status once + any_comp = mz_is_compressed || int_is_compressed + + if !use_cache + @debug "Parsing spectrum block from XML (this may take time for large files)..." + spectra_metadata = parse_imzml_spectrum_block(stream, ibd_handles[1], param_groups, width, height, num_spectra, + default_mz_format, default_intensity_format, + mz_is_compressed, int_is_compressed, global_mode) + + @debug "Metadata parsing complete. Saving cache for next time..." + # We create a temporary MSIData just for save_metadata_cache + tmp_source = ImzMLSource(ibd_handles, default_mz_format, default_intensity_format, mmap_data, any_comp) + tmp_msi = MSIData(tmp_source, spectra_metadata, instrument_meta, (width, height), nothing, cache_size) + save_metadata_cache(tmp_msi, cache_path) end - - println("DEBUG: m/z format: $default_mz_format, Intensity format: $default_intensity_format") - println("DEBUG: m/z compressed: $mz_is_compressed, Intensity compressed: $int_is_compressed") - println("DEBUG: Global mode: $global_mode") - - local spectra_metadata = parse_imzml_spectrum_block(stream, ts_hIbd, param_groups, width, height, num_spectra, - default_mz_format, default_intensity_format, - mz_is_compressed, int_is_compressed, global_mode) - - println("DEBUG: Metadata parsing complete.") - # Build coordinate map for imzML files - println("DEBUG: Building coordinate map...") + # Build coordinate map ... coordinate_map = zeros(Int, width, height) for (idx, meta) in enumerate(spectra_metadata) - if idx == 1 - println("DIAGNOSTIC_WRITE: For index 1, attempting to write to coordinate_map[$(meta.x), $(meta.y)]") - end if 1 <= meta.x <= width && 1 <= meta.y <= height coordinate_map[meta.x, meta.y] = idx end end - println("DEBUG: Coordinate map built.") - # --- NEW: Update acquisition mode based on spectrum parsing --- - acq_mode_symbol = if global_mode == CENTROID - :centroid - elseif global_mode == PROFILE - :profile - else - :unknown - end - - final_instrument_meta = InstrumentMetadata( - instrument_meta.resolution, - acq_mode_symbol, # Update with parsed mode - instrument_meta.mz_axis_type, - instrument_meta.calibration_status, - instrument_meta.instrument_model, - instrument_meta.mass_accuracy_ppm, - instrument_meta.laser_settings, - instrument_meta.polarity, - instrument_meta.vendor_preprocessing_steps # Add this new field - ) + source = ImzMLSource(ibd_handles, default_mz_format, default_intensity_format, mmap_data, any_comp) + @debug "Creating MSIData object." + msi_data = MSIData(source, spectra_metadata, instrument_meta, (width, height), coordinate_map, cache_size) - source = ImzMLSource(ts_hIbd, default_mz_format, default_intensity_format) - println("DEBUG: Creating MSIData object.") - msi_data = MSIData(source, spectra_metadata, final_instrument_meta, (width, height), coordinate_map, cache_size) + # NOTE: Do NOT set analytics_ready here even though the cache provides fast metadata loading. + # The cache only stores binary offsets (SpectrumMetadataBinary). It does NOT populate + # msi_data.spectrum_stats_df (TIC, BPI, BasePeakMZ, MinMZ, MaxMZ), which requires a + # streaming pass via precompute_analytics(). Setting the flag prematurely causes + # get_mz_slice to skip that pass, leaving stats_df=nothing and all min/max bounds at 0.0, + # resulting in zero pixels populated in every image slice. + @debug "Metadata loaded from cache — analytics scan deferred until first use." - # Close the XML stream as it's no longer needed close(stream) - return msi_data catch e close(stream) - close(ts_hIbd) # Ensure IBD handle is closed on error + # Check if handles exist before closing + if @isdefined(ibd_handles) + for h in ibd_handles + isopen(h) && close(h) + end + end rethrow(e) end end @@ -820,8 +828,8 @@ function parse_compressed(stream::IO, hIbd::Union{IO, ThreadSafeFileHandle}, par if length(mz_data) != 1 || length(int_data) != 1 println("DEBUG: Spectrum $k is empty or invalid - creating placeholder metadata") - mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, Int64(0), 0, :mz) - int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, Int64(0), 0, :intensity) + mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, Int64(0), 0, :mz, 0.0, 0.0) + int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, Int64(0), 0, :intensity, 0.0, 0.0) else mz_info = mz_data[1] int_info = int_data[1] @@ -835,9 +843,9 @@ function parse_compressed(stream::IO, hIbd::Union{IO, ThreadSafeFileHandle}, par end mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, mz_info.offset, - mz_is_compressed ? mz_info.encoded_length : mz_info.array_length, :mz) + mz_is_compressed ? mz_info.encoded_length : mz_info.array_length, :mz, 0.0, 0.0) int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, int_info.offset, - int_is_compressed ? int_info.encoded_length : int_info.array_length, :intensity) + int_is_compressed ? int_info.encoded_length : int_info.array_length, :intensity, 0.0, 0.0) end spectra_metadata[k] = SpectrumMetadata(x, y, "", :sample, spectrum_mode, mz_asset, int_asset) @@ -869,7 +877,7 @@ This optimized version uses binary search for efficiency. # Returns - The intensity (`Float64`) of the peak if found, otherwise `0.0`. """ -function find_mass(mz_array::AbstractVector{<:Real}, intensity_array::AbstractVector{<:Real}, +@inline function find_mass(mz_array::AbstractVector{<:Real}, intensity_array::AbstractVector{<:Real}, target_mass::Real, tolerance::Real) # Fast-path rejection: if the array is empty or the target is out of range if isempty(mz_array) || target_mass + tolerance < first(mz_array) || target_mass - tolerance > last(mz_array) @@ -930,59 +938,91 @@ function get_mz_slice(data::MSIData, mass::Real, tolerance::Real; mask_path::Uni precompute_analytics(data) end - println("Using high-performance sequential iterator...") target_min = mass - tolerance target_max = mass + tolerance + + # PERFORMANCE: Access raw vectors from metadata to avoid DataFrame dependency + # If stats_df exists, use it; otherwise, use the pre-computed bounds in SpectrumAsset stats_df = get_spectrum_stats(data) - bloom_filters = get_bloom_filters(data) + + n_total = length(data.spectra_metadata) + # Instantiate thread-local buffer locally since global pools are deprecated + candidate_indices = Vector{Int}(undef, n_total) + + # We'll use local views of min/max if stats_df is missing to represent zero-allocation fallback + # But for extreme performance, we avoid list comprehensions [m... for m in ...] as they allocate. + local min_mzs::Vector{Float64} + local max_mzs::Vector{Float64} - # 1. Find all candidate spectra first for efficient filtering - candidate_indices = Set{Int}() - indices_to_check = masked_indices === nothing ? (1:length(data.spectra_metadata)) : masked_indices + if stats_df !== nothing + min_mzs = stats_df.MinMZ + max_mzs = stats_df.MaxMZ + else + # Fallback path: extract to local buffers or use metadata directly in loop + # For now, let's assume stats_df is usually populated by precompute_analytics. + # If not, we'll access it directly inside the filter loop. + end + candidate_count = 0 + indices_to_check = masked_indices === nothing ? (1:n_total) : masked_indices + discretization_factor = 100.0 + bloom_filters = get_bloom_filters(data) for i in indices_to_check - # NEW: Bloom filter check with discretization - if bloom_filters !== nothing && !is_empty(bloom_filters[i]) - discretization_factor = 100.0 - min_mass_int = round(Int, (mass - tolerance) * discretization_factor) - max_mass_int = round(Int, (mass + tolerance) * discretization_factor) - - found = false - for mass_int in min_mass_int:max_mass_int - if mass_int in bloom_filters[i] - found = true - break + # Range check first (cheapest) + # Access metadata directly if stats_df is missing to ensure zero-allocation + @inbounds meta = data.spectra_metadata[i] + s_min, s_max = (stats_df !== nothing) ? (min_mzs[i], max_mzs[i]) : (meta.mz_asset.min_val, meta.mz_asset.max_val) + + if target_max < s_min || target_min > s_max + continue + end + + # Bloom filter rejection (very fast) + if bloom_filters !== nothing + bf = bloom_filters[i] + if !is_empty(bf) + min_mass_int = round(Int, (mass - tolerance) * discretization_factor) + max_mass_int = round(Int, (mass + tolerance) * discretization_factor) + + found = false + @inbounds for mass_int in min_mass_int:max_mass_int + if mass_int in bf + found = true + break + end end - end - - if !found - continue # Definitely not in this spectrum + !found && continue end end - spec_min_mz = stats_df.MinMZ[i] - spec_max_mz = stats_df.MaxMZ[i] - if target_max >= spec_min_mz && target_min <= spec_max_mz - push!(candidate_indices, i) - end + candidate_count += 1 + @inbounds candidate_indices[candidate_count] = i end - println("Found $(length(candidate_indices)) candidate spectra (filtered from $(length(indices_to_check)) initial spectra)") - + # Use a view of the pre-allocated vector to avoid collect() allocations + valid_candidates = view(candidate_indices, 1:candidate_count) + # 2. Iterate using the optimized, low-allocation iterator - results_count = 0 - _iterate_spectra_fast(data, collect(candidate_indices)) do idx, mz_array, intensity_array - meta = data.spectra_metadata[idx] - intensity = find_mass(mz_array, intensity_array, mass, tolerance) - if intensity > 0.0 - if 1 <= meta.x <= width && 1 <= meta.y <= height - slice_matrix[meta.y, meta.x] = intensity - results_count += 1 + # Use Atomic for thread-safe increment and avoid Ref-boxing + results_count = Base.Threads.Atomic{Int}(0) + + # Use let block to ensure closure captures are optimized (avoid boxing) + let slice_matrix=slice_matrix, results_count=results_count, width=width, height=height, + spectra_metadata=data.spectra_metadata, mass=mass, tolerance=tolerance + + _iterate_spectra_fast(data, valid_candidates) do idx, mz_array, intensity_array + @inbounds meta = spectra_metadata[idx] + intensity = find_mass(mz_array, intensity_array, mass, tolerance) + if intensity > 0.0 + if 1 <= meta.x <= width && 1 <= meta.y <= height + @inbounds slice_matrix[meta.y, meta.x] = intensity + Base.Threads.atomic_add!(results_count, 1) + end end end end - println("Populated $results_count pixels with intensity data") + println("Populated $(results_count[]) pixels with intensity data") replace!(slice_matrix, NaN => 0.0) return slice_matrix end @@ -1007,13 +1047,14 @@ This is a highly performant function that iterates through the full dataset only function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, tolerance::Real; mask_path::Union{String, Nothing}=nothing) width, height = data.image_dims - # Sort masses to improve cache locality during search + # Sort masses to improve cache locality and allow binary search sorted_masses = sort(masses) + n_masses = length(sorted_masses) - # 1. Initialize a dictionary to hold the output slice matrices - slice_dict = Dict{Real, Matrix{Float64}}() + # 1. Initialize a dictionary to hold the output slice matrices (using Float32 for 50% RAM savings) + slice_dict = Dict{Real, Matrix{Float32}}() for mass in sorted_masses - slice_dict[mass] = zeros(Float64, height, width) + slice_dict[mass] = zeros(Float32, height, width) end local masked_indices::Union{Set{Int}, Nothing} = nothing @@ -1029,42 +1070,50 @@ function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, t precompute_analytics(data) end - println("Filtering candidate spectra for $(length(masses)) m/z values...") + println("Filtering candidate spectra for $n_masses m/z values...") stats_df = get_spectrum_stats(data) bloom_filters = get_bloom_filters(data) - candidate_indices = Set{Int}() + + # Use a BitSet for faster index tracking + candidate_indices = BitSet() indices_to_check = masked_indices === nothing ? (1:length(data.spectra_metadata)) : masked_indices - # 3. Find all spectra that could contain *any* of the requested masses. - for mass in sorted_masses - target_min = mass - tolerance - target_max = mass + tolerance - for i in indices_to_check - # If already a candidate, no need to check again - if i in candidate_indices - continue - end - # NEW: Bloom filter check with discretization + # 3. Optimized filtering: Iterate through spectra ONCE and check against all masses + # This changes complexity from O(M*N) to O(N * log M) or O(N + M) depending on range overlap + discretization_factor = 100.0 + for i in indices_to_check + spec_min = stats_df.MinMZ[i] + spec_max = stats_df.MaxMZ[i] + + # Binary search to find masses that might overlap with this spectrum's range + # target_min = mass - tolerance => mass = target_min + tolerance + # We need mass such that mass + tolerance >= spec_min => mass >= spec_min - tolerance + # and mass - tolerance <= spec_max => mass <= spec_max + tolerance + + m_start_idx = searchsortedfirst(sorted_masses, spec_min - tolerance) + m_end_idx = searchsortedlast(sorted_masses, spec_max + tolerance) + + if m_start_idx <= m_end_idx + # Range overlap found, now check Bloom filter if available if bloom_filters !== nothing && !is_empty(bloom_filters[i]) - discretization_factor = 100.0 - min_mass_int = round(Int, (mass - tolerance) * discretization_factor) - max_mass_int = round(Int, (mass + tolerance) * discretization_factor) - - found = false - for mass_int in min_mass_int:max_mass_int - if mass_int in bloom_filters[i] - found = true - break + found_any = false + @inbounds for m_idx in m_start_idx:m_end_idx + mass = sorted_masses[m_idx] + min_mass_int = round(Int, (mass - tolerance) * discretization_factor) + max_mass_int = round(Int, (mass + tolerance) * discretization_factor) + + for mass_int in min_mass_int:max_mass_int + if mass_int in bloom_filters[i] + found_any = true + break + end end + found_any && break end - - if !found - continue # Definitely not in this spectrum + if found_any + push!(candidate_indices, i) end - end - spec_min_mz = stats_df.MinMZ[i] - spec_max_mz = stats_df.MaxMZ[i] - if target_max >= spec_min_mz && target_min <= spec_max_mz + else push!(candidate_indices, i) end end @@ -1073,29 +1122,40 @@ function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, t println("Found $(length(candidate_indices)) total candidate spectra.") # 4. Iterate through the data a single time using the optimized iterator. + # We collect candidate_indices to pass to parallel iterator _iterate_spectra_fast(data, collect(candidate_indices)) do idx, mz_array, intensity_array meta = data.spectra_metadata[idx] - # For this single spectrum, check all masses of interest - for mass in sorted_masses - # Check if this spectrum's range actually covers the current mass - # This is a finer-grained check than the initial filtering - if !isempty(mz_array) && (mass + tolerance) >= first(mz_array) && (mass - tolerance) <= last(mz_array) - intensity = find_mass(mz_array, intensity_array, mass, tolerance) - if intensity > 0.0 - if 1 <= meta.x <= width && 1 <= meta.y <= height - slice_dict[mass][meta.y, meta.x] = intensity - end + if isempty(mz_array) + return + end + + # Spectrum-level boundaries + spec_first = first(mz_array) + spec_last = last(mz_array) + + # Find which of our target masses fall within this specific spectrum's actual range + m_start_idx = searchsortedfirst(sorted_masses, spec_first - tolerance) + m_end_idx = searchsortedlast(sorted_masses, spec_last + tolerance) + + @inbounds for m_idx in m_start_idx:m_end_idx + mass = sorted_masses[m_idx] + intensity = find_mass(mz_array, intensity_array, mass, tolerance) + if intensity > 0.0 + if 1 <= meta.x <= width && 1 <= meta.y <= height + # Note: Concurrent writes to different matrices/coordinates are safe. + # Dictionary access is safe because it's read-only after initialization. + slice_dict[mass][meta.y, meta.x] = Float32(intensity) end end end end - # 5. Clean up and return + # 5. Clean up - replaces NaNs with 0.0 directly in Float32 matrices for mass in sorted_masses - replace!(slice_dict[mass], NaN => 0.0) + replace!(slice_dict[mass], NaN32 => 0.0f0) end - println("Finished generating $(length(masses)) slices in a single pass.") + println("Finished generating $n_masses slices in a single pass.") return slice_dict end @@ -1733,7 +1793,7 @@ Generates a colorbar image for a given slice of data. - `fig::Figure`: A figure with the colorbar. """ function generate_colorbar_image(slice_data::AbstractMatrix, color_levels::Int, output_path::String, - bounds::Tuple{Float64, Float64}; + bounds::Tuple{Real, Real}; use_triq::Bool=false, triq_prob::Float64=0.98, mask_path::Union{String, Nothing}=nothing) # Use the provided bounds instead of recalculating diff --git a/src/mzML.jl b/src/mzML.jl index 6864b2a..f9904b8 100644 --- a/src/mzML.jl +++ b/src/mzML.jl @@ -231,7 +231,7 @@ function get_spectrum_asset_metadata(stream::IO) #println("DEBUG: Exiting get_spectrum_asset_metadata.") # Create SpectrumAsset directly from the variables - return SpectrumAsset(data_format, compression_flag, binary_offset, encoded_length, axis) + return SpectrumAsset(data_format, compression_flag, binary_offset, encoded_length, axis, 0.0, 0.0) end # This function is updated to return the generic SpectrumMetadata struct @@ -394,38 +394,52 @@ then parses the metadata for each spectrum without loading the binary data. """ function load_mzml_lazy(file_path::String; cache_size::Int=100) println("DEBUG: Opening file stream for $file_path") - ts_stream = ThreadSafeFileHandle(file_path, "r") + + # --- Handle Pool Optimization --- + # Open multiple handles to the .mzML file to avoid lock contention + num_handles = Threads.nthreads() + mzml_handles = [open(file_path, "r") for _ in 1:num_handles] + + # Use the first handle for initial parsing + primary_handle = mzml_handles[1] try # --- NEW: Parse instrument metadata from header --- println("DEBUG: Parsing instrument metadata from header...") - instrument_meta = parse_instrument_metadata_mzml(ts_stream.handle) + instrument_meta = parse_instrument_metadata_mzml(primary_handle) - println("--- Extracted Instrument Metadata ---") - println("Resolution: ", instrument_meta.resolution) - println("Acquisition Mode (pre-check): ", instrument_meta.acquisition_mode) - println("Calibration Status: ", instrument_meta.calibration_status) - println("Instrument Model: ", instrument_meta.instrument_model) - println("Mass Accuracy (ppm): ", instrument_meta.mass_accuracy_ppm) - println("Laser Settings: ", instrument_meta.laser_settings) - println("Polarity: ", instrument_meta.polarity) - println("------------------------------------") - - seekstart(ts_stream.handle) # Reset stream after header parsing + seekstart(primary_handle) # Reset stream after header parsing println("DEBUG: Finding index offset...") - index_offset = find_index_offset(ts_stream.handle) + index_offset = find_index_offset(primary_handle) + + # --- NEW: Mmap Optimization with RAM Safety --- + mmap_data = nothing + try + file_size = filesize(file_path) + free_ram = Sys.free_memory() + if file_size > free_ram * 0.8 + @warn "Dataset size ($(round(file_size/1e9, digits=2)) GB) exceeds 80% of free RAM. Mmap will still work via 'Streaming' mode." + end + + println("DEBUG: Memory mapping .mzML file...") + seekstart(primary_handle) # Anchor Mmap to the beginning of the file to prevent overflow + mmap_data = Mmap.mmap(primary_handle, Vector{UInt8}, (file_size,)) + println("DEBUG: .mzML file mmapped successfully.") + catch e + @warn "Memory mapping failed for mzML, falling back to standard I/O: $e" + end + println("DEBUG: Seeking to index list at offset $index_offset.") - seek(ts_stream.handle, index_offset) + seek(primary_handle, index_offset) println("DEBUG: Searching for ''.") - if find_tag(ts_stream.handle, r" m == CENTROID, modes) - num_profile = count(m -> m == PROFILE, modes) + # Determine overall acquisition mode ... + num_centroid = count(m -> m.mode == CENTROID, spectra_metadata) + num_profile = count(m -> m.mode == PROFILE, spectra_metadata) acq_mode_symbol = if num_centroid > 0 && num_profile == 0 :centroid @@ -468,26 +476,28 @@ function load_mzml_lazy(file_path::String; cache_size::Int=100) else :unknown end - println("DEBUG: Inferred overall acquisition mode: $acq_mode_symbol (Centroid: $num_centroid, Profile: $num_profile)") final_instrument_meta = InstrumentMetadata( instrument_meta.resolution, - acq_mode_symbol, # Update with parsed mode + acq_mode_symbol, instrument_meta.mz_axis_type, instrument_meta.calibration_status, instrument_meta.instrument_model, instrument_meta.mass_accuracy_ppm, instrument_meta.laser_settings, instrument_meta.polarity, - instrument_meta.vendor_preprocessing_steps # Add this new field + instrument_meta.vendor_preprocessing_steps ) - source = MzMLSource(ts_stream, mz_format, intensity_format) + source = MzMLSource(mzml_handles, mz_format, intensity_format, mmap_data) println("DEBUG: Creating MSIData object.") return MSIData(source, spectra_metadata, final_instrument_meta, (0, 0), nothing, cache_size) catch e - close(ts_stream) # Ensure stream is closed on error + # Close all handles in the pool if initialization fails + for h in mzml_handles + isopen(h) && close(h) + end rethrow(e) end end diff --git a/start_MSI_GUI.jl b/start_MSI_GUI.jl index d523fc8..5bf4bae 100644 --- a/start_MSI_GUI.jl +++ b/start_MSI_GUI.jl @@ -19,6 +19,52 @@ end using Genie +# --- Cross-Platform Startup Cleanup --- +# Remove orphaned GenieSessionFileSession directories from previous runs. +# These accumulate in the OS temp directory as jl_XXXXXX folders containing +# serialized session files (64-char hex filenames). Over long sessions or +# after crashes, they can consume gigabytes of disk space. +function cleanup_orphaned_sessions() + tmp = Base.tempdir() + cleaned_count = 0 + cleaned_bytes = 0 + + for entry in readdir(tmp; join=false) + # Only target directories matching Julia's temp naming pattern + startswith(entry, "jl_") || continue + full_path = joinpath(tmp, entry) + isdir(full_path) || continue + + # Validate: a Genie session dir contains files with 64-char hex names + try + contents = readdir(full_path) + isempty(contents) && continue + + # Check if at least one file matches the 64-char hex session ID pattern + is_session_dir = any(contents) do f + length(f) == 64 && all(c -> c in "0123456789abcdef", f) + end + is_session_dir || continue + + # Safe to remove — this is an orphaned Genie session directory + dir_size = sum(filesize(joinpath(full_path, f)) for f in contents; init=0) + rm(full_path; recursive=true, force=true) + cleaned_count += 1 + cleaned_bytes += dir_size + catch e + @debug "Skipping $entry during cleanup: $e" + end + end + + if cleaned_count > 0 + size_mb = round(cleaned_bytes / (1024^2), digits=1) + @info "Startup cleanup: removed $cleaned_count orphaned session dir(s), freed $(size_mb) MB" + end +end + +cleanup_orphaned_sessions() + + # Load and configure Genie Genie.loadapp() diff --git a/test/new_benchmark_mmap.jl b/test/new_benchmark_mmap.jl new file mode 100644 index 0000000..53aaca7 --- /dev/null +++ b/test/new_benchmark_mmap.jl @@ -0,0 +1,48 @@ +using BenchmarkTools +using MSI_src +using Statistics +using DataFrames + +# =================================================================== +# HIGH-PRECISION COMPARATIVE SUITE +# =================================================================== + +function run_advanced_benchmark(path, mz, tol) + println("\n" * "="^40) + println("TARGET: $(basename(path))") + println("="^40) + + # 1. NEW LIBRARY: Metadata Load (The "Control Tower" startup) + # This measures how fast the Mmap and Cache system works + t_load_new = @belapsed OpenMSIData($path) + + # 2. NEW LIBRARY: Slice Generation (The "Streaming" speed) + msi_new = OpenMSIData(path) + # We use @benchmark to get a distribution (min, mean, max) + b_slice_new = @benchmark get_mz_slice($msi_new, $mz, $tol) + + # --- Metrics Table --- + results = DataFrame( + Metric = ["Metadata Load", "Slice Gen (Min)", "Slice Gen (Mean)", "Allocations"], + JuliaMSI = [ + "$(round(t_load_new * 1000, digits=2)) ms", + "$(round(minimum(b_slice_new.times)/1e6, digits=2)) ms", + "$(round(mean(b_slice_new.times)/1e6, digits=2)) ms", + "$(b_slice_new.allocs) allocs" + ] + ) + + println(results) + + # --- The "Throughput" Test --- + # How many slices per second can we handle? + throughput_new = 1.0 / mean(b_slice_new.times/1e9) + println("\nThroughput: $(round(throughput_new, digits=1)) slices/sec") + + return results +end + +# Example Run +@time run_advanced_benchmark("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML", 716.053, 0.1) + +# For multithread: julia --threads auto --project=. test/new_benchmark_mmap.jl \ No newline at end of file diff --git a/test/run_precalculation_example.jl b/test/run_precalculation_example.jl index 05f7cec..53d1de9 100644 --- a/test/run_precalculation_example.jl +++ b/test/run_precalculation_example.jl @@ -16,8 +16,9 @@ using MSI_src const TEST_MZML_FILE = "" # const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.imzML" #const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML" -const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Thricoderma_etc/Imaging_interaccion_trichoderma_vs_streptomyces.imzML" +# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Thricoderma_etc/Imaging_interaccion_trichoderma_vs_streptomyces.imzML" #const MASK_ROUTE = "/home/pixel/Documents/Cinvestav_2025/JuliaMSI/public/css/masks/Stomach_DHB_uncompressed.png" +const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML" const MASK_ROUTE = "" #= diff --git a/test/run_preprocessing.jl b/test/run_preprocessing.jl index e244f8a..68354cf 100644 --- a/test/run_preprocessing.jl +++ b/test/run_preprocessing.jl @@ -19,11 +19,11 @@ using MSI_src # const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Thricoderma_etc/Imaging_interaccion_trichoderma_vs_streptomyces.imzML" # const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/set de datos MS/Atropina_tuneo_fraq_20ev.mzML" -const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML" -# const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML" +# const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML" +const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML" -const MASK_ROUTE = "/home/pixel/Documents/Cinvestav_2025/JuliaMSI/public/css/masks/Stomach_DHB_uncompressed.png" -# const MASK_ROUTE = "" +# const MASK_ROUTE = "/home/pixel/Documents/Cinvestav_2025/JuliaMSI/public/css/masks/Stomach_DHB_uncompressed.png" +const MASK_ROUTE = "" const OUTPUT_DIR = "./test/results/preprocessing_results" diff --git a/test/run_tests.jl b/test/run_tests.jl index 6de7a38..0504fc7 100644 --- a/test/run_tests.jl +++ b/test/run_tests.jl @@ -31,9 +31,10 @@ using MSI_src # --- Test Case 1: Standard .mzML file --- # A regular, non-imaging mzML file. # const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/mzML/T9_A1.mzML" -const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML" +# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML" # const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging_paper_spray/Imaging_paper_spray.mzML" # const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging prueba Roya 1/Roya.mzML" +const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/mzML" const SPECTRUM_TO_PLOT = 1 # Which spectrum to plot from the file # --- Test Case 2: .mzML + Sync File for Conversion --- @@ -58,16 +59,17 @@ const CONVERSION_TARGET_IMZML = "test/results/converted_mzml.imzML" # const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging prueba Roya 1/royaimg.imzML" # const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/ltpmsi-chilli.imzML" # centroid aparently? # const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_compressed.imzML" # centroid compressed -const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML" # centroid +# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML" # centroid +const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML" # The m/z value to use for creating an image slice. # const MZ_VALUE_FOR_SLICE = 309.06 # BF -# const MZ_VALUE_FOR_SLICE = 896.0 # HR2MSI +const MZ_VALUE_FOR_SLICE = 896.0 # HR2MSI # const MZ_VALUE_FOR_SLICE = 76.03 # I PS # const MZ_VALUE_FOR_SLICE = 313 # ROYA -const MZ_VALUE_FOR_SLICE = 100 # advanced processing +# const MZ_VALUE_FOR_SLICE = 100 # advanced processing # const MZ_TOLERANCE = 0.1 # const MZ_TOLERANCE = 1 -const MZ_TOLERANCE = 0.1 +const MZ_TOLERANCE = 0.2 # Coordinates to plot a specific spectrum from imzML const COORDS_TO_PLOT = (50, 50) # Example coordinates (X, Y) diff --git a/test/test_streaming_pipeline.jl b/test/test_streaming_pipeline.jl new file mode 100644 index 0000000..7e72cea --- /dev/null +++ b/test/test_streaming_pipeline.jl @@ -0,0 +1,124 @@ +#!/usr/bin/env julia +# test/test_streaming_pipeline.jl +# ============================================================================ +# Validation test for Sprint 2: The Streaming Pipeline +# +# This test exercises process_dataset! against the HR2MSI mouse bladder +# dataset and verifies: +# 1. Correct sparse matrix creation +# 2. Non-zero peak population +# 3. RAM savings vs dense equivalent +# 4. Allocation count and throughput +# ============================================================================ + +using Pkg +Pkg.activate(".") + +using MSI_src +using SparseArrays + +# ============================================================================= +# Configuration +# ============================================================================= + +const IMZML_PATH = "/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML" + +function main() + println("=" ^ 60) + println("SPRINT 2: Streaming Pipeline Validation") + println("=" ^ 60) + + if !isfile(IMZML_PATH) + println("SKIPPED: Dataset not found at $IMZML_PATH") + return + end + + # --- 1. Load dataset --- + println("\n--- Step 1: Loading dataset ---") + data = OpenMSIData(IMZML_PATH) + println("Loaded: $(length(data.spectra_metadata)) spectra") + + # --- 2. Configure the streaming pipeline --- + println("\n--- Step 2: Configuring pipeline ---") + config = PipelineConfig( + steps = [ + StreamingStep(:baseline_correction, Dict{Symbol,Any}(:method => :snip, :iterations => 50)), + StreamingStep(:normalization, Dict{Symbol,Any}(:method => :tic)), + StreamingStep(:peak_picking, Dict{Symbol,Any}( + :method => :profile, + :snr_threshold => 3.0, + :half_window => 10, + :min_peak_prominence => 0.1, + :merge_peaks_tolerance => 0.002 + )), + ], + num_bins = 2000, + frequency_threshold = 0.01 # Bins must appear in at least 1% of spectra + ) + println("Steps: $(join([s.name for s in config.steps], " → "))") + println("Bins: $(config.num_bins), Frequency threshold: $(config.frequency_threshold)") + + # --- 3. Run the streaming pipeline --- + println("\n--- Step 3: Running streaming pipeline ---") + stats = @timed begin + feature_matrix, bin_centers = process_dataset!(data, config) + end + + feature_matrix = stats.value[1] + bin_centers = stats.value[2] + + println("\n--- Results ---") + println(" Feature matrix size: $(size(feature_matrix))") + println(" Non-zeros: $(nnz(feature_matrix))") + println(" Bin centers: $(length(bin_centers))") + println(" Time: $(round(stats.time, digits=2))s") + println(" Allocations: $(stats.bytes ÷ 1_000_000) MB") + println(" GC time: $(round(stats.gctime, digits=2))s") + + # --- 4. Validate --- + println("\n--- Step 4: Validation ---") + + passed = true + + # Check matrix dimensions + if size(feature_matrix, 1) > 0 && size(feature_matrix, 2) > 0 + println(" ✓ Matrix has valid dimensions") + else + println(" ✗ Matrix has invalid dimensions: $(size(feature_matrix))") + passed = false + end + + # Check non-zeros + if nnz(feature_matrix) > 0 + println(" ✓ Matrix has $(nnz(feature_matrix)) non-zero entries") + else + println(" ✗ Matrix is completely empty") + passed = false + end + + # Check sparsity savings + dense_mb = size(feature_matrix, 1) * size(feature_matrix, 2) * 8 / 1e6 + sparse_mb = nnz(feature_matrix) * 16 / 1e6 # index + value per entry + if dense_mb > 0 + savings = (1.0 - sparse_mb / dense_mb) * 100 + println(" ✓ RAM savings: $(round(savings, digits=1))% ($(round(sparse_mb, digits=1)) MB vs $(round(dense_mb, digits=1)) MB dense)") + end + + # Check bin centers alignment + if length(bin_centers) == size(feature_matrix, 1) + println(" ✓ Bin centers match matrix rows") + else + println(" ✗ Bin center count ($(length(bin_centers))) != matrix rows ($(size(feature_matrix, 1)))") + passed = false + end + + println("\n" * "=" ^ 60) + if passed + println("ALL VALIDATIONS PASSED ✓") + else + println("SOME VALIDATIONS FAILED ✗") + end + println("=" ^ 60) +end + +@time main() -- 2.43.0 From 05cb79253fe469680cdc222b7edf2d988a818cea Mon Sep 17 00:00:00 2001 From: Pixelguy14 Date: Thu, 9 Apr 2026 13:46:30 -0600 Subject: [PATCH 5/6] cleaned the repository from dead code, simplified some functions, users can now create a system image of the JuliaMSI repo that runs the core dependencies faster without compiler time, alongside a docker environment for running scripts, the instructions are on the readme --- .github/workflows/ci.yml | 32 ++++ Dockerfile.headless | 31 ++++ Project.toml | 9 + README.md | 45 +++-- app.jl.html | 20 ++- build_sysimage.jl | 40 +++-- precompile_script.jl | 16 ++ src/BloomFilters.jl | 23 --- src/Common.jl | 76 +------- src/FusedPipeline.jl | 93 ---------- src/MSIData.jl | 370 +++++---------------------------------- src/MSI_src.jl | 27 ++- src/ResourcePool.jl | 16 -- src/imzML.jl | 253 +------------------------- test/runtests.jl | 58 ++++++ 15 files changed, 281 insertions(+), 828 deletions(-) create mode 100644 .github/workflows/ci.yml create mode 100644 Dockerfile.headless delete mode 100644 src/FusedPipeline.jl create mode 100644 test/runtests.jl diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..bebd45f --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,32 @@ +name: CI + +on: + push: + branches: + - main + pull_request: + +jobs: + test: + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + - '1.10' + os: + - ubuntu-latest + - windows-latest + - macOS-latest + arch: + - x64 + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: julia-actions/cache@v2 + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-runtest@v1 diff --git a/Dockerfile.headless b/Dockerfile.headless new file mode 100644 index 0000000..788e600 --- /dev/null +++ b/Dockerfile.headless @@ -0,0 +1,31 @@ +# Dockerfile.headless +# This Dockerfile creates a minimal, high-speed container solely intended +# for running JuliaMSI scripts via the pre-compiled headless sysimage. +# It explicitly excludes the Genie UI and graphical routing to optimize size and start speed. + +FROM julia:1.11 + +# System dependencies for core MSI math libraries +RUN apt-get update && apt-get install -y \ + build-essential \ + zlib1g-dev \ + && rm -rf /var/lib/apt/lists/* + +WORKDIR /app + +# Only copy essential library files (no UI or app.jl required) +COPY Project.toml Manifest.toml ./ +COPY src/ ./src/ +COPY build_sysimage.jl precompile_script.jl ./ +COPY test/ ./test/ + +# Instantiate the environment +RUN julia --project=. -e 'import Pkg; Pkg.instantiate(); Pkg.precompile()' + +# Build the headless sysimage +# This command strips out Genie and plot libraries, creating a pure data-engine .so +RUN julia --project=. build_sysimage.jl --headless + +# The resulting image is heavily optimized for zero-JIT execution. +# Default command simply enters a fast REPL. +CMD ["julia", "--project=.", "--threads", "auto", "-J", "sys_msi_headless.so"] diff --git a/Project.toml b/Project.toml index e171608..8225540 100644 --- a/Project.toml +++ b/Project.toml @@ -1,4 +1,5 @@ name = "MSI_src" +uuid = "8f395f85-3b9c-4f7f-bf7e-12345678abcd" authors = ["JJSA"] version = "0.1.0" @@ -47,6 +48,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StipplePlotly = "ec984513-233d-481d-95b0-a3b58b97af2b" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" [compat] @@ -92,5 +94,12 @@ Setfield = "1.1" Statistics = "1.11" StatsBase = "0.34" StipplePlotly = "0.13" +Test = "1.11.0" UUIDs = "1.11" julia = "1.11, 1.12" + +[extras] +Test = "8dfed614-e60c-5e00-aba5-a33c2d43236e" + +[targets] +test = ["Test"] diff --git a/README.md b/README.md index e83f4b9..517c8cd 100644 --- a/README.md +++ b/README.md @@ -35,35 +35,44 @@ Minimum system requirements: 4 core processor, 8 GB RAM
JuliaMSI is a Graphical User Interface for a library of MSI tools in Julia: https://github.com/CINVESTAV-LABI/julia_mzML_imzML ## Build the System Image (One-time) -Alternatively, you can generate the .so file by running the next build script in your directory: +Alternatively, you can generate custom pre-compiled `.so` / `.dll` system images by running the build script in your directory. This bakes the heavy mass-spectrometry kernels and UI dependencies into a single machine-code binary for ultra-fast startup times. +### 1. Build the GUI Sysimage (for full App usage) ```bash julia --project=. build_sysimage.jl ``` -This may take 5–10 minutes as it merges all dependencies and JuliaMSI functions into a single binary. The resulting .so/.dll file will be around **300MB–600MB** because it contains the pre-compiled machine code for your entire environment. -A sysimage built on Linux (.so) will not work on Windows. -You must run the build_sysimage.jl script once on each target operating system. +This may take 5–15 minutes. The resulting `sys_msi_gui.so` (or `.dll` on Windows) file will be around **300MB–600MB** because it contains the pre-compiled machine code for your entire graphical environment. +*Note: A sysimage built on Linux (.so) will not work on Windows. You must run the build script once on each target operating system.* -Once MSI_sysimage.so is created in your directory, adapt your command like this: +Once `sys_msi_gui.so` is created in your directory, adapt your launch command: ```bash julia --project=. -e 'using Pkg; Pkg.precompile()' -# This will find MSI_sysimage.so, .dll, or .dylib automatically: -julia --threads auto --project=. --sysimage MSI_sysimage* start_MSI_GUI.jl +# Adjust the .so extension to .dll if on Windows or .dylib if on macOS +julia --threads auto --project=. --sysimage sys_msi_gui.so start_MSI_GUI.jl ``` -The command above loads a pre-compiled version of your environment (note, you have to use either the sysimage command or the normal start command), speeding up the booting delay. Your GUI or scripts using JuliaMSI should start almost instantly and process files significantly faster than before. -Always run the build script **in the same project root** where you intend to run the app. -In scripts like test/run_tests.jl, you no longer need to manually include("../src/MSI_src.jl"). You can simply treat it like a globally installed package to load the precompiled binary: - -```julia -using MSI_src -# ... rest of your script -``` -### Run with the --sysimage flag -Launch your tests using the same flag to bypass all compilation overhead: +### 2. Build the Headless Sysimage (for Scripts & Data Scientists) +If you only want to run scripts or the core `MSI_src` engine without the overhead of the Genie web server and Plotly, you can build a stripped-down, high-speed system image: ```bash -julia --threads auto --project=. --sysimage MSI_sysimage.so test/run_tests.jl +julia --project=. build_sysimage.jl --headless +``` + +Launch your automated tests or custom processing scripts using the headless image to bypass virtually all compilation overhead: +```bash +julia --threads auto --project=. --sysimage sys_msi_headless.so test/test_streaming_pipeline.jl +``` + +### 3. Docker Option for High-Speed Processing +For maximum portability and execution speed without clutter, we provide a `Dockerfile.headless`. This option creates a minimal, ultra-fast container that completely omits `app.jl` and the web interface. + +It automatically builds and bundles the headless sysimage, giving researchers an instant-start engine ready for heavy-duty `.imzML` batch pipelines with zero JIT latency. + +To build and run the Docker image: +```bash +docker build -t juliamsi-headless -f Dockerfile.headless . +# Run a script mapped from your local data folder: +docker run -it -v /my/local/data:/data juliamsi-headless julia -J sys_msi_headless.so /data/my_processing_script.jl ``` ## License diff --git a/app.jl.html b/app.jl.html index 8dd5d49..907fc89 100644 --- a/app.jl.html +++ b/app.jl.html @@ -546,13 +546,19 @@ -
- - - +
+
+ Feature matrices are exported in highly optimized MatrixMarket Coordinate format (.mtx). +
Python: scipy.io.mmread('filtered_data.mtx') | R: Matrix::readMM('filtered_data.mtx') +
+
+ + + +
diff --git a/build_sysimage.jl b/build_sysimage.jl index f1c8e43..30d1a97 100644 --- a/build_sysimage.jl +++ b/build_sysimage.jl @@ -1,27 +1,35 @@ # build_sysimage.jl using PackageCompiler +headless = "--headless" in ARGS + println("--- Starting Custom System Image Build ---") -println("This process can take 5-10 minutes.") +println("Mode: ", headless ? "Headless (Scripting/Docker)" : "Full GUI (Genie App)") +println("This process can take 5-15 minutes.") -# Define the image path with OS-specific extension ext = Sys.iswindows() ? ".dll" : Sys.isapple() ? ".dylib" : ".so" -sysimage_path = "MSI_sysimage" * ext +sysimage_name = headless ? "sys_msi_headless" * ext : "sys_msi_gui" * ext +sysimage_path = joinpath(@__DIR__, sysimage_name) -# List dependencies to include in the sysimage for maximum stability -# We have removed graphical heavy-lifters (CairoMakie, PlotlyBase) to ensure -# the build completes successfully on systems with standard RAM. -packages_to_include = [ - :Genie, - :GenieFramework, - :PlotlyBase, - :Libz, +packages_to_include = Symbol[ + :MSI_src, # <--- Bake our core library :DataFrames, + :SparseArrays, + :Mmap, + :Libz, :Serialization, :Printf, :JSON ] +if !headless + append!(packages_to_include, [ + :Genie, + :GenieFramework, + :PlotlyBase + ]) +end + create_sysimage( packages_to_include, sysimage_path = sysimage_path, @@ -31,5 +39,11 @@ create_sysimage( ) println("--- System Image Build Successful! ---") -println("To start Julia with this image, use:") -println("julia --threads auto --project=. --sysimage $(sysimage_path) start_MSI_GUI.jl") +println("Created: $(sysimage_path)") +if headless + println("To use in a headles script:") + println("julia --threads auto --project=. --sysimage $(sysimage_name) my_script.jl") +else + println("To start Julia with the GUI image, use:") + println("julia --threads auto --project=. --sysimage $(sysimage_name) start_MSI_GUI.jl") +end diff --git a/precompile_script.jl b/precompile_script.jl index 050df9a..d1f7d6a 100644 --- a/precompile_script.jl +++ b/precompile_script.jl @@ -133,6 +133,22 @@ try imzml_data = MSI_src.load_imzml_lazy(imzml_path; use_mmap=true) MSI_src.get_multiple_mz_slices(imzml_data, [10.0, 20.0], 0.1) + # Exercise Sprint 2 Streaming Pipeline + config = MSI_src.PipelineConfig( + steps=[ + MSI_src.StreamingStep(:baseline_correction, Dict(:method => :snip, :iterations => 5)), + MSI_src.StreamingStep(:normalization, Dict(:method => :tic)), + MSI_src.StreamingStep(:peak_picking, Dict(:method => :profile, :snr_threshold => 3.0)) + ], + adaptive_binning=true, + bin_tolerance_ppm=50.0 + ) + try + MSI_src.process_dataset!(imzml_data, config) + catch + # Ignore matrix mismatch errors from tiny random data + end + # Exercise mzML pathway mzml_data = MSI_src.load_mzml_lazy(mzml_path) MSI_src.GetSpectrum(mzml_data, 1) diff --git a/src/BloomFilters.jl b/src/BloomFilters.jl index 3b555af..3570da7 100644 --- a/src/BloomFilters.jl +++ b/src/BloomFilters.jl @@ -20,30 +20,7 @@ mutable struct BloomFilter{T} count::Int end -""" - StreamingBloomFilter -A memory-efficient Bloom filter that doesn't store all values at once. - -# Fields -- `bits::BitVector`: The underlying bit array -- `size::Int`: Number of bits in the filter -- `hash_count::Int`: Number of hash functions to use -- `seed::UInt64`: Random seed for hash functions -- `count::Int`: Number of elements added (for monitoring) -- `current_chunk::Vector{Int}`: Current chunk of data being processed -- `chunk_size::Int`: Size of each chunk -""" -mutable struct StreamingBloomFilter - bits::BitVector - size::Int - hash_count::Int - seed::UInt64 - count::Int - # Streaming state - current_chunk::Vector{Int} - chunk_size::Int -end """ BloomFilter{T}(expected_elements::Int, false_positive_rate::Float64=0.01; kwargs...) -> Return type diff --git a/src/Common.jl b/src/Common.jl index de1139c..629aed2 100644 --- a/src/Common.jl +++ b/src/Common.jl @@ -31,76 +31,7 @@ function posix_madvise(buffer::AbstractArray, advice::Integer) end # --- Buffer Pooling --- -""" - BufferPool -A thread-safe pool of `Vector{UInt8}` buffers, categorized by their size. -This helps reduce memory allocations and garbage collection overhead by reusing buffers. - -# Arguments: -- `lock`: A `ReentrantLock` to ensure thread-safe access to the pool. -- `buffers`: A dictionary mapping buffer size (in bytes) to a list of available buffers of that size. -""" -mutable struct BufferPool - lock::ReentrantLock - buffers::Dict{Int, Vector{Vector{UInt8}}} - - function BufferPool() - new(ReentrantLock(), Dict{Int, Vector{Vector{UInt8}}}()) - end -end - -""" - get_buffer!(pool::BufferPool, size::Int) -> Vector{UInt8} - -Retrieves a `Vector{UInt8}` buffer of at least `size` bytes from the pool. -If no suitable buffer is available, a new one is allocated. - -# Arguments: -- `pool`: The `BufferPool` to retrieve the buffer from. -- `size`: The minimum desired size of the buffer in bytes. - -# Returns: -- A `Vector{UInt8}` buffer. -""" -function get_buffer!(pool::BufferPool, size::Int)::Vector{UInt8} - lock(pool.lock) do - # Find an existing buffer that is large enough - for (buffer_size, buffer_list) in pool.buffers - if buffer_size >= size && !isempty(buffer_list) - buffer = pop!(buffer_list) - # Resize if necessary (and if it's significantly larger than needed) - if length(buffer) > size * 2 # Heuristic: if buffer is more than twice the requested size - return Vector{UInt8}(undef, size) # Allocate new smaller buffer - else - return buffer - end - end - end - # No suitable buffer found, allocate a new one - return Vector{UInt8}(undef, size) - end -end - -""" - release_buffer!(pool::BufferPool, buffer::Vector{UInt8}) - -Returns a `Vector{UInt8}` buffer to the pool for future reuse. - -# Arguments: -- `pool`: The `BufferPool` to return the buffer to. -- `buffer`: The `Vector{UInt8}` buffer to release. -""" -function release_buffer!(pool::BufferPool, buffer::Vector{UInt8}) - lock(pool.lock) do - buffer_size = length(buffer) - if !haskey(pool.buffers, buffer_size) - pool.buffers[buffer_size] = Vector{Vector{UInt8}}() - end - push!(pool.buffers[buffer_size], buffer) - end - return nothing -end """ @@ -140,11 +71,16 @@ If no suitable buffer is available, a new one is allocated. - A `Vector{UInt8}` buffer. """ function get_buffer!(pool::SimpleBufferPool, size::Int)::Vector{UInt8} - lock(pool.lock) do + buf = lock(pool.lock) do # Check for existing buffers of exact size first if haskey(pool.buffers, size) && !isempty(pool.buffers[size]) return pop!(pool.buffers[size]) end + return nothing + end + + if buf !== nothing + return buf end # No suitable buffer found, allocate new one (outside lock to reduce contention) diff --git a/src/FusedPipeline.jl b/src/FusedPipeline.jl deleted file mode 100644 index 47f835a..0000000 --- a/src/FusedPipeline.jl +++ /dev/null @@ -1,93 +0,0 @@ -# src/FusedPipeline.jl -# This file defines a high-performance in-place preprocessing pipeline -# that minimizes allocations by using reused buffers from ResourcePool. - -using .MSI_src # Ensure it can see the module's contents if needed - -""" - SpectralPipeline - -Holds a sequence of preprocessing steps and the necessary buffers to execute them in-place. -""" -struct SpectralPipeline - steps::Vector{AbstractPreprocessingStep} -end - -""" - apply_pipeline!(mz::Vector{Float64}, intensity::Vector{Float64}, pipeline::SpectralPipeline; data::MSIData) - -Applies all steps in the pipeline to the spectrum arrays in-place. -Uses internal buffers from data.resource_pool where needed. -""" -function apply_pipeline!(mz::Vector{Float64}, intensity::Vector{Float64}, pipeline::SpectralPipeline, data::MSIData) - # Process each step in sequence - for step in pipeline.steps - apply_step!(mz, intensity, step, data) - end - return mz, intensity -end - -# --- Basic in-place implementations of core preprocessing steps --- - -function apply_step!(mz, int, step::Normalization, data) - if step.method === :tic - s = sum(int) - if s > 0 - int ./= s - end - elseif step.method === :median - m = median(int) - if m > 0 - int ./= m - end - end - return int -end - -function apply_step!(mz, int, step::Smoothing, data) - if step.method === :savitzky_golay - # SavitzkyGolay.savitzky_golay currently allocates, but we can't easily fix that here - # without refactoring the library. However, we can use a pooled vector for its output - # then copy back to intensity. - # [Wait: For now we'll call smoothed_y = smooth_spectrum_core(int, ...)] - # We'll use a resource from the pool to avoid fresh allocation - temp_buf = acquire(data.resource_pool) - resize!(temp_buf, length(int)) - - # Call existing core which returns a new vector, unfortunately - # But we'll copy it back to 'int' to maintain in-place pipeline - smoothed = smooth_spectrum_core(int; method=step.method, window=step.window, order=step.order) - copyto!(int, smoothed) - - release!(data.resource_pool, temp_buf) - end - return int -end - -function apply_step!(mz, int, step::BaselineCorrection, data) - if step.method === :snip - # SNIP is easy to make in-place! - iterations = (step.iterations === nothing) ? 100 : step.iterations - snip_baseline_inplace!(int, iterations) - end - return int -end - -""" - snip_baseline_inplace!(y, iterations) - -In-place implementation of the Sensitive Nonlinear Iterative Peak clipping algorithm. -""" -function snip_baseline_inplace!(y::Vector{Float64}, iterations::Int) - n = length(y) - n < 3 && return y - - # We still need one temporary buffer for the SNIP iteration to read from the previous state - # Actually, we can just return the baseline and subtract it, but to BE in-place, - # we need a temporary to hold the baseline during calculation. - - # For now, we'll use the existing _snip_baseline_impl and subtract - baseline = _snip_baseline_impl(y, iterations=iterations) - y .-= baseline - return y -end diff --git a/src/MSIData.jl b/src/MSIData.jl index b40a7ef..1e88fa3 100644 --- a/src/MSIData.jl +++ b/src/MSIData.jl @@ -10,100 +10,6 @@ using Base64, Libz, Serialization, Printf, DataFrames, Base.Threads, StatsBase, const FILE_HANDLE_LOCK = ReentrantLock() -""" - ThreadSafeFileHandle - -Wraps a file handle with thread-safe operations. -# Fields - -- `path::String`: The path to the file. -- `handle::IO`: The file handle. -- `lock::ReentrantLock`: The lock for thread-safe operations. -""" -mutable struct ThreadSafeFileHandle - path::String - handle::IO - lock::ReentrantLock -end - -""" - ThreadSafeFileHandle(path::String, mode::String="r") - -Wraps a file handle with thread-safe operations. - -# Arguments - -- `path::String`: The path to the file. -- `mode::String`: The mode to open the file in. - -# Returns - -- A `ThreadSafeFileHandle`. -""" -function ThreadSafeFileHandle(path::String, mode::String="r") - handle = lock(FILE_HANDLE_LOCK) do - open(path, mode) - end - return ThreadSafeFileHandle(path, handle, ReentrantLock()) -end - -""" - read_at!(tsfh::ThreadSafeFileHandle, a::AbstractArray, pos::Integer) - -Atomically seeks to a position and reads data into an array. This is the thread-safe -way to read from a specific offset in the file. -""" -function read_at!(tsfh::ThreadSafeFileHandle, a::AbstractArray, pos::Integer) - if pos < 0 - throw(ArgumentError("Invalid seek position: $pos")) - end - lock(tsfh.lock) do - # Also check against file size if possible, though filesize() might be expensive to call repeatedly - # Let's rely on seek throwing if it goes way out of bounds, but catch the negative case which is definitely an error. - seek(tsfh.handle, pos) - read!(tsfh.handle, a) - end -end - -function Base.read(tsfh::ThreadSafeFileHandle, n::Integer) - lock(tsfh.lock) do - return read(tsfh.handle, n) - end -end - -function Base.close(tsfh::ThreadSafeFileHandle) - close(tsfh.handle) -end - -function Base.isopen(tsfh::ThreadSafeFileHandle) - isopen(tsfh.handle) -end - -function Base.position(tsfh::ThreadSafeFileHandle) - lock(tsfh.lock) do - return position(tsfh.handle) - end -end - -function Base.filesize(tsfh::ThreadSafeFileHandle) - lock(tsfh.lock) do - return filesize(tsfh.handle) - end -end - -# Overload Base.seek for ThreadSafeFileHandle -function Base.seek(tsfh::ThreadSafeFileHandle, pos::Integer) - lock(tsfh.lock) do - seek(tsfh.handle, pos) - end -end - -# Overload Base.read! for ThreadSafeFileHandle -function Base.read!(tsfh::ThreadSafeFileHandle, a::AbstractArray) - lock(tsfh.lock) do - read!(tsfh.handle, a) - end -end # Abstract type for different data sources (e.g., mzML, imzML) # This allows dispatching to the correct binary reading logic. @@ -331,7 +237,7 @@ mutable struct MSIData global_min_mz::Base.Threads.Atomic{Float64} global_max_mz::Base.Threads.Atomic{Float64} spectrum_stats_df::Union{DataFrame, Nothing} - bloom_filters::Union{Vector{<:BloomFilter}, Nothing} + bloom_filters::Union{Vector{Union{BloomFilter{Int}, Nothing}}, Nothing} analytics_ready::AtomicFlag preprocessing_hints::Union{Dict{Symbol, Any}, Nothing} # For auto-determined parameters @@ -487,65 +393,7 @@ function Base.close(data::MSIData) println("MSIData object closed and resources released.") end -""" - warm_cache(data::MSIData, indices::AbstractVector{Int}) -Pre-loads a specified list of spectra into the cache. This is useful for warming up -the cache with frequently accessed spectra to improve performance for subsequent -interactive analysis. - -# Arguments - -- `data::MSIData`: The MSIData object. -- `indices::AbstractVector{Int}`: The indices of the spectra to warm the cache with. - -# Returns -- `nothing` -""" -function warm_cache(data::MSIData, indices::AbstractVector{Int}) - println("Warming cache with $(length(indices)) spectra...") - for idx in indices - # Calling GetSpectrum will load the data and place it in the cache - GetSpectrum(data, idx) - end - println("Cache warming complete.") -end - -""" - get_memory_usage(data::MSIData) - -Calculates and returns a summary of the memory currently used by the MSIData object, -including the spectrum cache and metadata. - -# Arguments - -- `data::MSIData`: The MSIData object. - -# Returns -- `(total_mb, cache_mb, metadata_mb, coordinate_map_mb)`: A named tuple with memory usage in MB for different components. -""" -function get_memory_usage(data::MSIData) - cache_bytes = 0 - lock(data.cache_lock) do - # This is an approximation; `sizeof` on a dictionary is not fully representative, - # but it's much better to sum the size of the actual vectors. - for (mz, intensity) in values(data.cache) - cache_bytes += sizeof(mz) + sizeof(intensity) - end - end - - metadata_bytes = sizeof(data.spectra_metadata) - coordinate_map_bytes = data.coordinate_map === nothing ? 0 : sizeof(data.coordinate_map) - - total_bytes = cache_bytes + metadata_bytes + coordinate_map_bytes - - return ( - total_mb = total_bytes / 1024^2, - cache_mb = cache_bytes / 1024^2, - metadata_mb = metadata_bytes / 1024^2, - coordinate_map_mb = coordinate_map_bytes / 1024^2 - ) -end """ get_masked_spectrum_indices(data::MSIData, mask_matrix::BitMatrix) -> Set{Int} @@ -648,26 +496,7 @@ function read_binary_vector(data::MSIData, io::IO, asset::SpectrumAsset) return out_array end -""" - read_binary_vector(data::MSIData, ts_handle::ThreadSafeFileHandle, asset::SpectrumAsset) -Reads a single spectrum's m/z and intensity arrays directly from the `.ibd` -binary file for an `.imzML` dataset. It handles reading the raw binary data -and converting it from little-endian to the host's native byte order. - -# Arguments -- `data`: The `MSIData` object. -- `ts_handle`: The `ThreadSafeFileHandle` containing the file handle and data formats. -- `asset`: The `SpectrumAsset` containing metadata for the array. - -# Returns -- A `Vector` of the appropriate type containing the decoded data. -""" -function read_binary_vector(data::MSIData, ts_handle::ThreadSafeFileHandle, asset::SpectrumAsset) - lock(ts_handle.lock) do - return read_binary_vector(data, ts_handle.handle, asset) - end -end # Overload for different source types """ @@ -747,9 +576,10 @@ This is a high-level function that orchestrates the reading process by calling - A tuple `(mz, intensity)` containing the two requested data arrays. """ function read_spectrum_from_disk(data::MSIData, source::MzMLSource, meta::SpectrumMetadata) + handle = get_handle(source) # For mzML, data is Base64 encoded within the XML - mz = read_binary_vector(data, source.file_handle, meta.mz_asset) - intensity = read_binary_vector(data, source.file_handle, meta.int_asset) + mz = read_binary_vector(data, handle, meta.mz_asset) + intensity = read_binary_vector(data, handle, meta.int_asset) mz_f64, intensity_f64 = Float64.(mz), Float64.(intensity) @@ -1064,8 +894,9 @@ function precompute_analytics(msi_data::MSIData) max_mzs = Vector{Float64}(undef, num_spectra) modes = Vector{SpectrumMode}(undef, num_spectra) - # Don't precompute Bloom filters by default - create them lazily - bloom_filters = nothing + # Precompute Bloom filters during this pass + bloom_filters = Vector{Union{BloomFilter{Int}, Nothing}}(undef, num_spectra) + fill!(bloom_filters, nothing) # Process in chunks to reduce memory pressure chunk_size = min(10000, num_spectra) @@ -1105,6 +936,7 @@ function precompute_analytics(msi_data::MSIData) bpis[idx] = max_int bp_mzs[idx] = mz[max_idx] num_points[idx] = length(mz) + bloom_filters[idx] = create_bloom_filter_for_spectrum(mz) end # Force GC between chunks to control memory @@ -1163,17 +995,17 @@ end + """ - get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000) + get_total_spectrum_core(msi_data::MSIData; num_bins::Int=2000, masked_indices::Union{Set{Int}, Nothing}=nothing) -Internal function to calculate the total spectrum for an `.imzML` file. +Internal function to calculate the total spectrum. -It uses a fast, two-pass approach: -1. The first pass finds the global m/z range across all spectra. +It uses a two-pass approach: +1. The first pass finds the global m/z range by iterating through all spectra. 2. The second pass sums intensities into a pre-defined number of bins. -This function is highly optimized for `.imzML` by leveraging direct binary -reading and optimized binning logic. It is called by `get_total_spectrum`. +This function is called by `get_total_spectrum`. # Arguments @@ -1182,16 +1014,10 @@ reading and optimized binning logic. It is called by `get_total_spectrum`. - `masked_indices::Union{Set{Int}, Nothing}`: The indices of the spectra to mask. # Returns -- `(mz_range, total_spectrum)`: A tuple containing the m/z range and the total spectrum. - -# Example - -```julia -get_total_spectrum_imzml(msi_data) -``` +- `(mz_range, total_spectrum, num_spectra_processed)`: A tuple containing the m/z range, total spectrum, and spectrum count. """ -function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_indices::Union{Set{Int}, Nothing}=nothing) - println("Calculating total spectrum for imzML (2-pass method)...") +function get_total_spectrum_core(msi_data::MSIData; num_bins::Int=2000, masked_indices::Union{Set{Int}, Nothing}=nothing) + println("Calculating total spectrum (2-pass method)...") total_start_time = time_ns() local global_min_mz, global_max_mz @@ -1201,7 +1027,7 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_ global_min_mz = Base.Threads.atomic_add!(msi_data.global_min_mz, 0.0) global_max_mz = Base.Threads.atomic_add!(msi_data.global_max_mz, 0.0) else - # 1. First Pass: Find the global m/z range by reading the fast .ibd file + # 1. First Pass: Find the global m/z range pass1_start_time = time_ns() println(" Pass 1: Finding global m/z range...") g_min_mz = Inf @@ -1214,8 +1040,9 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_ end end pass1_duration = (time_ns() - pass1_start_time) / 1e9 + if !isfinite(g_min_mz) - @warn "Could not determine a valid m/z range for imzML. All spectra might be empty." + @warn "Could not determine a valid m/z range. All spectra might be empty." return (Float64[], Float64[], 0) end @@ -1226,6 +1053,7 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_ println(" Global m/z range found and cached: [$(global_min_mz), $(global_max_mz)]") @printf " (Pass 1 took %.2f seconds)\n" pass1_duration end + # 2. Define Bins and precompute constants mz_bins = range(global_min_mz, stop=global_max_mz, length=num_bins) bin_step = step(mz_bins) @@ -1272,127 +1100,12 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_ pass2_duration = (time_ns() - pass2_start_time) / 1e9 total_duration = (time_ns() - total_start_time) / 1e9 - println("\n--- imzML Profiling (Post-Optimization) ---") + println("\n--- Profiling (Post-Optimization) ---") @printf " Pass 2 (I/O + Binning): %.2f seconds\n" pass2_duration @printf " Total Function Time: %.2f seconds\n" total_duration println("----------------------------------------\n") - println("Total spectrum calculation complete for imzML.") - return (collect(mz_bins), intensity_sum, num_spectra_processed) -end - -""" - get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000) - -Internal function to calculate the total spectrum for an `.mzML` file. - -It uses a two-pass approach analogous to the `imzML` implementation: -1. The first pass finds the global m/z range by iterating through all spectra. -2. The second pass sums intensities into a pre-defined number of bins. - -This function is called by `get_total_spectrum`. - -# Arguments - -- `msi_data::MSIData`: The MSIData object. -- `num_bins::Int`: The number of bins to use for the total spectrum. -- `masked_indices::Union{Set{Int}, Nothing}`: The indices of the spectra to mask. - -# Returns -- `(mz_range, total_spectrum)`: A tuple containing the m/z range and the total spectrum. - -# Example - -```julia -get_total_spectrum_mzml(msi_data) -``` -""" -function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000, masked_indices::Union{Set{Int}, Nothing}=nothing) - println("Calculating total spectrum for mzML (2-pass method)...") - total_start_time = time_ns() - - local global_min_mz, global_max_mz - - if Base.Threads.atomic_add!(msi_data.global_min_mz, 0.0) !== Inf - println(" Using pre-computed m/z range.") - global_min_mz = Base.Threads.atomic_add!(msi_data.global_min_mz, 0.0) - global_max_mz = Base.Threads.atomic_add!(msi_data.global_max_mz, 0.0) - else - # --- Pass 1: Find m/z range and cache it --- - pass1_start_time = time_ns() - println(" Pass 1: Finding global m/z range...") - g_min_mz = Inf - g_max_mz = -Inf - _iterate_spectra_fast(msi_data, masked_indices === nothing ? nothing : collect(masked_indices)) do idx, mz, intensity - if !isempty(mz) - local_min, local_max = extrema(mz) - g_min_mz = min(g_min_mz, local_min) - g_max_mz = max(g_max_mz, local_max) - end - end - pass1_duration = (time_ns() - pass1_start_time) / 1e9 - - if !isfinite(g_min_mz) - @warn "Could not determine a valid m/z range for mzML. All spectra might be empty." - return (Float64[], Float64[], 0) - end - - # Use and cache the result - global_min_mz = g_min_mz - global_max_mz = g_max_mz - set_global_mz_range!(msi_data, global_min_mz, global_max_mz) - - println(" Global m/z range found and cached: [$(global_min_mz), $(global_max_mz)]") - @printf " (Pass 1 took %.2f seconds)\n" pass1_duration - end - - # --- Pass 2: Bin intensities --- - pass2_start_time = time_ns() - println(" Pass 2: Summing intensities into $num_bins bins...") - - mz_bins = range(global_min_mz, stop=global_max_mz, length=num_bins) - bin_step = step(mz_bins) - inv_bin_step = 1.0 / bin_step # Precompute reciprocal - min_mz = global_min_mz - - # Initialize thread-local intensity sums - thread_local_intensity_sums = [zeros(Float64, num_bins) for _ in 1:Base.Threads.nthreads()] - thread_local_spectra_processed = [0 for _ in 1:Base.Threads.nthreads()] - - _iterate_spectra_fast(msi_data, masked_indices === nothing ? nothing : collect(masked_indices)) do idx, mz, intensity - thread_id = Base.Threads.threadid() - thread_local_spectra_processed[thread_id] += 1 - if isempty(mz) - return - end - - @inbounds @simd for i in eachindex(mz) - # Calculate raw bin index - raw_index = (mz[i] - min_mz) * inv_bin_step + 1.0 - bin_index = trunc(Int, raw_index) - - # Branchless version using clamp - final_index = clamp(bin_index, 1, num_bins) - thread_local_intensity_sums[thread_id][final_index] += intensity[i] - end - end - - # Combine thread-local results - intensity_sum = zeros(Float64, num_bins) - num_spectra_processed = 0 - for i in 1:Base.Threads.nthreads() - intensity_sum .+= thread_local_intensity_sums[i] - num_spectra_processed += thread_local_spectra_processed[i] - end - pass2_duration = (time_ns() - pass2_start_time) / 1e9 - - total_duration = (time_ns() - total_start_time) / 1e9 - println("\n--- mzML Profiling (Post-Optimization) ---") - @printf " Pass 2 (I/O + Binning): %.2f seconds\n" pass2_duration - @printf " Total Function Time: %.2f seconds\n" total_duration - println("----------------------------------------\n") - - println("Total spectrum calculation complete for mzML.") + println("Total spectrum calculation complete.") return (collect(mz_bins), intensity_sum, num_spectra_processed) end @@ -1428,11 +1141,7 @@ function get_total_spectrum(msi_data::MSIData; num_bins::Int=2000, mask_path::Un println("Calculating total spectrum for masked region (mask from: $(mask_path))") end - if msi_data.source isa ImzMLSource - return get_total_spectrum_imzml(msi_data, num_bins=num_bins, masked_indices=masked_indices) - else # MzMLSource - return get_total_spectrum_mzml(msi_data, num_bins=num_bins, masked_indices=masked_indices) - end + return get_total_spectrum_core(msi_data, num_bins=num_bins, masked_indices=masked_indices) end """ @@ -1726,8 +1435,9 @@ function _iterate_compressed_fast(f::Function, data::MSIData, source::ImzMLSourc # For compressed data, we currently allocate new arrays per spectrum. # To truly minimize allocations, we'd need read_compressed_array! (in-place version). # For now, let's ensure we are at least using the optimized type-stable version. - mz_array = read_compressed_array(data, source.ibd_handle, meta.mz_asset, source.mz_format) - intensity_array = read_compressed_array(data, source.ibd_handle, meta.int_asset, source.intensity_format) + handle = get_handle(source) + mz_array = read_compressed_array(data, handle, meta.mz_asset, source.mz_format) + intensity_array = read_compressed_array(data, handle, meta.int_asset, source.intensity_format) mz_array .= ltoh.(mz_array) intensity_array .= ltoh.(intensity_array) @@ -2044,19 +1754,19 @@ it will be created on-demand and cached. - A `BloomFilter{Int}` for the specified spectrum. """ function get_bloom_filter(data::MSIData, index::Int)::BloomFilter{Int} - lock(data.cache_lock) do - if data.bloom_filters === nothing - # Initialize Bloom filters on-demand - data.bloom_filters = Vector{Union{BloomFilter{Int}, Nothing}}(undef, length(data.spectra_metadata)) - fill!(data.bloom_filters, nothing) # Initialize with nothing + # Double-checked locking for resilience if not precomputed + if data.bloom_filters === nothing || data.bloom_filters[index] === nothing + lock(data.cache_lock) do + if data.bloom_filters === nothing + data.bloom_filters = Vector{Union{BloomFilter{Int}, Nothing}}(undef, length(data.spectra_metadata)) + fill!(data.bloom_filters, nothing) + end + + if data.bloom_filters[index] === nothing + mz, _ = GetSpectrum(data, index) + data.bloom_filters[index] = create_bloom_filter_for_spectrum(mz) + end end - - if data.bloom_filters[index] === nothing - # Create Bloom filter lazily - mz, intensity = GetSpectrum(data, index) - data.bloom_filters[index] = create_bloom_filter_for_spectrum(mz) - end - - return data.bloom_filters[index] end + return data.bloom_filters[index] end diff --git a/src/MSI_src.jl b/src/MSI_src.jl index aab9265..95d1206 100644 --- a/src/MSI_src.jl +++ b/src/MSI_src.jl @@ -21,7 +21,17 @@ export OpenMSIData, _iterate_spectra_fast, validate_spectrum, get_mz_slice, - REGISTRY_LOCK + REGISTRY_LOCK, + SpectrumMetadata, + SpectrumAsset, + SpectrumMode, + CENTROID, + PROFILE, + MzMLSource, + ImzMLSource, + SimpleBufferPool, + get_buffer!, + release_buffer! # Define shared registry lock const REGISTRY_LOCK = ReentrantLock() @@ -82,7 +92,6 @@ include("mzML.jl") include("imzML.jl") include("MzmlConverter.jl") include("Preprocessing.jl") -include("FusedPipeline.jl") include("ImageProcessing.jl") include("Precalculations.jl") include("PreprocessingPipeline.jl") @@ -111,13 +120,17 @@ This is the main entry point for the new data access API. - An `MSIData` object. """ function OpenMSIData(filepath::String; cache_size=300, spectrum_type_map::Union{Dict{Int, Symbol}, Nothing}=nothing) + # Apply standard path normalization for cross-platform compatibility + # Ensure Windows backslashes are converted to OS-native separators + norm_filepath = normpath(replace(filepath, "\\" => "/")) + local msi_data - if endswith(lowercase(filepath), ".mzml") - msi_data = load_mzml_lazy(filepath, cache_size=cache_size) - elseif endswith(lowercase(filepath), ".imzml") - msi_data = load_imzml_lazy(filepath, cache_size=cache_size) + if endswith(lowercase(norm_filepath), ".mzml") + msi_data = load_mzml_lazy(norm_filepath, cache_size=cache_size) + elseif endswith(lowercase(norm_filepath), ".imzml") + msi_data = load_imzml_lazy(norm_filepath, cache_size=cache_size) else - error("Unsupported file type: $filepath. Please provide a .mzML or .imzML file.") + error("Unsupported file type: $norm_filepath. Please provide a .mzML or .imzML file.") end # Apply spectrum type map if provided diff --git a/src/ResourcePool.jl b/src/ResourcePool.jl index c0e0413..aaa0ed8 100644 --- a/src/ResourcePool.jl +++ b/src/ResourcePool.jl @@ -21,23 +21,7 @@ mutable struct ResourcePool{T} constructor::Function end -""" - aligned_vector(::Type{T}, n::Int; alignment::Int=64) where T -Creates a `Vector{T}` that is aligned to `alignment` bytes. -Note: In modern Julia, standard vectors are often 16 or 64 byte aligned, but for -HPC we ensure this by allocating slightly more and using a view, or using -specific pointers. For simplicity and performance, we use a small hack: -allocating a larger array and taking a 64-byte aligned view. -""" -function aligned_vector(::Type{T}, n::Int; alignment::Int=64) where T - # Allocate enough space to find an aligned starting point - raw = Vector{UInt8}(undef, n * sizeof(T) + alignment) - ptr = Int(pointer(raw)) - off = (alignment - (ptr % alignment)) % alignment - # Return a reinterpret view of the aligned segment - return reinterpret(T, view(raw, (off + 1):(off + n * sizeof(T)))) -end """ ResourcePool{T}(constructor::Function; max_size::Int=2 * nthreads()) diff --git a/src/imzML.jl b/src/imzML.jl index 7b319ea..101f6fd 100644 --- a/src/imzML.jl +++ b/src/imzML.jl @@ -194,86 +194,6 @@ function axes_config_img(stream::IO) return param_groups end -""" - get_spectrum_tag_offset(stream) - -Calculates the character offset within a `` tag, ignoring attribute values. - -# Arguments - -- `stream::IO`: The stream to parse. - -# Returns - -- `offset::Int`: The character offset. -""" -function get_spectrum_tag_offset(stream::IO) - offset = position(stream) - tag = find_tag(stream, r"^\s*") - spectrum_end_skip = 0 - - return [x_skip, y_skip, mz_first, array_len_skip, spectrum_end_skip] -end - - """ read_spectrum_block(stream::IO) @@ -454,7 +374,7 @@ end # Arguments - `stream::IO`: The stream to parse. -- `hIbd::Union{IO, ThreadSafeFileHandle}`: The IBD file handle. +- `hIbd::IO`: The IBD file handle. - `param_groups::Dict{String, SpecDim}`: The parameter groups. - `width::Int32`: The width of the image. - `height::Int32`: The height of the image. @@ -469,7 +389,7 @@ end - `spectra_metadata::Vector{SpectrumMetadata}`: The spectrum metadata. """ -function parse_imzml_spectrum_block(stream::IO, hIbd::Union{IO, ThreadSafeFileHandle}, param_groups::Dict{String, SpecDim}, +function parse_imzml_spectrum_block(stream::IO, hIbd::IO, param_groups::Dict{String, SpecDim}, width::Int32, height::Int32, num_spectra::Int32, default_mz_format::DataType, default_intensity_format::DataType, mz_is_compressed::Bool, int_is_compressed::Bool, global_mode::SpectrumMode) @@ -685,175 +605,6 @@ function load_imzml_lazy(file_path::String; cache_size::Int=100, use_mmap::Bool= end end -""" - parse_compressed(stream::IO, hIbd::Union{IO, ThreadSafeFileHandle}, param_groups::Dict{String, SpecDim}, - width::Int32, height::Int32, num_spectra::Int32, - default_mz_format::DataType, default_intensity_format::DataType, - mz_is_compressed::Bool, int_is_compressed::Bool, global_mode::SpectrumMode) - -Helper function to parse spectrum metadata from an imzML file. - -# Arguments - -- `stream::IO`: The stream to parse. -- `hIbd::Union{IO, ThreadSafeFileHandle}`: The IBD file handle. -- `param_groups::Dict{String, SpecDim}`: The parameter groups. -- `width::Int32`: The width of the image. -- `height::Int32`: The height of the image. -- `num_spectra::Int32`: The number of spectra. -- `default_mz_format::DataType`: The default data type for mz arrays. -- `default_intensity_format::DataType`: The default data type for intensity arrays. -- `mz_is_compressed::Bool`: Whether mz arrays are compressed. -- `int_is_compressed::Bool`: Whether intensity arrays are compressed. -- `global_mode::SpectrumMode`: The global mode. - -# Returns - -- `spectra_metadata::Vector{SpectrumMetadata}`: The spectrum metadata. -""" -function parse_compressed(stream::IO, hIbd::Union{IO, ThreadSafeFileHandle}, param_groups::Dict{String, SpecDim}, - width::Int32, height::Int32, num_spectra::Int32, - default_mz_format::DataType, default_intensity_format::DataType, - mz_is_compressed::Bool, int_is_compressed::Bool, global_mode::SpectrumMode) - spectra_metadata = Vector{SpectrumMetadata}(undef, num_spectra) - - array_data_type = @NamedTuple{is_mz::Bool, array_length::Int32, encoded_length::Int64, offset::Int64} - - for k in 1:num_spectra - # Initialize variables for this spectrum - x = Int32(0) - y = Int32(0) - spectrum_mode = global_mode - current_array_data = array_data_type[] - spectrum_start_line = "" - - # Find the start of the spectrum tag - line = "" - while !eof(stream) - line = readline(stream) - if occursin("", line) - break - end - - # Parse coordinates - x_match = match(r"IMS:1000050.*?value=\"(\d+)\"", line) - if x_match !== nothing - x = parse(Int32, x_match.captures[1]) - end - y_match = match(r"IMS:1000051.*?value=\"(\d+)\"", line) - if y_match !== nothing - y = parse(Int32, y_match.captures[1]) - end - - # Parse mode - if occursin("MS:1000127", line) - spectrum_mode = CENTROID - elseif occursin("MS:1000128", line) - spectrum_mode = PROFILE - end - - # More robust binaryDataArray detection - if occursin("", line) - while !eof(stream) - bda_line = readline(stream) - push!(bda_lines, bda_line) - if occursin("", bda_line) - break - end - end - end - bda_content = join(bda_lines, "\n") - - # Parse from bda_content - is_mz = occursin("MS:1000514", bda_content) || occursin("mzArray", bda_content) - - array_len_cv_match = match(r"IMS:1000103.*?value=\"(\d+)\"", bda_content) - array_length = Int32(0) - if array_len_cv_match !== nothing - array_length = parse(Int32, array_len_cv_match.captures[1]) - end - - # Fallback for defaultArrayLength - if array_length == 0 - default_match = match(r"defaultArrayLength=\"(\d+)\"", spectrum_start_line) - if default_match !== nothing - array_length = parse(Int32, default_match.captures[1]) - end - end - - encoded_len_cv_match = match(r"IMS:1000104.*?value=\"(\d+)\"", bda_content) - encoded_length = Int64(0) - if encoded_len_cv_match !== nothing - encoded_length = parse(Int64, encoded_len_cv_match.captures[1]) - else - encoded_len_attr_match = match(r"encodedLength=\"(\d+)\"", bda_content) - if encoded_len_attr_match !== nothing - encoded_length = parse(Int64, encoded_len_attr_match.captures[1]) - end - end - - offset_match = match(r"IMS:1000102.*?value=\"(\d+)\"", bda_content) - offset = Int64(0) - if offset_match !== nothing - offset = parse(Int64, offset_match.captures[1]) - end - - if array_length > 0 && offset > 0 - push!(current_array_data, (is_mz=is_mz, array_length=array_length, - encoded_length=encoded_length, offset=offset)) - end - end - - # Reset line to continue loop - line = "" - end - - # Separate m/z and intensity arrays - mz_data = filter(d -> d.is_mz, current_array_data) - int_data = filter(d -> !d.is_mz, current_array_data) - - if length(mz_data) != 1 || length(int_data) != 1 - println("DEBUG: Spectrum $k is empty or invalid - creating placeholder metadata") - mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, Int64(0), 0, :mz, 0.0, 0.0) - int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, Int64(0), 0, :intensity, 0.0, 0.0) - else - mz_info = mz_data[1] - int_info = int_data[1] - - if k == 1 - println("DEBUG First spectrum parsed:") - println(" Coordinates: x=$x, y=$y") - println(" Mode: $spectrum_mode") - println(" m/z array: array_length=$(mz_info.array_length), encoded_length=$(mz_info.encoded_length), offset=$(mz_info.offset)") - println(" intensity array: array_length=$(int_info.array_length), encoded_length=$(int_info.encoded_length), offset=$(int_info.offset)") - end - - mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, mz_info.offset, - mz_is_compressed ? mz_info.encoded_length : mz_info.array_length, :mz, 0.0, 0.0) - int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, int_info.offset, - int_is_compressed ? int_info.encoded_length : int_info.array_length, :intensity, 0.0, 0.0) - end - - spectra_metadata[k] = SpectrumMetadata(x, y, "", :sample, spectrum_mode, mz_asset, int_asset) - end - - return spectra_metadata -end - # ============================================================================= # diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..481ea7a --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,58 @@ +using Test + +# We need to simulate loading the MSI_src module directly +include("../src/MSI_src.jl") +using .MSI_src + +@testset "JuliaMSI Core Systems" begin + @testset "Basic Instantiation" begin + # 1. Test basic metadata structure initialization + meta = SpectrumMetadata( + Int32(1), Int32(2), + "test_id", + :sample, + CENTROID, + SpectrumAsset(Float64, false, Int64(0), 100, :mz, 0.0, 0.0), + SpectrumAsset(Float32, true, Int64(100), 50, :intensity, 0.0, 0.0) + ) + @test meta.x == 1 + @test meta.y == 2 + @test meta.mode == CENTROID + @test meta.mz_asset.format == Float64 + @test meta.int_asset.format == Float32 + @test meta.int_asset.is_compressed == true + + # 2. Test cache pool and MSIData structural stability + source = MzMLSource([], Float64, Float32, nothing) # Empty source for structural testing + + # Test constructor doesn't throw + # Constructor signature: (source, metadata, instrument_meta, dims, coordinate_map, cache_size) + msi_data = MSIData( + source, + [meta], + nothing, # instrument_meta + (100, 100), # dims + nothing, # coord map + 10 # cache size + ) + + @test msi_data.image_dims == (100, 100) + @test length(msi_data.spectra_metadata) == 1 + @test msi_data.cache_size == 10 + end + + @testset "Buffer & Cache Subsystems" begin + pool = SimpleBufferPool() + # Test basic retrieval + buf = get_buffer!(pool, 1024) + @test length(buf) == 1024 + + # Test release + release_buffer!(pool, buf) + @test length(pool.buffers[1024]) == 1 + + # Test reuse + buf2 = get_buffer!(pool, 1024) + @test buf === buf2 # Should return the EXACT same buffer object + end +end -- 2.43.0 From 89de5a56e09b6d370104cf7cf1a0f40d14bf5270 Mon Sep 17 00:00:00 2001 From: Pixelguy14 Date: Thu, 16 Apr 2026 16:39:31 -0600 Subject: [PATCH 6/6] Fixed small issues within the file load and the precompiler script, added a final test benchmark suit to compare old Julia MSI library with current --- precompile_script.jl | 15 ++- src/imzML.jl | 7 +- src/mzML.jl | 4 - test/benchmark.jl | 12 ++- test/benchmark_v3.jl | 215 +++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 231 insertions(+), 22 deletions(-) create mode 100644 test/benchmark_v3.jl diff --git a/precompile_script.jl b/precompile_script.jl index d1f7d6a..26d71c9 100644 --- a/precompile_script.jl +++ b/precompile_script.jl @@ -19,9 +19,9 @@ function create_warmup_datasets(dir, T_mz::Type, T_int::Type) mz_bytes = n_points * sizeof(T_mz) int_bytes = n_points * sizeof(T_int) - # Write some dummy binary data + # Write some dummy binary data (m/z must be sorted for validation) open(ibd_path, "w") do f - write(f, rand(T_mz, n_points)) + write(f, sort(rand(T_mz, n_points))) write(f, rand(T_int, n_points)) end @@ -83,7 +83,7 @@ function create_warmup_datasets(dir, T_mz::Type, T_int::Type) # 2. Create minimal mzML (Base64) mzml_path = joinpath(dir, "warmup_$(T_mz)_$(T_int).mzML") - b64_mz = base64encode(rand(T_mz, n_points)) + b64_mz = base64encode(sort(rand(T_mz, n_points))) b64_int = base64encode(rand(T_int, n_points)) mzml_content = """ @@ -140,8 +140,7 @@ try MSI_src.StreamingStep(:normalization, Dict(:method => :tic)), MSI_src.StreamingStep(:peak_picking, Dict(:method => :profile, :snr_threshold => 3.0)) ], - adaptive_binning=true, - bin_tolerance_ppm=50.0 + num_bins=2000 ) try MSI_src.process_dataset!(imzml_data, config) @@ -167,8 +166,8 @@ try p = PlotlyBase.Plot(PlotlyBase.scatter(x=1:2, y=[1,2])) end catch e - @warn "Global Warmup failed: $e" -catch - # Silent fail for interrupt + if !(e isa InterruptException) + @warn "Global Warmup failed: $e" + end end println("Precompilation workload completed.") diff --git a/src/imzML.jl b/src/imzML.jl index 101f6fd..a144893 100644 --- a/src/imzML.jl +++ b/src/imzML.jl @@ -499,11 +499,6 @@ function load_imzml_lazy(file_path::String; cache_size::Int=100, use_mmap::Bool= if use_mmap try file_size = filesize(ibd_path) - free_ram = Sys.free_memory() - - if file_size > free_ram * 0.8 - @warn "Dataset size ($(round(file_size/1e9, digits=2)) GB) exceeds 80% of free RAM. Mmap will still work via 'Streaming', but expect slight I/O overhead." - end @debug "Memory mapping .ibd file..." # We use the first handle for mmapping @@ -773,7 +768,7 @@ function get_mz_slice(data::MSIData, mass::Real, tolerance::Real; mask_path::Uni end end - println("Populated $(results_count[]) pixels with intensity data") + #println("Populated $(results_count[]) pixels with intensity data") replace!(slice_matrix, NaN => 0.0) return slice_matrix end diff --git a/src/mzML.jl b/src/mzML.jl index f9904b8..d3afd94 100644 --- a/src/mzML.jl +++ b/src/mzML.jl @@ -417,10 +417,6 @@ function load_mzml_lazy(file_path::String; cache_size::Int=100) mmap_data = nothing try file_size = filesize(file_path) - free_ram = Sys.free_memory() - if file_size > free_ram * 0.8 - @warn "Dataset size ($(round(file_size/1e9, digits=2)) GB) exceeds 80% of free RAM. Mmap will still work via 'Streaming' mode." - end println("DEBUG: Memory mapping .mzML file...") seekstart(primary_handle) # Anchor Mmap to the beginning of the file to prevent overflow diff --git a/test/benchmark.jl b/test/benchmark.jl index 8622a8d..a78e74f 100644 --- a/test/benchmark.jl +++ b/test/benchmark.jl @@ -49,12 +49,16 @@ const BENCHMARK_CASES = [ # BenchmarkCase("/path/to/your/small_file.imzML", 309.06, 0.1), # BenchmarkCase("/path/to/your/medium_file.imzML", 896.0, 1.0), # BenchmarkCase("/path/to/your/large_file.imzML", 100.0, 0.1), - BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_LA-ESI/180817_NEG_Thaliana_Leaf_bottom_1_0841.imzML",116.07,0.1, "Thaliana Leaf"), - BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_LTP/ltpmsi-chilli.imzML",420,0.1, "Chilli Pepper"), # Chilli - BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_DESI/ColAd_Individual/40TopL,10TopR,30BottomL,20BottomR/40TopL,10TopR,30BottomL,20BottomR-centroid.imzML",885.5,0.1, "Colon Cancer Human"), # Human Cancer + # BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_LA-ESI/180817_NEG_Thaliana_Leaf_bottom_1_0841.imzML",116.07,0.1, "Thaliana Leaf"), + BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/Chilli/ltpmsi-chilli.imzML",420,0.1, "Chilli Pepper"), # Chilli + BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_DESI/40TopL,10TopR,30BottomL,20BottomR/40TopL,10TopR,30BottomL,20BottomR-centroid.imzML",885.5,0.1, "Colon Cancer Human"), # Human Cancer BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML", 716.053,0.1, "Mouse Urinary Bladder"), # Mouse bladder + BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/Leafs/CE1_Leaf_R3.imzML",306.1,0.1, "Leaf"), BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/Liv2_imzML_TIMSConvert-selected/Liv2.imzML",796.18,0.1, "Liver Cut"), #Lib2 - BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML",804.3,0.1, "Mouse Stomach"), # Mouse Stomach + BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML",804.3,0.1, "Mouse Stomach 2GB"), # Mouse Stomach + BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/salida_Tims/Stomach_DHB.imzML",804.3,0.1, "Mouse Stomach 4GB"), # Mouse Stomach + + ] const NUM_REPETITIONS = 50 # Number of times to generate the image for averaging diff --git a/test/benchmark_v3.jl b/test/benchmark_v3.jl new file mode 100644 index 0000000..f668587 --- /dev/null +++ b/test/benchmark_v3.jl @@ -0,0 +1,215 @@ +# test/benchmark_v3.jl + +# =================================================================== +# High-Precision Performance Benchmark Suite (v3) +# =================================================================== +# This script integrates the robust statistical sampling of `BenchmarkTools` +# with the visual comparative mechanics against the legacy library. +# It captures the paradigm shift from "Time per slice" to "Pipeline Throughput". +# =================================================================== + +using Pkg +Pkg.activate(joinpath(@__DIR__, "..")) + +using BenchmarkTools +using DataFrames +using CSV +using CairoMakie +using Statistics +using MSI_src +using julia_mzML_imzML + +struct BenchmarkCase + filepath::String + mz_value::Float64 + mz_tolerance::Float64 + name::String +end + +const RESULTS_DIR = joinpath(@__DIR__, "results") + +# List to benchmark +const BENCHMARK_CASES = [ + BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/Chilli/ltpmsi-chilli.imzML", 420.0, 0.1, "Chilli Pepper"), + BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_DESI/40TopL,10TopR,30BottomL,20BottomR/40TopL,10TopR,30BottomL,20BottomR-centroid.imzML", 885.5, 0.1, "Colon Cancer Human"), + BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML", 716.053, 0.1, "Mouse Urinary Bladder"), + BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/Leafs/CE1_Leaf_R3.imzML", 306.1, 0.1, "Leaf"), + BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/Liv2_imzML_TIMSConvert-selected/Liv2.imzML", 796.18, 0.1, "Liver Cut"), + BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/salida_Tims/Stomach_DHB.imzML", 804.3, 0.1, "Mouse Stomach 4GB") +] + +function get_total_file_size_mb(filepath::String) + imzml_size = isfile(filepath) ? filesize(filepath) : 0 + ibd_path = replace(filepath, r"\.(imzML|imzml)$" => ".ibd") + ibd_size = isfile(ibd_path) ? filesize(ibd_path) : 0 + return round((imzml_size + ibd_size) / 1024^2, digits=2) +end + +function flush_memory() + GC.gc(true) + if Sys.islinux() + # Force the OS to reclaim memory from the glibc allocator + ccall(:malloc_trim, Int32, (Int32,), 0) + end +end + +function run_v3_benchmarks() + mkpath(RESULTS_DIR) + results = DataFrame() + + println("="^60) + println("STARTING ENTERPRISE BENCHMARK SUITE (v3)") + println("="^60) + + for case in BENCHMARK_CASES + if !isfile(case.filepath) + @warn "File not found: $(case.filepath). Skipping." + continue + end + + println("\n--- Target: $(case.name) ---") + file_size = get_total_file_size_mb(case.filepath) + + # ------------------------------------------------------------ + # 1. JuliaMSI (New Architecture) + # ------------------------------------------------------------ + println("[JuliaMSI - New Engine]") + flush_memory() + + # Load Phase (Metadata Only / Memory Mapping) + load_stats_new = @timed OpenMSIData(case.filepath) + msi_data = load_stats_new.value + load_time_new_s = load_stats_new.time + mem_load_new_mb = load_stats_new.bytes / 1024^2 + + # Slicing Phase (High Precision) + b_slice_new = @benchmark get_mz_slice($msi_data, $(case.mz_value), $(case.mz_tolerance)) samples=10 seconds=5 + mean_time_new_s = mean(b_slice_new.times) / 1e9 + + # Calculate Amortized Throughput (10 slices) + # Includes the 'Loading Wall' penalty + total_time_10_new = load_time_new_s + (10 * mean_time_new_s) + amortized_sps_new = 10.0 / total_time_10_new + + close(msi_data) + flush_memory() + + # ------------------------------------------------------------ + # 2. julia_mzML_imzML (Legacy Architecture) + # ------------------------------------------------------------ + println("[julia_mzML_imzML - Legacy]") + + load_stats_old = try + @timed LoadImzml(case.filepath) + catch e + @warn "Legacy load failed: $e" + (time=Inf, bytes=Inf, value=nothing) + end + + old_data = load_stats_old.value + mem_load_old_mb = load_stats_old.bytes / 1024^2 + + load_time_old_s = Inf + mean_time_old_s = Inf + amortized_sps_old = 0.0 + + if old_data !== nothing + b_slice_old = try + # Legacy can be extremely slow, constrain it heavily + @benchmark GetSlice($old_data, $(case.mz_value), $(case.mz_tolerance)) samples=3 seconds=10 + catch e + @warn "Legacy slice failed: $e" + nothing + end + + if b_slice_old !== nothing + mean_time_old_s = mean(b_slice_old.times) / 1e9 + + # Accurately reflect the massive load time in the throughput + load_time_old_s = load_stats_old.time + total_time_10_old = load_time_old_s + (10 * mean_time_old_s) + amortized_sps_old = 10.0 / total_time_10_old + end + end + + flush_memory() + + # ------------------------------------------------------------ + # Record + # ------------------------------------------------------------ + push!(results, ( + Dataset = case.name, + FileSize_MB = file_size, + + # Legacy Metrics + Legacy_LoadMem_MB = mem_load_old_mb, + Legacy_LoadTime_s = load_stats_old.time, + Legacy_Amortized_SPS_10 = amortized_sps_old, + + # New Metrics + New_LoadMem_MB = mem_load_new_mb, + New_LoadTime_s = load_time_new_s, + New_Amortized_SPS_10 = amortized_sps_new, + + # Competitive Deltas + RAM_Reduction_Pct = isfinite(mem_load_old_mb) ? ((mem_load_old_mb - mem_load_new_mb) / mem_load_old_mb)*100 : NaN, + UX_Speedup_Factor = isfinite(amortized_sps_old) ? (amortized_sps_new / amortized_sps_old) : NaN + )) + + println(" > RAM Reduction: $(round(results[end, :RAM_Reduction_Pct], digits=2))%") + println(" > Amortized Throughput (10 Slices): $(round(amortized_sps_old, digits=2)) -> $(round(amortized_sps_new, digits=2)) slices/sec") + end + + csv_path = joinpath(RESULTS_DIR, "v3_enterprise_benchmarks.csv") + CSV.write(csv_path, results) + println("\nData saved to $csv_path") + + plot_v3_results(results) + return results +end + +function plot_v3_results(df::DataFrame) + if isempty(df) return end + + fig = Figure(size=(1800, 1200), fontsize=24) + x_pos = 1:nrow(df) + labels = df.Dataset + + # ---------- Plot 1: The RAM Revolution ---------- + ax1 = Axis(fig[1, 1], + title="Initial RAM Cost (Loading & Mmap)", + ylabel="Memory (MB)", + xticks=(x_pos, labels), xticklabelrotation=π/8) + + barplot!(ax1, x_pos .- 0.2, df.New_LoadMem_MB, color="#10b981", width=0.4, label="JuliaMSI (Mmap Lazy)") + barplot!(ax1, x_pos .+ 0.2, df.Legacy_LoadMem_MB, color="#ef4444", width=0.4, label="Legacy (Dense Matrix)") + axislegend(ax1, position=:lt) + + # ---------- Plot 2: Throughput Leap ---------- + ax2 = Axis(fig[1, 2], + title="Amortized Throughput (10 Slices - Higher is Better)", + ylabel="Slices / Second (Inc. Load Time)", + xticks=(x_pos, labels), xticklabelrotation=π/8) + + barplot!(ax2, x_pos .- 0.2, df.New_Amortized_SPS_10, color="#10b981", width=0.4, label="JuliaMSI") + barplot!(ax2, x_pos .+ 0.2, df.Legacy_Amortized_SPS_10, color="#ef4444", width=0.4, label="Legacy") + axislegend(ax2, position=:lt) + + # ---------- Plot 3: The Scaling Wall (File Size vs Throughput) ---------- + ax3 = Axis(fig[2, 1:2], + title="Performance Scaling: UX Throughput vs Dataset Size", + xlabel="Dataset File Size (MB)", + ylabel="Amortized Throughput (Slices/sec)") + + scatterlines!(ax3, df.FileSize_MB, df.New_Amortized_SPS_10, color="#10b981", markersize=15, linewidth=4, label="JuliaMSI Engine") + scatterlines!(ax3, df.FileSize_MB, df.Legacy_Amortized_SPS_10, color="#ef4444", markersize=15, linewidth=4, label="Legacy Engine") + axislegend(ax3, position=:rt) + + save(joinpath(RESULTS_DIR, "v3_enterprise_dashboard.png"), fig) + println("\nDashboards saved successfully.") +end + +# Ensure we process if run as the main script +if abspath(PROGRAM_FILE) == @__FILE__ + run_v3_benchmarks() +end -- 2.43.0