From edfbb77f58eab7eb571b2e0b544ba2ffb18fe2f1 Mon Sep 17 00:00:00 2001 From: Pixelguy14 Date: Sun, 9 Nov 2025 13:41:29 -0600 Subject: [PATCH] Added bloom filters custom API, thread safety for multithreading several important functions, streaming processing to reduce GC and Memory allocations calls, buffer pool changed with lock safety, and debug memory monitoring to track leaks, alongside better loading times in general --- Manifest.toml | 12 +- Project.toml | 42 +++ app.jl | 91 +++--- julia_imzML_visual.jl | 4 +- src/BloomFilters.jl | 205 +++++++++++++ src/Common.jl | 185 +++++++++++- src/MSIData.jl | 664 +++++++++++++++++++++++++++++++----------- src/MSI_src.jl | 10 +- src/imzML.jl | 493 +++++++++++++++---------------- src/mzML.jl | 40 ++- start_MSI_GUI.jl | 1 + 11 files changed, 1256 insertions(+), 491 deletions(-) create mode 100644 src/BloomFilters.jl diff --git a/Manifest.toml b/Manifest.toml index 4ec7d85..62323b8 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.11.7" manifest_format = "2.0" -project_hash = "f5ab935406d0b7b7e657445ecfa45944ea71f247" +project_hash = "38f13d997585a9af185f1863ff626d584a7e791e" [[deps.ATK_jll]] deps = ["Artifacts", "Glib_jll", "JLLWrappers", "Libdl"] @@ -651,6 +651,12 @@ git-tree-sha1 = "05882d6995ae5c12bb5f36dd2ed3f61c98cbb172" uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" version = "0.8.5" +[[deps.FlameGraphs]] +deps = ["AbstractTrees", "Colors", "FileIO", "FixedPointNumbers", "IndirectArrays", "LeftChildRightSiblingTrees", "Profile"] +git-tree-sha1 = "0166baf81babb91cf78bfcc771d8e87c43d568df" +uuid = "08572546-2f56-4bcf-ba4e-bab62c3a3f89" +version = "1.1.0" + [[deps.Fontconfig_jll]] deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Zlib_jll"] git-tree-sha1 = "f85dac9a96a01087df6e3a749840015a0ca3817d" @@ -1779,6 +1785,10 @@ deps = ["Unicode"] uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" version = "1.11.0" +[[deps.Profile]] +uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79" +version = "1.11.0" + [[deps.ProgressMeter]] deps = ["Distributed", "Printf"] git-tree-sha1 = "fbb92c6c56b34e1a2c4c36058f68f332bec840e7" diff --git a/Project.toml b/Project.toml index 5d7b777..242e031 100644 --- a/Project.toml +++ b/Project.toml @@ -12,6 +12,7 @@ Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" +FlameGraphs = "08572546-2f56-4bcf-ba4e-bab62c3a3f89" GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" Genie = "c43c736e-a2d1-11e8-161f-af95117fbd1e" GenieFramework = "a59fdf5c-6bf0-4f5d-949c-a137c9e2f353" @@ -19,6 +20,7 @@ HistogramThresholding = "2c695a8d-9458-5d45-9878-1b8a99cf7853" ImageBinarization = "cbc4b850-ae4b-5111-9e64-df94c024a13d" ImageComponentAnalysis = "d9b9e9a0-1569-11e9-2cb5-bbca914b0e89" ImageContrastAdjustment = "f332f351-ec65-5f6a-b3d1-319c6670881a" +ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534" ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5" ImageMorphology = "787d08f9-d448-5407-9aad-5290dd7ab264" ImageSegmentation = "80713f31-8817-5129-9cf8-209ff8fb23e1" @@ -39,3 +41,43 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StipplePlotly = "ec984513-233d-481d-95b0-a3b58b97af2b" UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" + +[compat] +Accessors = "0.1" +Base64 = "1.11" +CSV = "0.10" +CairoMakie = "0.13" +ColorSchemes = "3.31" +Colors = "0.12" +DataFrames = "1.7" +Dates = "1.11" +FileIO = "1.17" +GLMakie = "0.11" +Genie = "5.31" +GenieFramework = "2.8" +HistogramThresholding = "0.3" +ImageBinarization = "0.3" +ImageComponentAnalysis = "0.2" +ImageContrastAdjustment = "0.3" +ImageCore = "0.10.5" +ImageFiltering = "0.7" +ImageMorphology = "0.4" +ImageSegmentation = "1.9" +Images = "0.26" +JSON = "0.21" +Libz = "1.0" +LinearAlgebra = "1.11" +Loess = "0.6" +Mmap = "1.11" +NativeFileDialog = "0.2" +NaturalSort = "1.0" +PlotlyBase = "0.8" +Printf = "1.11" +ProgressMeter = "1.11" +SavitzkyGolay = "0.9" +Serialization = "1.11" +Statistics = "1.11" +StatsBase = "0.34" +StipplePlotly = "0.13" +UUIDs = "1.11" +julia = "1.11" diff --git a/app.jl b/app.jl index 1e96ee7..5e4a5ef 100644 --- a/app.jl +++ b/app.jl @@ -19,9 +19,10 @@ using Base.Filesystem: mv # To rename files in the system using Printf # Required for @sprintf macro in colorbar generation using JSON using Dates +using Base.Threads # Bring MSIData into App module's scope -using .MSI_src: MSIData, OpenMSIData, process_spectrum, IterateSpectra, ImzMLSource, _iterate_spectra_fast, MzMLSource, find_mass, ViridisPalette, get_mz_slice, get_multiple_mz_slices, quantize_intensity, save_bitmap, median_filter, save_bitmap, downsample_spectrum, TrIQ, precompute_analytics, ImportMzmlFile, generate_colorbar_image, load_and_prepare_mask +using .MSI_src: MSIData, OpenMSIData, process_spectrum, IterateSpectra, ImzMLSource, _iterate_spectra_fast, MzMLSource, find_mass, ViridisPalette, get_mz_slice, get_multiple_mz_slices, quantize_intensity, save_bitmap, median_filter, save_bitmap, downsample_spectrum, TrIQ, precompute_analytics, ImportMzmlFile, generate_colorbar_image, load_and_prepare_mask, set_global_mz_range! if !@isdefined(increment_image) include("./julia_imzML_visual.jl") @@ -837,7 +838,7 @@ end end end - @onbutton createMeanPlot begin + @onbutton createMeanPlot @time begin if isempty(selected_folder_main) msg = "No dataset selected. Please process a file and select a folder first." warning_msg = true @@ -862,12 +863,18 @@ end end if msi_data === nothing || full_route != target_path + if msi_data !== nothing + close(msi_data) + end msg = "Reloading $(basename(target_path)) for analysis..." full_route = target_path msi_data = OpenMSIData(target_path) if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing - msi_data.global_min_mz = entry["metadata"]["global_min_mz"] - msi_data.global_max_mz = entry["metadata"]["global_max_mz"] + 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)) else precompute_analytics(msi_data) end @@ -904,7 +911,7 @@ end end end - @onbutton createSumPlot begin + @onbutton createSumPlot @time begin if isempty(selected_folder_main) msg = "No dataset selected. Please process a file and select a folder first." warning_msg = true @@ -914,32 +921,37 @@ end progressSpectraPlot = true btnPlotDisable = true btnStartDisable = true - msg = "Loading total spectrum plot for $(selected_folder_main)..." - - try - sTime = time() - registry = load_registry(registry_path) - entry = registry[selected_folder_main] - target_path = entry["source_path"] - - if target_path == "unknown (manually added)" - msg = "Dataset selected contained no route." - warning_msg = true - return - end - - if msi_data === nothing || full_route != target_path - msg = "Reloading $(basename(target_path)) for analysis..." - full_route = target_path - msi_data = OpenMSIData(target_path) - if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing - msi_data.global_min_mz = entry["metadata"]["global_min_mz"] - msi_data.global_max_mz = entry["metadata"]["global_max_mz"] - else - precompute_analytics(msi_data) - end - end - + msg = "Loading total spectrum plot for $(selected_folder_main)..." + + try + sTime = time() + registry = load_registry(registry_path) + entry = registry[selected_folder_main] + target_path = entry["source_path"] + + if target_path == "unknown (manually added)" + msg = "Dataset selected contained no route." + warning_msg = true + return + end + + if msi_data === nothing || full_route != target_path + if msi_data !== nothing + close(msi_data) + end + msg = "Reloading $(basename(target_path)) for analysis..." + full_route = target_path + 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)) + else + precompute_analytics(msi_data) + end + end local mask_path_for_plot::Union{String, Nothing} = nothing if maskEnabled && get(entry, "has_mask", false) mask_path_for_plot = get(entry, "mask_path", "") @@ -971,7 +983,7 @@ end end end - @onbutton createXYPlot begin + @onbutton createXYPlot @time begin if isempty(selected_folder_main) msg = "No dataset selected. Please process a file and select a folder first." warning_msg = true @@ -1005,12 +1017,18 @@ end end if msi_data === nothing || full_route != target_path + if msi_data !== nothing + close(msi_data) + end msg = "Reloading $(basename(target_path)) for analysis..." full_route = target_path msi_data = OpenMSIData(target_path) if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing - msi_data.global_min_mz = entry["metadata"]["global_min_mz"] - msi_data.global_max_mz = entry["metadata"]["global_max_mz"] + 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)) else precompute_analytics(msi_data) end @@ -1351,6 +1369,9 @@ end # This handler will now correctly load the first image from the newly selected folder. @onchange selected_folder_main begin + if msi_data !== nothing + close(msi_data) + end msi_data = nothing log_memory_usage("Folder Changed (msi_data cleared)", msi_data) GC.gc() @@ -1900,7 +1921,7 @@ end @mounted watchplots() - @onchange isready begin + @onchange isready @time begin if isready && !registry_init_done warmup_init() try diff --git a/julia_imzML_visual.jl b/julia_imzML_visual.jl index e18f6a5..b9b40c4 100644 --- a/julia_imzML_visual.jl +++ b/julia_imzML_visual.jl @@ -847,8 +847,8 @@ function extract_metadata(msi_data::MSIData, source_path::String) Dict("parameter" => "File Name", "value" => basename(source_path)), Dict("parameter" => "Number of Spectra", "value" => length(msi_data.spectra_metadata)), Dict("parameter" => "Image Dimensions", "value" => "$(msi_data.image_dims[1]) x $(msi_data.image_dims[2])"), - Dict("parameter" => "Global Min m/z", "value" => @sprintf("%.4f", msi_data.global_min_mz)), - Dict("parameter" => "Global Max m/z", "value" => @sprintf("%.4f", msi_data.global_max_mz)), + Dict("parameter" => "Global Min m/z", "value" => @sprintf("%.4f", Threads.atomic_add!(msi_data.global_min_mz, 0.0))), + Dict("parameter" => "Global Max m/z", "value" => @sprintf("%.4f", Threads.atomic_add!(msi_data.global_max_mz, 0.0))), Dict("parameter" => "Mean TIC", "value" => @sprintf("%.2e", mean(df.TIC))), Dict("parameter" => "Mean BPI", "value" => @sprintf("%.2e", mean(df.BPI))), Dict("parameter" => "Mean # Points", "value" => @sprintf("%.1f", mean(df.NumPoints))), diff --git a/src/BloomFilters.jl b/src/BloomFilters.jl new file mode 100644 index 0000000..2d95ed2 --- /dev/null +++ b/src/BloomFilters.jl @@ -0,0 +1,205 @@ +# src/BloomFilter.jl + +""" + BloomFilter{T} + +A simple, efficient Bloom filter implementation for probabilistic set membership testing. + +# Fields +- `bits::BitVector`: The underlying bit array +- `size::Int`: Number of bits in the filter +- `hash_count::Int`: Number of hash functions to use +- `seed::UInt64`: Random seed for hash functions +- `count::Int`: Number of elements added (for monitoring) +""" +mutable struct BloomFilter{T} + bits::BitVector + size::Int + hash_count::Int + seed::UInt64 + count::Int +end + +""" + StreamingBloomFilter + +A memory-efficient Bloom filter that doesn't store all values at once. +""" +mutable struct StreamingBloomFilter + bits::BitVector + size::Int + hash_count::Int + seed::UInt64 + count::Int + # Streaming state + current_chunk::Vector{Int} + chunk_size::Int +end + +""" + BloomFilter(expected_elements::Int, false_positive_rate::Float64=0.01; seed::UInt64=0x12345678) + +Creates a Bloom filter optimized for the expected number of elements and desired false positive rate. +""" +function BloomFilter{T}(expected_elements::Int, false_positive_rate::Float64=0.01; seed::Union{UInt32,UInt64}=0x12345678) where T + # Convert seed to UInt64 for consistency + seed_uint64 = UInt64(seed) + + # Calculate optimal parameters using standard formulas + size = optimal_bit_size(expected_elements, false_positive_rate) + hash_count = optimal_hash_count(expected_elements, size) + + bits = falses(size) + return BloomFilter{T}(bits, size, hash_count, seed_uint64, 0) +end + +""" + optimal_bit_size(n::Int, p::Float64) -> Int + +Calculates the optimal number of bits for a Bloom filter given: +- n: expected number of elements +- p: desired false positive rate +""" +function optimal_bit_size(n::Int, p::Float64)::Int + if p <= 0.0 || p >= 1.0 + throw(ArgumentError("False positive rate must be between 0 and 1")) + end + # m = - (n * ln(p)) / (ln(2)^2) + m = ceil(Int, - (n * log(p)) / (log(2)^2)) + return max(m, 1) +end + +""" + optimal_hash_count(n::Int, m::Int) -> Int + +Calculates the optimal number of hash functions for a Bloom filter. +""" +function optimal_hash_count(n::Int, m::Int)::Int + if n <= 0 || m <= 0 + return 1 + end + # k = (m / n) * ln(2) + k = max(1, round(Int, (m / n) * log(2))) + return min(k, 8) # Practical limit to avoid too many hashes +end + +""" + hash_functions(item::T, count::Int, size::Int, seed::UInt64) -> Vector{Int} + +Generates multiple hash values for an item using double hashing technique. +""" +function hash_functions(item::T, count::Int, size::Int, seed::UInt64) where T + # Use Julia's built-in hash with different seeds + hashes = Vector{Int}(undef, count) + + # First hash with the main seed + h1 = hash(item, seed) + h2 = hash(item, seed + 1) + + for i in 1:count + # Double hashing: h_i = h1 + i * h2 + combined_hash = UInt64(h1) + UInt64(i) * UInt64(h2) + hashes[i] = (combined_hash % UInt64(size)) + 1 # 1-based indexing + end + + return hashes +end + +""" + Base.push!(bf::BloomFilter{T}, item::T) + +Adds an element to the Bloom filter. +""" +function Base.push!(bf::BloomFilter{T}, item::T) where T + hashes = hash_functions(item, bf.hash_count, bf.size, bf.seed) + + for h in hashes + bf.bits[h] = true + end + bf.count += 1 + + return bf +end + +""" + Base.in(item::T, bf::BloomFilter{T}) -> Bool + +Checks whether an element is possibly in the Bloom filter. +Returns `true` if the element might be in the set, `false` if it's definitely not. +""" +function Base.in(item::T, bf::BloomFilter{T})::Bool where T + hashes = hash_functions(item, bf.hash_count, bf.size, bf.seed) + + for h in hashes + if !bf.bits[h] + return false # Definitely not in the set + end + end + return true # Possibly in the set +end + +""" + contains(bf::BloomFilter{T}, item::T) -> Bool + +Alias for `in()` for compatibility with your existing code. +""" +contains(bf::BloomFilter{T}, item::T) where T = item in bf + +""" + add!(bf::BloomFilter{T}, item::T) + +Alias for `push!()` for compatibility with your existing code. +""" +add!(bf::BloomFilter{T}, item::T) where T = push!(bf, item) + +""" + false_positive_rate(bf::BloomFilter) -> Float64 + +Estimates the current false positive rate of the Bloom filter. +""" +function false_positive_rate(bf::BloomFilter)::Float64 + if bf.count == 0 + return 0.0 + end + # Theoretical false positive rate: (1 - e^(-k * n / m)) ^ k + k = bf.hash_count + n = bf.count + m = bf.size + return (1 - exp(-k * n / m)) ^ k +end + +""" + fill_ratio(bf::BloomFilter) -> Float64 + +Returns the fraction of bits that are set to 1. +""" +function fill_ratio(bf::BloomFilter)::Float64 + return count(bf.bits) / length(bf.bits) +end + +""" + is_empty(bf::BloomFilter) -> Bool + +Checks if the Bloom filter is empty (no elements added). +""" +function is_empty(bf::BloomFilter)::Bool + return bf.count == 0 +end + +""" + reset!(bf::BloomFilter) + +Clears the Bloom filter, removing all elements. +""" +function reset!(bf::BloomFilter) + fill!(bf.bits, false) + bf.count = 0 + return bf +end + +# Specialized constructor for empty Bloom filters +function BloomFilter{T}(;size::Int=100, hash_count::Int=3, seed::Union{UInt32,UInt64}=0x12345678) where T + seed_uint64 = UInt64(seed) + bits = falses(size) + return BloomFilter{T}(bits, size, hash_count, seed_uint64, 0) +end diff --git a/src/Common.jl b/src/Common.jl index 2ecea27..8de5929 100644 --- a/src/Common.jl +++ b/src/Common.jl @@ -1,4 +1,114 @@ -# src/Common.jl +# src/Common.jl - Updated with BloomFilter +using Base.Threads + +# --- Buffer Pooling --- +""" + BufferPool + +A thread-safe pool of `Vector{UInt8}` buffers, categorized by their size. +This helps reduce memory allocations and garbage collection overhead by reusing buffers. + +# Fields +- `lock`: A `ReentrantLock` to ensure thread-safe access to the pool. +- `buffers`: A dictionary mapping buffer size (in bytes) to a list of available buffers of that size. +""" +mutable struct BufferPool + lock::ReentrantLock + buffers::Dict{Int, Vector{Vector{UInt8}}} + + function BufferPool() + new(ReentrantLock(), Dict{Int, Vector{Vector{UInt8}}}()) + end +end + +""" + get_buffer!(pool::BufferPool, size::Int) -> Vector{UInt8} + +Retrieves a `Vector{UInt8}` buffer of at least `size` bytes from the pool. +If no suitable buffer is available, a new one is allocated. + +# Arguments +- `pool`: The `BufferPool` to retrieve the buffer from. +- `size`: The minimum desired size of the buffer in bytes. + +# Returns +- A `Vector{UInt8}` buffer. +""" +function get_buffer!(pool::BufferPool, size::Int)::Vector{UInt8} + lock(pool.lock) do + # Find an existing buffer that is large enough + for (buffer_size, buffer_list) in pool.buffers + if buffer_size >= size && !isempty(buffer_list) + buffer = pop!(buffer_list) + # Resize if necessary (and if it's significantly larger than needed) + if length(buffer) > size * 2 # Heuristic: if buffer is more than twice the requested size + return Vector{UInt8}(undef, size) # Allocate new smaller buffer + else + return buffer + end + end + end + # No suitable buffer found, allocate a new one + return Vector{UInt8}(undef, size) + end +end + +""" + release_buffer!(pool::BufferPool, buffer::Vector{UInt8}) + +Returns a `Vector{UInt8}` buffer to the pool for future reuse. + +# Arguments +- `pool`: The `BufferPool` to return the buffer to. +- `buffer`: The `Vector{UInt8}` buffer to release. +""" +function release_buffer!(pool::BufferPool, buffer::Vector{UInt8}) + lock(pool.lock) do + buffer_size = length(buffer) + if !haskey(pool.buffers, buffer_size) + pool.buffers[buffer_size] = Vector{Vector{UInt8}}() + end + push!(pool.buffers[buffer_size], buffer) + end + return nothing +end + + +""" + SimpleBufferPool + +A dramatically simplified buffer pool that avoids complex locking. +""" +mutable struct SimpleBufferPool + buffers::Dict{Int, Vector{Vector{UInt8}}} + max_pool_size::Int +end + +SimpleBufferPool() = SimpleBufferPool(Dict{Int, Vector{Vector{UInt8}}}(), 50) + +function get_buffer!(pool::SimpleBufferPool, size::Int)::Vector{UInt8} + # Check for existing buffers of exact size first + if haskey(pool.buffers, size) && !isempty(pool.buffers[size]) + return pop!(pool.buffers[size]) + end + + # No suitable buffer found, allocate new one + return Vector{UInt8}(undef, size) +end + +function release_buffer!(pool::SimpleBufferPool, buffer::Vector{UInt8}) + size = length(buffer) + + if !haskey(pool.buffers, size) + pool.buffers[size] = Vector{Vector{UInt8}}() + end + + # Limit pool size to prevent memory bloat + if length(pool.buffers[size]) < pool.max_pool_size + push!(pool.buffers[size], buffer) + end + # If pool is full, let buffer get GC'd +end # --- Unified Error Types --- abstract type MSIError <: Exception end @@ -19,19 +129,13 @@ end # --- Shared Constants --- const DEFAULT_CACHE_SIZE = 100 const DEFAULT_NUM_BINS = 2000 - +const DEFAULT_BLOOM_FILTER_SIZE = 1000 # Default size for empty spectra +const DEFAULT_FALSE_POSITIVE_RATE = 0.01 """ validate_spectrum_data(mz, intensity, id) Checks a spectrum for common data integrity issues. - -- Checks for `NaN` or `Inf` values in intensity and m/z arrays. -- Verifies that the m/z array is sorted in ascending order. -- Ensures that m/z and intensity arrays have the same length. - -# Throws -- `InvalidSpectrumError` if any of the checks fail. """ function validate_spectrum_data(mz::AbstractVector, intensity::AbstractVector, id) if length(mz) != length(intensity) @@ -48,3 +152,66 @@ function validate_spectrum_data(mz::AbstractVector, intensity::AbstractVector, i return true end + +""" + create_bloom_filter_for_spectrum(mz::AbstractVector{<:Real}) -> BloomFilter{Int} + +Memory-optimized version that processes data in chunks to reduce temporary allocations +and avoid holding large intermediate arrays. +""" +function create_bloom_filter_for_spectrum(mz::AbstractVector{<:Real})::BloomFilter{Int} + if isempty(mz) + return empty_bloom_filter() + end + + # Size Bloom filter appropriately + n_points = length(mz) + bf = BloomFilter{Int}(n_points, DEFAULT_FALSE_POSITIVE_RATE) + + discretization_factor = 100.0 + chunk_size = 10000 # Process in chunks to reduce memory pressure + + # Process spectrum data in chunks to avoid large temporary arrays + for chunk_start in 1:chunk_size:n_points + chunk_end = min(chunk_start + chunk_size - 1, n_points) + + @inbounds for i in chunk_start:chunk_end + discretized_val = round(Int, mz[i] * discretization_factor) + push!(bf, discretized_val) + end + end + + return bf +end + +""" + empty_bloom_filter() -> BloomFilter{Float64} + +Creates an empty Bloom filter for spectra with no data. +""" +function empty_bloom_filter()::BloomFilter{Int} + return BloomFilter{Int}(size=DEFAULT_BLOOM_FILTER_SIZE, hash_count=3) +end + +""" + AtomicFlag + +A thread-safe boolean flag using atomic operations. +""" +struct AtomicFlag + value::Base.Threads.Atomic{Int} +end + +AtomicFlag() = AtomicFlag(Base.Threads.Atomic{Int}(0)) + +function set!(flag::AtomicFlag) + Base.Threads.atomic_add!(flag.value, 1) +end + +function is_set(flag::AtomicFlag)::Bool + return Base.Threads.atomic_add!(flag.value, 0) > 0 +end + +function reset!(flag::AtomicFlag) + Base.Threads.atomic_xchg!(flag.value, 0) +end diff --git a/src/MSIData.jl b/src/MSIData.jl index 65970aa..ed47a8b 100644 --- a/src/MSIData.jl +++ b/src/MSIData.jl @@ -8,6 +8,64 @@ efficiently. using Base64, Libz, Serialization, Printf, DataFrames, Base.Threads # For reading binary data +const FILE_HANDLE_LOCK = ReentrantLock() + +""" + ThreadSafeFileHandle + +Wraps a file handle with thread-safe operations. +""" +mutable struct ThreadSafeFileHandle + path::String + handle::IO + lock::ReentrantLock +end + +function ThreadSafeFileHandle(path::String, mode::String="r") + handle = lock(FILE_HANDLE_LOCK) do + open(path, mode) + end + return ThreadSafeFileHandle(path, handle, ReentrantLock()) +end + +function Base.seek(tsfh::ThreadSafeFileHandle, pos::Integer) + lock(tsfh.lock) do + seek(tsfh.handle, pos) + end +end + +function Base.read(tsfh::ThreadSafeFileHandle, n::Integer) + lock(tsfh.lock) do + return read(tsfh.handle, n) + end +end + +function Base.read!(tsfh::ThreadSafeFileHandle, a::AbstractArray) + lock(tsfh.lock) do + return read!(tsfh.handle, a) + end +end + +function Base.close(tsfh::ThreadSafeFileHandle) + close(tsfh.handle) +end + +function Base.isopen(tsfh::ThreadSafeFileHandle) + isopen(tsfh.handle) +end + +function Base.position(tsfh::ThreadSafeFileHandle) + lock(tsfh.lock) do + return position(tsfh.handle) + end +end + +function Base.filesize(tsfh::ThreadSafeFileHandle) + lock(tsfh.lock) do + return filesize(tsfh.handle) + end +end + # Abstract type for different data sources (e.g., mzML, imzML) # This allows dispatching to the correct binary reading logic. abstract type MSDataSource end @@ -20,7 +78,7 @@ 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::IO + ibd_handle::Union{IO, ThreadSafeFileHandle} mz_format::Type intensity_format::Type end @@ -32,7 +90,7 @@ 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::IO + file_handle::Union{IO, ThreadSafeFileHandle} mz_format::Type intensity_format::Type end @@ -122,24 +180,31 @@ efficient repeated access to spectra. mutable struct MSIData source::MSDataSource spectra_metadata::Vector{SpectrumMetadata} - image_dims::Tuple{Int, Int} # (width, height) for imaging data - coordinate_map::Union{Matrix{Int}, Nothing} # Maps (x,y) to linear index for imzML + image_dims::Tuple{Int, Int} + coordinate_map::Union{Matrix{Int}, Nothing} # LRU Cache for GetSpectrum cache::Dict{Int, Tuple{Vector{Float64}, Vector{Float64}}} - cache_order::Vector{Int} # Stores indices, with most recently used at the end - cache_size::Int # Max number of spectra in cache - cache_lock::ReentrantLock # To make cache access thread-safe + cache_order::Vector{Int} + cache_size::Int + cache_lock::ReentrantLock - # Pre-computed analytics/metadata - global_min_mz::Union{Float64, Nothing} - global_max_mz::Union{Float64, Nothing} + # Buffer Pool for binary data operations + buffer_pool::SimpleBufferPool + + # Pre-computed analytics/metadata - use Threads.Atomic for compatibility + global_min_mz::Threads.Atomic{Float64} + global_max_mz::Threads.Atomic{Float64} spectrum_stats_df::Union{DataFrame, Nothing} + bloom_filters::Union{Vector{<:BloomFilter}, Nothing} + analytics_ready::AtomicFlag function MSIData(source, metadata, dims, coordinate_map, cache_size) obj = new(source, metadata, dims, coordinate_map, Dict(), [], cache_size, ReentrantLock(), - nothing, nothing, nothing) # Initialize new fields to nothing + SimpleBufferPool(), + Threads.Atomic{Float64}(Inf), Threads.Atomic{Float64}(-Inf), + nothing, nothing, AtomicFlag()) # Ensure file handles are closed when the object is garbage collected finalizer(obj) do o @@ -153,6 +218,38 @@ mutable struct MSIData end end +# Thread-safe state management +function set_global_mz_range!(data::MSIData, min_mz::Float64, max_mz::Float64) + Threads.atomic_xchg!(data.global_min_mz, min_mz) + Threads.atomic_xchg!(data.global_max_mz, max_mz) +end + +function get_global_mz_range(data::MSIData) + # Atomic read using atomic_add! with zero + min_mz = Threads.atomic_add!(data.global_min_mz, 0.0) + max_mz = Threads.atomic_add!(data.global_max_mz, 0.0) + return (min_mz, max_mz) +end + +function set_analytics_data!(data::MSIData, stats_df::DataFrame, bloom_filters::Vector) + lock(data.cache_lock) do + data.spectrum_stats_df = stats_df + data.bloom_filters = bloom_filters + end +end + +function get_bloom_filters(data::MSIData) + lock(data.cache_lock) do + return data.bloom_filters + end +end + +function get_spectrum_stats(data::MSIData) + lock(data.cache_lock) do + return data.spectrum_stats_df + end +end + """ Base.close(data::MSIData) @@ -286,22 +383,40 @@ by this function and is assumed to be handled by the caller if necessary. # Returns - A `Vector` of the appropriate type containing the decoded data. """ -function read_binary_vector(io::IO, asset::SpectrumAsset) +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))")) end + seek(io, asset.offset) raw_b64 = read(io, asset.encoded_length) - decoded_bytes = base64decode(strip(String(raw_b64))) - bytes_io = IOBuffer(asset.is_compressed ? Libz.inflate(decoded_bytes) : decoded_bytes) - n_elements = Int(bytes_io.size / sizeof(asset.format)) - out_array = Array{asset.format}(undef, n_elements) - read!(bytes_io, out_array) - - # FIX: little-endian to host byte order for mzML data - out_array .= ltoh.(out_array) - + # Use String directly to avoid intermediate allocations + b64_string = String(raw_b64) + + 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) + out_array = Vector{asset.format}(undef, n_elements) + + # Use unsafe_wrap to avoid copy - this is the critical optimization + temp_array = unsafe_wrap(Vector{asset.format}, pointer(decoded_bytes), n_elements) + + # Convert byte order in-place + @inbounds for i in 1:n_elements + out_array[i] = ltoh(temp_array[i]) + end + return out_array end @@ -344,11 +459,9 @@ function read_spectrum_from_disk(source::ImzMLSource, meta::SpectrumMetadata) mz .= ltoh.(mz) intensity .= ltoh.(intensity) # FIX: Added missing conversion for intensity - mz_f64, intensity_f64 = Float64.(mz), Float64.(intensity) + validate_spectrum_data(mz, intensity, meta.id) - validate_spectrum_data(mz_f64, intensity_f64, meta.id) - - return mz_f64, intensity_f64 + return mz, intensity end """ @@ -365,16 +478,16 @@ This is a high-level function that orchestrates the reading process by calling # Returns - A tuple `(mz, intensity)` containing the two requested data arrays. """ -function read_spectrum_from_disk(source::MzMLSource, meta::SpectrumMetadata) +function read_spectrum_from_disk(data::MSIData, source::MzMLSource, meta::SpectrumMetadata) # For mzML, data is Base64 encoded within the XML - mz = read_binary_vector(source.file_handle, meta.mz_asset) - intensity = read_binary_vector(source.file_handle, meta.int_asset) + mz = read_binary_vector(data, source.file_handle, meta.mz_asset) + intensity = read_binary_vector(data, source.file_handle, meta.int_asset) mz_f64, intensity_f64 = Float64.(mz), Float64.(intensity) - validate_spectrum_data(mz_f64, intensity_f64, meta.id) + validate_spectrum_data(mz, intensity, meta.id) - return mz_f64, intensity_f64 + return mz, intensity end # --- Public API --- # @@ -414,7 +527,13 @@ function GetSpectrum(data::MSIData, index::Int) # Phase 2: Cache Miss - Read from disk (no lock) meta = data.spectra_metadata[index] - spectrum = read_spectrum_from_disk(data.source, meta) + + local spectrum + if data.source isa ImzMLSource + spectrum = read_spectrum_from_disk(data.source, meta) + else # MzMLSource + spectrum = read_spectrum_from_disk(data, data.source, meta) + end # Phase 3: Update cache (with lock) and get final value return lock(data.cache_lock) do @@ -527,122 +646,118 @@ as they can use this cached data. This function modifies the `MSIData` object in and is idempotent. """ function precompute_analytics(msi_data::MSIData) - # Idempotency check - if msi_data.spectrum_stats_df !== nothing && hasproperty(msi_data.spectrum_stats_df, :MinMZ) - println("Analytics have already been pre-computed.") + # Double-checked locking pattern + if is_set(msi_data.analytics_ready) return end - println("Pre-computing analytics (single pass)...") - start_time = time_ns() - - num_spectra = length(msi_data.spectra_metadata) - - # DEBUG: Check the first spectrum's metadata - if num_spectra > 0 - first_meta = msi_data.spectra_metadata[1] - println("DEBUG First spectrum metadata:") - println(" mz_asset: format=$(first_meta.mz_asset.format), compressed=$(first_meta.mz_asset.is_compressed)") - println(" mz_asset: offset=$(first_meta.mz_asset.offset), encoded_length=$(first_meta.mz_asset.encoded_length)") - println(" int_asset: format=$(first_meta.int_asset.format), compressed=$(first_meta.int_asset.is_compressed)") - println(" int_asset: offset=$(first_meta.int_asset.offset), encoded_length=$(first_meta.int_asset.encoded_length)") - println(" mode: $(first_meta.mode)") - end - - # Initialize variables for global stats - g_min_mz = Inf - g_max_mz = -Inf - - # Initialize vectors for per-spectrum stats - tics = Vector{Float64}(undef, num_spectra) - bpis = Vector{Float64}(undef, num_spectra) - bp_mzs = Vector{Float64}(undef, num_spectra) - num_points = Vector{Int}(undef, num_spectra) - min_mzs = Vector{Float64}(undef, num_spectra) - max_mzs = Vector{Float64}(undef, num_spectra) - modes = Vector{SpectrumMode}(undef, num_spectra) - is_compressed = Vector{Bool}(undef, num_spectra) - - # DEBUG: Add counter to see how many spectra have data - spectra_with_data = 0 - empty_spectra = 0 - - _iterate_spectra_fast(msi_data) do idx, mz, intensity - # Store metadata - modes[idx] = msi_data.spectra_metadata[idx].mode - is_compressed[idx] = msi_data.spectra_metadata[idx].mz_asset.is_compressed || - msi_data.spectra_metadata[idx].int_asset.is_compressed - - if isempty(mz) - empty_spectra += 1 - tics[idx] = 0.0 - bpis[idx] = 0.0 - bp_mzs[idx] = 0.0 - num_points[idx] = 0 - min_mzs[idx] = Inf - max_mzs[idx] = -Inf + lock(msi_data.cache_lock) do + if is_set(msi_data.analytics_ready) return - else - spectra_with_data += 1 + end + + println("Pre-computing analytics (streaming mode)...") + start_time = time_ns() + + num_spectra = length(msi_data.spectra_metadata) + + # Initialize thread-local variables for global stats + thread_local_min_mz = [Inf for _ in 1:Threads.nthreads()] + thread_local_max_mz = [-Inf for _ in 1:Threads.nthreads()] + + # Initialize vectors for per-spectrum stats + tics = Vector{Float64}(undef, num_spectra) + bpis = Vector{Float64}(undef, num_spectra) + bp_mzs = Vector{Float64}(undef, num_spectra) + num_points = Vector{Int}(undef, num_spectra) + min_mzs = Vector{Float64}(undef, num_spectra) + max_mzs = Vector{Float64}(undef, num_spectra) + modes = Vector{SpectrumMode}(undef, num_spectra) + + # Don't precompute Bloom filters by default - create them lazily + bloom_filters = nothing + + # Process in chunks to reduce memory pressure + chunk_size = min(1000, num_spectra) + + for chunk_start in 1:chunk_size:num_spectra + chunk_end = min(chunk_start + chunk_size - 1, num_spectra) + chunk_range = chunk_start:chunk_end + + println("Processing chunk $chunk_start - $chunk_end / $num_spectra") + + # Process current chunk + _iterate_spectra_fast(msi_data, collect(chunk_range)) do idx, mz, intensity + # Store metadata + modes[idx] = msi_data.spectra_metadata[idx].mode + + if isempty(mz) + tics[idx] = 0.0 + bpis[idx] = 0.0 + bp_mzs[idx] = 0.0 + num_points[idx] = 0 + min_mzs[idx] = Inf + max_mzs[idx] = -Inf + return + end + + # Update thread-local global m/z range + local_min, local_max = extrema(mz) + thread_id = Threads.threadid() + thread_local_min_mz[thread_id] = min(thread_local_min_mz[thread_id], local_min) + thread_local_max_mz[thread_id] = max(thread_local_max_mz[thread_id], local_max) + min_mzs[idx] = local_min + max_mzs[idx] = local_max + + # Calculate per-spectrum stats + tics[idx] = sum(intensity) + max_int, max_idx = findmax(intensity) + bpis[idx] = max_int + bp_mzs[idx] = mz[max_idx] + num_points[idx] = length(mz) + end + + # Force GC between chunks to control memory + if chunk_end % 5000 == 0 + GC.gc() + end end - # Update global m/z range - local_min, local_max = extrema(mz) - g_min_mz = min(g_min_mz, local_min) - g_max_mz = max(g_max_mz, local_max) - min_mzs[idx] = local_min - max_mzs[idx] = local_max + # Combine thread-local global stats + g_min_mz = minimum(thread_local_min_mz) + g_max_mz = maximum(thread_local_max_mz) - # Calculate per-spectrum stats - tics[idx] = sum(intensity) - max_int, max_idx = findmax(intensity) - bpis[idx] = max_int - bp_mzs[idx] = mz[max_idx] - num_points[idx] = length(mz) + # Create results + computed_stats_df = DataFrame( + SpectrumID = 1:num_spectra, + TIC = tics, + BPI = bpis, + BasePeakMZ = bp_mzs, + NumPoints = num_points, + MinMZ = min_mzs, + MaxMZ = max_mzs, + Mode = modes + ) + + # Set atomic values + Threads.atomic_xchg!(msi_data.global_min_mz, g_min_mz) + Threads.atomic_xchg!(msi_data.global_max_mz, g_max_mz) + + # Set other data + msi_data.spectrum_stats_df = computed_stats_df + msi_data.bloom_filters = bloom_filters # Will be created on-demand + + set!(msi_data.analytics_ready) + + duration = (time_ns() - start_time) / 1e9 + @printf "Streaming analytics complete in %.2f seconds.\n" duration end - - # Add mode statistics - centroid_count = count(==(CENTROID), modes) - profile_count = count(==(PROFILE), modes) - unknown_count = count(==(UNKNOWN), modes) - - println("DEBUG Mode Statistics:") - println(" Centroid spectra: $centroid_count") - println(" Profile spectra: $profile_count") - println(" Unknown mode: $unknown_count") - - # DEBUG: Print summary - println("DEBUG Analytics Summary:") - println(" Total spectra: $num_spectra") - println(" Spectra with data: $spectra_with_data") - println(" Empty spectra: $empty_spectra") - println(" Global m/z range: [$g_min_mz, $g_max_mz]") - println(" Centroid spectra: $(count(==(CENTROID), modes))") - println(" Profile spectra: $(count(==(PROFILE), modes))") - println(" Compressed spectra: $(sum(is_compressed))") - println(" Average points per spectrum: $(mean(num_points))") - - # Populate the MSIData object - msi_data.global_min_mz = g_min_mz - msi_data.global_max_mz = g_max_mz - msi_data.spectrum_stats_df = DataFrame( - SpectrumID = 1:num_spectra, - TIC = tics, - BPI = bpis, - BasePeakMZ = bp_mzs, - NumPoints = num_points, - MinMZ = min_mzs, - MaxMZ = max_mzs, - Mode = modes, - IsCompressed = is_compressed - ) - - duration = (time_ns() - start_time) / 1e9 - @printf "Analytics pre-computation complete in %.2f seconds.\n" duration return end + + """ get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000) @@ -661,10 +776,10 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_ local global_min_mz, global_max_mz - if msi_data.global_min_mz !== nothing + if Threads.atomic_add!(msi_data.global_min_mz, 0.0) !== Inf println(" Using pre-computed m/z range.") - global_min_mz = msi_data.global_min_mz - global_max_mz = msi_data.global_max_mz + global_min_mz = Threads.atomic_add!(msi_data.global_min_mz, 0.0) + global_max_mz = Threads.atomic_add!(msi_data.global_max_mz, 0.0) else # 1. First Pass: Find the global m/z range by reading the fast .ibd file pass1_start_time = time_ns() @@ -687,24 +802,26 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_ # Use and cache the result global_min_mz = g_min_mz global_max_mz = g_max_mz - msi_data.global_min_mz = global_min_mz - msi_data.global_max_mz = global_max_mz + set_global_mz_range!(msi_data, global_min_mz, global_max_mz) println(" Global m/z range found and cached: [$(global_min_mz), $(global_max_mz)]") @printf " (Pass 1 took %.2f seconds)\n" pass1_duration end # 2. Define Bins and precompute constants mz_bins = range(global_min_mz, stop=global_max_mz, length=num_bins) - intensity_sum = zeros(Float64, num_bins) bin_step = step(mz_bins) inv_bin_step = 1.0 / bin_step # Precompute reciprocal to avoid division min_mz = global_min_mz + # Initialize thread-local intensity sums + thread_local_intensity_sums = [zeros(Float64, num_bins) for _ in 1:Threads.nthreads()] + thread_local_spectra_processed = [0 for _ in 1:Threads.nthreads()] + # 3. Second Pass: Optimized binning pass2_start_time = time_ns() println(" Pass 2: Summing intensities into $num_bins bins...") - num_spectra_processed = 0 _iterate_spectra_fast(msi_data, masked_indices === nothing ? nothing : collect(masked_indices)) do idx, mz, intensity - num_spectra_processed += 1 + thread_id = Threads.threadid() + thread_local_spectra_processed[thread_id] += 1 if isempty(mz) return end @@ -721,9 +838,17 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_ # Branchless version using clamp final_index = clamp(bin_index, 1, num_bins) - intensity_sum[final_index] += intensity_view[i] + thread_local_intensity_sums[thread_id][final_index] += intensity_view[i] end end + + # Combine thread-local results + intensity_sum = zeros(Float64, num_bins) + num_spectra_processed = 0 + for i in 1:Threads.nthreads() + intensity_sum .+= thread_local_intensity_sums[i] + num_spectra_processed += thread_local_spectra_processed[i] + end pass2_duration = (time_ns() - pass2_start_time) / 1e9 total_duration = (time_ns() - total_start_time) / 1e9 @@ -753,10 +878,10 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000, masked_i local global_min_mz, global_max_mz - if msi_data.global_min_mz !== nothing + if Threads.atomic_add!(msi_data.global_min_mz, 0.0) !== Inf println(" Using pre-computed m/z range.") - global_min_mz = msi_data.global_min_mz - global_max_mz = msi_data.global_max_mz + global_min_mz = Threads.atomic_add!(msi_data.global_min_mz, 0.0) + global_max_mz = Threads.atomic_add!(msi_data.global_max_mz, 0.0) else # --- Pass 1: Find m/z range and cache it --- pass1_start_time = time_ns() @@ -780,8 +905,7 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000, masked_i # Use and cache the result global_min_mz = g_min_mz global_max_mz = g_max_mz - msi_data.global_min_mz = global_min_mz - msi_data.global_max_mz = global_max_mz + set_global_mz_range!(msi_data, global_min_mz, global_max_mz) println(" Global m/z range found and cached: [$(global_min_mz), $(global_max_mz)]") @printf " (Pass 1 took %.2f seconds)\n" pass1_duration @@ -792,14 +916,17 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000, masked_i println(" Pass 2: Summing intensities into $num_bins bins...") mz_bins = range(global_min_mz, stop=global_max_mz, length=num_bins) - intensity_sum = zeros(Float64, num_bins) bin_step = step(mz_bins) inv_bin_step = 1.0 / bin_step # Precompute reciprocal min_mz = global_min_mz - num_spectra_processed = 0 + # Initialize thread-local intensity sums + thread_local_intensity_sums = [zeros(Float64, num_bins) for _ in 1:Threads.nthreads()] + thread_local_spectra_processed = [0 for _ in 1:Threads.nthreads()] + _iterate_spectra_fast(msi_data, masked_indices === nothing ? nothing : collect(masked_indices)) do idx, mz, intensity - num_spectra_processed += 1 + thread_id = Threads.threadid() + thread_local_spectra_processed[thread_id] += 1 if isempty(mz) return end @@ -811,9 +938,17 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000, masked_i # Branchless version using clamp final_index = clamp(bin_index, 1, num_bins) - intensity_sum[final_index] += intensity[i] + thread_local_intensity_sums[thread_id][final_index] += intensity[i] end end + + # Combine thread-local results + intensity_sum = zeros(Float64, num_bins) + num_spectra_processed = 0 + for i in 1:Threads.nthreads() + intensity_sum .+= thread_local_intensity_sums[i] + num_spectra_processed += thread_local_spectra_processed[i] + end pass2_duration = (time_ns() - pass2_start_time) / 1e9 total_duration = (time_ns() - total_start_time) / 1e9 @@ -956,7 +1091,7 @@ the file stream `io` is already positioned at the correct offset. - An error if zlib decompression fails, which can indicate corrupt data or an incorrect offset in the `.imzML` metadata. """ -function read_compressed_array(io::IO, asset::SpectrumAsset, format::Type) +function read_compressed_array(data::MSIData, io::IO, asset::SpectrumAsset, format::Type) # Add validation before seeking if asset.offset < 0 || asset.offset >= filesize(io) throw(FileFormatError("Invalid asset offset: $(asset.offset) for file size $(filesize(io))")) @@ -964,30 +1099,43 @@ function read_compressed_array(io::IO, asset::SpectrumAsset, format::Type) seek(io, asset.offset) if asset.is_compressed - # Read compressed bytes - use encoded_length as the number of compressed bytes - compressed_bytes = read(io, asset.encoded_length) + # 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))") + println("DEBUG: Decompressing data - offset=$(asset.offset), compressed_bytes=$(length(compressed_bytes_buffer))") - local decompressed_bytes + local decompressed_bytes_buffer try - decompressed_bytes = Libz.inflate(compressed_bytes) - println("DEBUG: Decompression successful - decompressed_bytes=$(length(decompressed_bytes))") + # 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) + + 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)) + 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, 1:bytes_to_print)) + println(stderr, view(compressed_bytes_buffer, 1:bytes_to_print)) rethrow(e) + finally + release_buffer!(data.buffer_pool, compressed_bytes_buffer) end # Use an IOBuffer to safely read the data - bytes_io = IOBuffer(decompressed_bytes) + bytes_io = IOBuffer(decompressed_bytes_buffer) n_elements = bytes_io.size ÷ sizeof(format) array = Array{format}(undef, n_elements) read!(bytes_io, array) + + release_buffer!(data.buffer_pool, decompressed_bytes_buffer) return array else # Read uncompressed data directly @@ -1082,8 +1230,8 @@ function _iterate_compressed_fast(f::Function, data::MSIData, source::ImzMLSourc end # Read and decompress each array - mz_array = read_compressed_array(source.ibd_handle, meta.mz_asset, source.mz_format) - intensity_array = read_compressed_array(source.ibd_handle, meta.int_asset, source.intensity_format) + 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) mz_array .= ltoh.(mz_array) intensity_array .= ltoh.(intensity_array) @@ -1143,15 +1291,15 @@ function _iterate_spectra_fast_impl(f::Function, data::MSIData, source::MzMLSour for (i, _) in indices_with_offsets meta = data.spectra_metadata[i] - mz = read_binary_vector(source.file_handle, meta.mz_asset) - intensity = read_binary_vector(source.file_handle, meta.int_asset) + mz = read_binary_vector(data, source.file_handle, meta.mz_asset) + intensity = read_binary_vector(data, source.file_handle, meta.int_asset) f(i, mz, intensity) end end """ - _iterate_spectra_fast(f::Function, data::MSIData, indices_to_iterate=nothing) + _iterate_spectra_fast_serial(f::Function, data::MSIData, indices_to_iterate=nothing) An internal, high-performance iterator for bulk processing that bypasses the `GetSpectrum` cache. It provides direct, sequential access to spectral data, making it the most efficient @@ -1173,7 +1321,183 @@ This is a low-level function that bypasses the public `GetSpectrum` API. It does use the cache and is not guaranteed to be thread-safe if the same `MSIData` object is accessed from multiple threads concurrently. """ -function _iterate_spectra_fast(f::Function, data::MSIData, indices_to_iterate::Union{AbstractVector{Int}, Nothing}=nothing) +function _iterate_spectra_fast_serial(f::Function, data::MSIData, indices_to_iterate::Union{AbstractVector{Int}, Nothing}=nothing) # Dispatch to the correct implementation based on the source type _iterate_spectra_fast_impl(f, data, data.source, indices_to_iterate) -end \ No newline at end of file +end + +""" + _iterate_spectra_fast(f::Function, data::MSIData, indices_to_iterate=nothing) + +A smart dispatcher for internal, high-performance iteration over spectra. +It automatically chooses between serial and parallel execution based on the +number of available threads (`Threads.nthreads()`). + +# Arguments +- `f`: A function to call for each spectrum, with the signature `f(index, mz, intensity)`. +- `data`: The `MSIData` object. +- `indices_to_iterate`: An optional `AbstractVector{Int}` specifying a subset of + spectrum indices to iterate over. If `nothing`, it iterates over all spectra. +""" +function _iterate_spectra_fast(f::Function, data::MSIData, indices_to_iterate::Union{AbstractVector{Int}, Nothing}=nothing) + # Default to all indices if not specified + all_indices = (indices_to_iterate === nothing) ? (1:length(data.spectra_metadata)) : indices_to_iterate + + if Threads.nthreads() > 1 && length(all_indices) > 100 # Heuristic: only parallelize for enough work + _iterate_spectra_fast_parallel(f, data, all_indices) + else + _iterate_spectra_fast_serial(f, data, all_indices) + end +end + +""" + _iterate_spectra_fast_threadsafe(f::Function, data::MSIData, indices::AbstractVector{Int}) + +Thread-safe version of the fast iterator that uses the cache lock to ensure exclusive file access. +This is safe to use when multiple threads might access the same MSIData object concurrently. + +# Arguments +- `f`: Function to call for each spectrum with signature `f(index, mz, intensity)` +- `data`: The MSIData object +- `indices`: Specific spectrum indices to iterate over + +# Performance +- Slower than non-thread-safe version due to locking +- Safe for concurrent access to same MSIData object +- Good for small to medium datasets where parallel overhead isn't justified +""" +function _iterate_spectra_fast_threadsafe(f::Function, data::MSIData, indices::AbstractVector{Int}) + # Use the cache lock to ensure exclusive access to file + lock(data.cache_lock) do + _iterate_spectra_fast(f, data, indices) + end +end + +""" + _iterate_spectra_fast_parallel(f::Function, data::MSIData, indices::AbstractVector{Int}) + +Parallel version that creates thread-local file handles for maximum performance. +Each thread gets its own file handle, eliminating contention. + +# Arguments +- `f`: Function to call for each spectrum with signature `f(index, mz, intensity)` +- `data`: The MSIData object +- `indices`: Specific spectrum indices to iterate over + +# Performance +- Fastest for large datasets on multi-core systems +- Eliminates file I/O contention between threads +- Higher memory usage due to multiple file handles +- 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 = Threads.nthreads() + chunk_size = ceil(Int, length(indices) / n_chunks) + chunks = collect(Iterators.partition(indices, chunk_size)) + + 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 + end + end +end + +""" + _iterate_spectra_fast_parallel(f::Function, data::MSIData) + +Parallel version that processes all spectra in the dataset. +Convenience wrapper for the indexed version. +""" +function _iterate_spectra_fast_parallel(f::Function, data::MSIData) + all_indices = collect(1:length(data.spectra_metadata)) + _iterate_spectra_fast_parallel(f, data, all_indices) +end + +""" + _iterate_spectra_fast_threadsafe(f::Function, data::MSIData) + +Thread-safe version that processes all spectra in the dataset. +Convenience wrapper for the indexed version. +""" +function _iterate_spectra_fast_threadsafe(f::Function, data::MSIData) + all_indices = collect(1:length(data.spectra_metadata)) + _iterate_spectra_fast_threadsafe(f, data, all_indices) +end + +""" + get_bloom_filter(data::MSIData, index::Int)::BloomFilter{Int} + +Retrieves the Bloom filter for a specific spectrum by its index. +If the Bloom filter has not yet been created (due to lazy initialization), +it will be created on-demand and cached. + +# Arguments +- `data`: The `MSIData` object. +- `index`: The linear index of the spectrum. + +# Returns +- A `BloomFilter{Int}` for the specified spectrum. +""" +function get_bloom_filter(data::MSIData, index::Int)::BloomFilter{Int} + lock(data.cache_lock) do + if data.bloom_filters === nothing + # Initialize Bloom filters on-demand + data.bloom_filters = Vector{Union{BloomFilter{Int}, Nothing}}(undef, length(data.spectra_metadata)) + fill!(data.bloom_filters, nothing) # Initialize with nothing + end + + if data.bloom_filters[index] === nothing + # Create Bloom filter lazily + mz, intensity = GetSpectrum(data, index) + data.bloom_filters[index] = create_bloom_filter_for_spectrum(mz) + end + + return data.bloom_filters[index] + end +end + +""" + monitor_memory_usage(phase::String) + +Utility function to print current memory usage (allocated bytes and RSS) +at different phases of execution. + +# Arguments +- `phase`: A string describing the current phase (e.g., "Initialization", "Processing chunk 1"). +""" +function monitor_memory_usage(phase::String) + # Use Julia's internal memory tracking + bytes_allocated = Base.gc_bytes() + mb_allocated = bytes_allocated / 1024^2 + + # Get process memory from system + if Sys.islinux() + proc_status = read("/proc/self/status", String) + vm_rss_match = match(r"VmRSS:\s*(\d+)\s*kB", proc_status) + if vm_rss_match !== nothing + vm_rss_mb = parse(Int, vm_rss_match[1]) / 1024 + @printf "[%s] Memory: %.2f MB allocated, %.2f MB RSS\n" phase mb_allocated vm_rss_mb + end + else + @printf "[%s] Memory: %.2f MB allocated\n" phase mb_allocated + end +end diff --git a/src/MSI_src.jl b/src/MSI_src.jl index f855f16..bf5df59 100644 --- a/src/MSI_src.jl +++ b/src/MSI_src.jl @@ -6,18 +6,17 @@ module MSI_src export OpenMSIData, GetSpectrum, IterateSpectra, - ImportMzmlFile, - load_slices, + ImportMzmlFile, plot_slices, plot_slice, get_total_spectrum, - get_average_spectrum, - LoadMzml, + get_average_spectrum, precompute_analytics, process_spectrum, generate_colorbar_image, process_image_pipeline, - load_and_prepare_mask + load_and_prepare_mask, + set_global_mz_range! # Export the public Preprocessing API export FeatureMatrix, @@ -40,6 +39,7 @@ export FeatureMatrix, get_common_calibration_standards # Include all source files directly into the main module +include("BloomFilters.jl") include("Common.jl") include("MSIData.jl") include("ParserHelpers.jl") diff --git a/src/imzML.jl b/src/imzML.jl index 10dccc4..6d66485 100644 --- a/src/imzML.jl +++ b/src/imzML.jl @@ -105,45 +105,6 @@ end Reads metadata to determine the byte offsets and data types for reading spectra. """ function get_spectrum_attributes(stream::IO, hIbd::IO) - skip = Vector{UInt32}(undef, 8) - offset = get_spectrum_tag_offset(stream) - - tag = find_tag(stream, r" accession=\"IMS:100005(\d)\"(.+)") - skip[1] = tag.captures[1][1] - '0' + 1 - skip[2] = xor(skip[1], 3) - - value = match(r" value=\"(\d+)\".+", tag.captures[2]) - skip[5] = position(stream) - offset - length(value.match) - 2 - value = match(r" value=\"\d+\".+", readline(stream)) - skip[6] = value.offset - - offset = position(stream) - tag = find_tag(stream, r"^\s*") - skip[8] = position(stream) - offset - return skip -end - -function get_spectrum_attributes_v2(stream::IO, hIbd::IO) # Look for position x readuntil(stream, "IMS:1000050") x_skip = 0 @@ -239,7 +200,7 @@ 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") - hIbd = open(ibd_path, "r") + ts_hIbd = ThreadSafeFileHandle(ibd_path) try println("DEBUG: Configuring axes...") @@ -287,7 +248,7 @@ function load_imzml_lazy(file_path::String; cache_size::Int=100) local spectra_metadata if parser_type == :neofx println("DEBUG: Using neofx parser.") - spectra_metadata = parse_neofx(stream, hIbd, param_groups, width, height, num_spectra, + spectra_metadata = parse_neofx(stream, ts_hIbd, param_groups, width, height, num_spectra, default_mz_format, default_intensity_format, mz_is_compressed, int_is_compressed, global_mode) #= @@ -305,7 +266,7 @@ function load_imzml_lazy(file_path::String; cache_size::Int=100) =# else println("DEBUG: Using compressed parser.") - spectra_metadata = parse_compressed(stream, hIbd, param_groups, width, height, num_spectra, + spectra_metadata = parse_compressed(stream, ts_hIbd, param_groups, width, height, num_spectra, default_mz_format, default_intensity_format, mz_is_compressed, int_is_compressed, global_mode) end @@ -325,19 +286,23 @@ function load_imzml_lazy(file_path::String; cache_size::Int=100) end println("DEBUG: Coordinate map built.") + source = ImzMLSource(ts_hIbd, default_mz_format, default_intensity_format) + println("DEBUG: Creating MSIData object.") + msi_data = MSIData(source, spectra_metadata, (width, height), coordinate_map, cache_size) + + # Close the XML stream as it's no longer needed close(stream) - source = ImzMLSource(hIbd, default_mz_format, default_intensity_format) - println("DEBUG: Creating MSIData object.") - return MSIData(source, spectra_metadata, (width, height), coordinate_map, cache_size) + return msi_data catch e close(stream) - close(hIbd) + close(ts_hIbd) # Ensure IBD handle is closed on error rethrow(e) end end +#= function parse_uncompressed(stream::IO, hIbd::IO, param_groups::Dict{String, SpecDim}, width::Int32, height::Int32, num_spectra::Int32, mz_format::DataType, intensity_format::DataType, @@ -350,7 +315,7 @@ function parse_uncompressed(stream::IO, hIbd::IO, param_groups::Dict{String, Spe seek(hIbd, 16) current_ibd_offset = position(hIbd) - attr = get_spectrum_attributes_v2(stream, hIbd) + attr = get_spectrum_attributes(stream, hIbd) seek(stream, start_of_spectra_xml) println("DEBUG: Initial IBD offset (after UUID): $current_ibd_offset") @@ -363,7 +328,7 @@ function parse_uncompressed(stream::IO, hIbd::IO, param_groups::Dict{String, Spe # Find X coordinate line = readuntil(stream, "IMS:1000050") if eof(stream) - error("Unexpected EOF while looking for X coordinate") + throw(FileFormatError("Unexpected EOF while looking for X coordinate")) end val_tag_x = find_tag(stream, r"value=\"(\d+)\"") x = parse(Int32, val_tag_x.captures[1]) @@ -428,128 +393,139 @@ function parse_uncompressed(stream::IO, hIbd::IO, param_groups::Dict{String, Spe return spectra_metadata end +=# -function parse_compressed(stream::IO, hIbd::IO, param_groups::Dict{String, SpecDim}, +function parse_compressed(stream::IO, hIbd::Union{IO, ThreadSafeFileHandle}, param_groups::Dict{String, SpecDim}, width::Int32, height::Int32, num_spectra::Int32, default_mz_format::DataType, default_intensity_format::DataType, mz_is_compressed::Bool, int_is_compressed::Bool, global_mode::SpectrumMode) spectra_metadata = Vector{SpectrumMetadata}(undef, num_spectra) - # Use concrete type for array_data array_data_type = @NamedTuple{is_mz::Bool, array_length::Int32, encoded_length::Int64, offset::Int64} for k in 1:num_spectra - # Read the full spectrum XML block - spectrum_buffer = IOBuffer() + # Initialize variables for this spectrum + x = Int32(0) + y = Int32(0) + spectrum_mode = global_mode + current_array_data = array_data_type[] + spectrum_start_line = "" + + # Find the start of the spectrum tag line = "" while !eof(stream) line = readline(stream) if occursin("", line) break end - end - spectrum_xml = String(take!(spectrum_buffer)) - # Parse coordinates - x_match = match(r"IMS:1000050.*?value=\"(\d+)\"", spectrum_xml) - y_match = match(r"IMS:1000051.*?value=\"(\d+)\"", spectrum_xml) - x = x_match !== nothing ? parse(Int32, x_match.captures[1]) : Int32(0) - y = y_match !== nothing ? parse(Int32, y_match.captures[1]) : Int32(0) - - # Parse mode - use SpectrumMode type directly - spectrum_mode = global_mode - if occursin("MS:1000127", spectrum_xml) - spectrum_mode = CENTROID - elseif occursin("MS:1000128", spectrum_xml) - spectrum_mode = PROFILE - end - - # Parse binary data arrays with concrete typing - array_data = array_data_type[] - - # Find all binaryDataArray blocks - array_matches = eachmatch(r""s, spectrum_xml) - for array_match in array_matches - array_xml = array_match.match - - # Determine if this is m/z or intensity array - is_mz = occursin("MS:1000514", array_xml) || occursin("mzArray", array_xml) - - # Parse external data parameters - # Get array_length (nPoints) - array_len_cv_match = match(r"IMS:1000103.*?value=\"(\d+)\"", array_xml) - array_length = Int32(0) - if array_len_cv_match !== nothing - array_length = parse(Int32, array_len_cv_match.captures[1]) + # Parse coordinates + x_match = match(r"IMS:1000050.*?value=\"(\d+)\"", line) + if x_match !== nothing + x = parse(Int32, x_match.captures[1]) end - if array_length == 0 - nPoints_match = match(r"defaultArrayLength=\"(\d+)\"", spectrum_xml) - if nPoints_match !== nothing - array_length = parse(Int32, nPoints_match.captures[1]) + y_match = match(r"IMS:1000051.*?value=\"(\d+)\"", line) + if y_match !== nothing + y = parse(Int32, y_match.captures[1]) + end + + # Parse mode + if occursin("MS:1000127", line) + spectrum_mode = CENTROID + elseif occursin("MS:1000128", line) + spectrum_mode = PROFILE + end + + # More robust binaryDataArray detection + if occursin("", line) + while !eof(stream) + bda_line = readline(stream) + push!(bda_lines, bda_line) + if occursin("", bda_line) + break + end + end + end + bda_content = join(bda_lines, "\n") + + # Parse from bda_content + is_mz = occursin("MS:1000514", bda_content) || occursin("mzArray", bda_content) + + array_len_cv_match = match(r"IMS:1000103.*?value=\"(\d+)\"", bda_content) + array_length = Int32(0) + if array_len_cv_match !== nothing + array_length = parse(Int32, array_len_cv_match.captures[1]) + end + + # Fallback for defaultArrayLength + if array_length == 0 + default_match = match(r"defaultArrayLength=\"(\d+)\"", spectrum_start_line) + if default_match !== nothing + array_length = parse(Int32, default_match.captures[1]) + end + end + + encoded_len_cv_match = match(r"IMS:1000104.*?value=\"(\d+)\"", bda_content) + encoded_length = Int64(0) + if encoded_len_cv_match !== nothing + encoded_length = parse(Int64, encoded_len_cv_match.captures[1]) + else + encoded_len_attr_match = match(r"encodedLength=\"(\d+)\"", bda_content) + if encoded_len_attr_match !== nothing + encoded_length = parse(Int64, encoded_len_attr_match.captures[1]) + end + end + + offset_match = match(r"IMS:1000102.*?value=\"(\d+)\"", bda_content) + offset = Int64(0) + if offset_match !== nothing + offset = parse(Int64, offset_match.captures[1]) + end + + if array_length > 0 && offset > 0 + push!(current_array_data, (is_mz=is_mz, array_length=array_length, + encoded_length=encoded_length, offset=offset)) end end - - # Get encoded_length - encoded_len_cv_match = match(r"IMS:1000104.*?value=\"(\d+)\"", array_xml) - encoded_length = Int64(0) - if encoded_len_cv_match !== nothing - encoded_length = parse(Int64, encoded_len_cv_match.captures[1]) - else - encoded_len_attr_match = match(r"encodedLength=\"(\d+)\"", array_xml) - if encoded_len_attr_match !== nothing - encoded_length = parse(Int64, encoded_len_attr_match.captures[1]) - end - end - - # Get offset - offset_match = match(r"IMS:1000102.*?value=\"(\d+)\"", array_xml) - offset = Int64(0) - if offset_match !== nothing - offset = parse(Int64, offset_match.captures[1]) - end - if array_length > 0 && offset > 0 - push!(array_data, (is_mz=is_mz, array_length=array_length, - encoded_length=encoded_length, offset=offset)) - end + # Reset line to continue loop + line = "" end # Separate m/z and intensity arrays - mz_data = filter(d -> d.is_mz, array_data) - int_data = filter(d -> !d.is_mz, array_data) - - # FIX: Handle empty spectra gracefully instead of throwing errors - current_ibd_pos = position(hIbd) + mz_data = filter(d -> d.is_mz, current_array_data) + int_data = filter(d -> !d.is_mz, current_array_data) if length(mz_data) != 1 || length(int_data) != 1 - # Create placeholder metadata for empty/invalid spectrum println("DEBUG: Spectrum $k is empty or invalid - creating placeholder metadata") - mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, current_ibd_pos, 0, :mz) - int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, current_ibd_pos, 0, :intensity) + 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) else mz_info = mz_data[1] int_info = int_data[1] - # DEBUG: Print first spectrum details if k == 1 println("DEBUG First spectrum parsed:") println(" Coordinates: x=$x, y=$y") println(" Mode: $spectrum_mode") println(" m/z array: array_length=$(mz_info.array_length), encoded_length=$(mz_info.encoded_length), offset=$(mz_info.offset)") println(" intensity array: array_length=$(int_info.array_length), encoded_length=$(int_info.encoded_length), offset=$(int_info.offset)") - println(" Expected m/z bytes: $(mz_info.array_length * sizeof(default_mz_format))") - println(" Expected intensity bytes: $(int_info.array_length * sizeof(default_intensity_format))") end - # Create SpectrumAsset objects mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, mz_info.offset, mz_is_compressed ? mz_info.encoded_length : mz_info.array_length, :mz) int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, int_info.offset, @@ -562,126 +538,137 @@ function parse_compressed(stream::IO, hIbd::IO, param_groups::Dict{String, SpecD return spectra_metadata end -function parse_neofx(stream::IO, hIbd::IO, param_groups::Dict{String, SpecDim}, +function parse_neofx(stream::IO, hIbd::Union{IO, ThreadSafeFileHandle}, param_groups::Dict{String, SpecDim}, width::Int32, height::Int32, num_spectra::Int32, default_mz_format::DataType, default_intensity_format::DataType, mz_is_compressed::Bool, int_is_compressed::Bool, global_mode::SpectrumMode) spectra_metadata = Vector{SpectrumMetadata}(undef, num_spectra) - + array_data_type = @NamedTuple{is_mz::Bool, array_length::Int32, encoded_length::Int64, offset::Int64} for k in 1:num_spectra - # Read the full spectrum XML block - spectrum_buffer = IOBuffer() + # Initialize variables for this spectrum + x = Int32(0) + y = Int32(0) + spectrum_mode = global_mode + current_array_data = array_data_type[] + spectrum_start_line = "" + + # Find the start of the spectrum tag line = "" while !eof(stream) line = readline(stream) if occursin("", line) break end - end - spectrum_xml = String(take!(spectrum_buffer)) - # Parse coordinates - x_match = match(r"IMS:1000050.*?value=\"(\d+)\"", spectrum_xml) - y_match = match(r"IMS:1000051.*?value=\"(\d+)\"", spectrum_xml) - x = x_match !== nothing ? parse(Int32, x_match.captures[1]) : Int32(0) - y = y_match !== nothing ? parse(Int32, y_match.captures[1]) : Int32(0) - - # Parse mode - spectrum_mode = global_mode - if occursin("MS:1000127", spectrum_xml) - spectrum_mode = CENTROID - elseif occursin("MS:1000128", spectrum_xml) - spectrum_mode = PROFILE - end - - # Parse binary data arrays - array_data = array_data_type[] - - # Find all binaryDataArray blocks - array_matches = eachmatch(r""s, spectrum_xml) - for array_match in array_matches - array_xml = array_match.match - - # Determine if this is m/z or intensity array - is_mz = occursin("MS:1000514", array_xml) || occursin("mzArray", array_xml) - - # Parse external data parameters - # Get array_length (nPoints) - array_len_cv_match = match(r"IMS:1000103.*?value=\"(\d+)\"", array_xml) - array_length = Int32(0) - if array_len_cv_match !== nothing - array_length = parse(Int32, array_len_cv_match.captures[1]) + # Parse coordinates + x_match = match(r"IMS:1000050.*?value=\"(\d+)\"", line) + if x_match !== nothing + x = parse(Int32, x_match.captures[1]) end - if array_length == 0 - nPoints_match = match(r"defaultArrayLength=\"(\d+)\"", spectrum_xml) - if nPoints_match !== nothing - array_length = parse(Int32, nPoints_match.captures[1]) + y_match = match(r"IMS:1000051.*?value=\"(\d+)\"", line) + if y_match !== nothing + y = parse(Int32, y_match.captures[1]) + end + + # Parse mode + if occursin("MS:1000127", line) + spectrum_mode = CENTROID + elseif occursin("MS:1000128", line) + spectrum_mode = PROFILE + end + + # More robust binaryDataArray detection + if occursin("", line) + while !eof(stream) + bda_line = readline(stream) + push!(bda_lines, bda_line) + if occursin("", bda_line) + break + end + end + end + bda_content = join(bda_lines, "\n") + + # Parse from bda_content + is_mz = occursin("MS:1000514", bda_content) || occursin("mzArray", bda_content) + + array_len_cv_match = match(r"IMS:1000103.*?value=\"(\d+)\"", bda_content) + array_length = Int32(0) + if array_len_cv_match !== nothing + array_length = parse(Int32, array_len_cv_match.captures[1]) + end + + # Fallback for defaultArrayLength + if array_length == 0 + default_match = match(r"defaultArrayLength=\"(\d+)\"", spectrum_start_line) + if default_match !== nothing + array_length = parse(Int32, default_match.captures[1]) + end + end + + encoded_len_cv_match = match(r"IMS:1000104.*?value=\"(\d+)\"", bda_content) + encoded_length = Int64(0) + if encoded_len_cv_match !== nothing + encoded_length = parse(Int64, encoded_len_cv_match.captures[1]) + else + encoded_len_attr_match = match(r"encodedLength=\"(\d+)\"", bda_content) + if encoded_len_attr_match !== nothing + encoded_length = parse(Int64, encoded_len_attr_match.captures[1]) + end + end + + offset_match = match(r"IMS:1000102.*?value=\"(\d+)\"", bda_content) + offset = Int64(0) + if offset_match !== nothing + offset = parse(Int64, offset_match.captures[1]) + end + + if array_length > 0 && offset > 0 + push!(current_array_data, (is_mz=is_mz, array_length=array_length, + encoded_length=encoded_length, offset=offset)) end end - - # Get encoded_length - encoded_len_cv_match = match(r"IMS:1000104.*?value=\"(\d+)\"", array_xml) - encoded_length = Int64(0) - if encoded_len_cv_match !== nothing - encoded_length = parse(Int64, encoded_len_cv_match.captures[1]) - else - encoded_len_attr_match = match(r"encodedLength=\"(\d+)\"", array_xml) - if encoded_len_attr_match !== nothing - encoded_length = parse(Int64, encoded_len_attr_match.captures[1]) - end - end - - # Get offset - offset_match = match(r"IMS:1000102.*?value=\"(\d+)\"", array_xml) - offset = Int64(0) - if offset_match !== nothing - offset = parse(Int64, offset_match.captures[1]) - end - if array_length > 0 && offset > 0 - push!(array_data, (is_mz=is_mz, array_length=array_length, - encoded_length=encoded_length, offset=offset)) - end + # Reset line to continue loop + line = "" end # Separate m/z and intensity arrays - mz_data = filter(d -> d.is_mz, array_data) - int_data = filter(d -> !d.is_mz, array_data) - - # FIX: Handle empty spectra gracefully instead of throwing errors - current_ibd_pos = position(hIbd) + mz_data = filter(d -> d.is_mz, current_array_data) + int_data = filter(d -> !d.is_mz, current_array_data) if length(mz_data) != 1 || length(int_data) != 1 - # Create placeholder metadata for empty/invalid spectrum println("DEBUG: Spectrum $k is empty or invalid - creating placeholder metadata") - mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, current_ibd_pos, 0, :mz) - int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, current_ibd_pos, 0, :intensity) + 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) else mz_info = mz_data[1] int_info = int_data[1] - # DEBUG: Print first spectrum details if k == 1 println("DEBUG First spectrum parsed:") println(" Coordinates: x=$x, y=$y") println(" Mode: $spectrum_mode") println(" m/z array: array_length=$(mz_info.array_length), encoded_length=$(mz_info.encoded_length), offset=$(mz_info.offset)") println(" intensity array: array_length=$(int_info.array_length), encoded_length=$(int_info.encoded_length), offset=$(int_info.offset)") - println(" Expected m/z bytes: $(mz_info.array_length * sizeof(default_mz_format))") - println(" Expected intensity bytes: $(int_info.array_length * sizeof(default_intensity_format))") end - # Create SpectrumAsset objects mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, mz_info.offset, mz_is_compressed ? mz_info.encoded_length : mz_info.array_length, :mz) int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, int_info.offset, @@ -739,6 +726,7 @@ function find_mass(mz_array::AbstractVector{<:Real}, intensity_array::AbstractVe return max_intensity end +#= """ load_slices(folder, masses, tolerance) @@ -795,7 +783,7 @@ function load_slices(folder::String, masses::AbstractVector{<:Real}, tolerance:: return (img_list, names) end - +=# """ get_mz_slice(data::MSIData, mass::Real, tolerance::Real) @@ -818,7 +806,7 @@ function get_mz_slice(data::MSIData, mass::Real, tolerance::Real; mask_path::Uni end # INTELLIGENT LOADING: Ensure analytics are computed for filtering. - if data.spectrum_stats_df === nothing || !hasproperty(data.spectrum_stats_df, :MinMZ) + if !is_set(data.analytics_ready) println("Per-spectrum metadata not found. Running one-time analytics computation...") precompute_analytics(data) end @@ -826,13 +814,33 @@ function get_mz_slice(data::MSIData, mass::Real, tolerance::Real; mask_path::Uni println("Using high-performance sequential iterator...") target_min = mass - tolerance target_max = mass + tolerance - stats_df = data.spectrum_stats_df + stats_df = get_spectrum_stats(data) + bloom_filters = get_bloom_filters(data) # 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 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 + end + end + + if !found + continue # Definitely not in this spectrum + 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 @@ -890,13 +898,14 @@ function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, t end # 2. Ensure analytics are computed for filtering. - if data.spectrum_stats_df === nothing || !hasproperty(data.spectrum_stats_df, :MinMZ) + if !is_set(data.analytics_ready) println("Per-spectrum metadata not found. Running one-time analytics computation...") precompute_analytics(data) end println("Filtering candidate spectra for $(length(masses)) m/z values...") - stats_df = data.spectrum_stats_df + stats_df = get_spectrum_stats(data) + bloom_filters = get_bloom_filters(data) candidate_indices = Set{Int}() indices_to_check = masked_indices === nothing ? (1:length(data.spectra_metadata)) : masked_indices @@ -909,6 +918,24 @@ function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, t if i in candidate_indices continue end + # 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 + end + end + + if !found + continue # Definitely not in this spectrum + 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 @@ -1272,7 +1299,7 @@ end # Analysis and Visualization # # ============================================================================ - +#= """ display_statistics(slices::AbstractArray{<:AbstractMatrix{<:Real}, 2}, names::AbstractVector{String}, masses::AbstractVector{<:Real}) @@ -1314,6 +1341,7 @@ function display_statistics(slices::AbstractArray{<:AbstractMatrix{<:Real}, 2}, end return all_dfs end +=# """ plot_slices(slices, names, masses, output_dir; stage_name, bins=256, dpi=150, global_bounds=nothing) @@ -1537,47 +1565,6 @@ function generate_palette(colorscheme, n_colors=256) return palette end -# ============================================================================ -# -# -# Additional Performance Optimizations -# -# ============================================================================ - -""" - precompute_spectrum_bounds(msi_data::MSIData) - -Precomputes min/max m/z values for each spectrum to enable fast filtering. -This avoids repeated computation during slice extraction. -""" -function precompute_spectrum_bounds(msi_data::MSIData) - n_spectra = length(msi_data.spectra_metadata) - min_mz = Vector{Float64}(undef, n_spectra) - max_mz = Vector{Float64}(undef, n_spectra) - - _iterate_spectra_fast(msi_data) do idx, mz_array, intensity_array - min_mz[idx] = first(mz_array) - max_mz[idx] = last(mz_array) - end - - return (min_mz, max_mz) -end - -""" - process_columns_first!(A::AbstractMatrix) - -Process a matrix column-by-column for better cache performance. -In Julia, arrays are stored in column-major order. -""" -function process_columns_first!(A::AbstractMatrix) - nrows, ncols = size(A) - @inbounds for col in 1:ncols, row in 1:nrows - # Process A[row, col] here - # This order is ~70% faster due to cache locality - end - return A -end - # Viridis color palette (256 colors) - defined as constant const ViridisPalette = generate_palette(ColorSchemes.viridis) diff --git a/src/mzML.jl b/src/mzML.jl index 1709780..1550a08 100644 --- a/src/mzML.jl +++ b/src/mzML.jl @@ -228,22 +228,22 @@ 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") - stream = open(file_path, "r") + ts_stream = ThreadSafeFileHandle(file_path, "r") try println("DEBUG: Finding index offset...") - index_offset = find_index_offset(stream) + index_offset = find_index_offset(ts_stream.handle) println("DEBUG: Seeking to index list at offset $index_offset.") - seek(stream, index_offset) + seek(ts_stream.handle, index_offset) println("DEBUG: Searching for ''.") - if find_tag(stream, r"(\d+)" ) +=# diff --git a/start_MSI_GUI.jl b/start_MSI_GUI.jl index f8211fa..bb7be41 100644 --- a/start_MSI_GUI.jl +++ b/start_MSI_GUI.jl @@ -4,6 +4,7 @@ Pkg.activate(".") # Only instantiate in development mode if get(ENV, "GENIE_ENV", "dev") != "prod" @info "Development environment detected. Instantiating packages..." + Pkg.resolve() Pkg.instantiate() end Pkg.gc()