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