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