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()