Compare commits
No commits in common. "main" and "dev/Adjustments" have entirely different histories.
main
...
dev/Adjust
32
.github/workflows/ci.yml
vendored
32
.github/workflows/ci.yml
vendored
@ -1,32 +0,0 @@
|
||||
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
|
||||
@ -1,31 +0,0 @@
|
||||
# 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"]
|
||||
10
Project.toml
10
Project.toml
@ -1,5 +1,4 @@
|
||||
name = "MSI_src"
|
||||
uuid = "8f395f85-3b9c-4f7f-bf7e-12345678abcd"
|
||||
authors = ["JJSA"]
|
||||
version = "0.1.0"
|
||||
|
||||
@ -44,11 +43,9 @@ 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]
|
||||
@ -94,12 +91,5 @@ 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"]
|
||||
|
||||
45
README.md
45
README.md
@ -35,44 +35,35 @@ Minimum system requirements: 4 core processor, 8 GB RAM<br>
|
||||
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.
|
||||
Alternatively, you can generate the .so file by running the next build script in your directory:
|
||||
|
||||
### 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.*
|
||||
This may take 5–10 minutes as it merges all dependencies and JuliaMSI functions into a single binary. The resulting .so/.dll file will be around **300MB–600MB** because it contains the pre-compiled machine code for your entire environment.
|
||||
A sysimage built on Linux (.so) will not work on Windows.
|
||||
You must run the build_sysimage.jl script once on each target operating system.
|
||||
|
||||
Once `sys_msi_gui.so` is created in your directory, adapt your launch command:
|
||||
Once MSI_sysimage.so is created in your directory, adapt your command like this:
|
||||
```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
|
||||
# This will find MSI_sysimage.so, .dll, or .dylib automatically:
|
||||
julia --threads auto --project=. --sysimage MSI_sysimage* start_MSI_GUI.jl
|
||||
```
|
||||
The command above loads a pre-compiled version of your environment (note, you have to use either the sysimage command or the normal start command), speeding up the booting delay. Your GUI or scripts using JuliaMSI should start almost instantly and process files significantly faster than before.
|
||||
Always run the build script **in the same project root** where you intend to run the app.
|
||||
|
||||
### 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:
|
||||
In scripts like test/run_tests.jl, you no longer need to manually include("../src/MSI_src.jl"). You can simply treat it like a globally installed package to load the precompiled binary:
|
||||
|
||||
```julia
|
||||
using MSI_src
|
||||
# ... rest of your script
|
||||
```
|
||||
### Run with the --sysimage flag
|
||||
Launch your tests using the same flag to bypass all compilation overhead:
|
||||
|
||||
```bash
|
||||
julia --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
|
||||
julia --threads auto --project=. --sysimage MSI_sysimage.so test/run_tests.jl
|
||||
```
|
||||
|
||||
## License
|
||||
|
||||
163
app.jl
163
app.jl
@ -29,9 +29,6 @@ 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()
|
||||
@ -63,7 +60,7 @@ if get(ENV, "GENIE_ENV", "dev") != "prod"
|
||||
println("--- MEMORY LOG [$(context)] ---")
|
||||
println(" Timestamp: $(now())")
|
||||
println(" Process RSS: $(rss_mb) MB")
|
||||
println(" global_msi_data[] size: $(msi_data_size_mb) MB")
|
||||
println(" msi_data size: $(msi_data_size_mb) MB")
|
||||
println(" Cumulative GC time: $(gc_time_s) s")
|
||||
println("--------------------------")
|
||||
end
|
||||
@ -196,15 +193,6 @@ 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 ==
|
||||
@ -488,7 +476,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
|
||||
# == DATA MANAGEMENT VARIABLES ==
|
||||
# Centralized MSIData object
|
||||
# global_msi_data[] is now global to avoid Genie Session memory leaks
|
||||
@out msi_data::Union{MSIData, Nothing} = nothing
|
||||
|
||||
# Image file management
|
||||
@out text_nmass="" # For specific mass charge image creation
|
||||
@ -673,7 +661,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
|
||||
try
|
||||
# 1. Clear large data objects explicitly
|
||||
global_msi_data[] = nothing
|
||||
msi_data = nothing
|
||||
feature_matrix_result = nothing
|
||||
bin_info_result = nothing
|
||||
|
||||
@ -729,41 +717,29 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
sure the file can be processed by later steps like mainProcess
|
||||
=#
|
||||
@onbutton btnSearch begin
|
||||
# 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
|
||||
is_processing = true
|
||||
push!(__model__)
|
||||
picked_route = pick_file(; filterlist="imzML,imzml,mzML,mzml")
|
||||
|
||||
if isnothing(picked_route) || isempty(picked_route)
|
||||
is_processing = false
|
||||
return
|
||||
end
|
||||
|
||||
# 2. Update reactive state synchronously
|
||||
is_processing = true
|
||||
msg = "Opening file: $(basename(picked_route))..."
|
||||
SpectraEnabled = false
|
||||
btnMetadataDisable = true
|
||||
push!(__model__)
|
||||
|
||||
# 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
|
||||
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
|
||||
|
||||
msg = "Opening file: $(basename(picked_route))..."
|
||||
|
||||
try
|
||||
dataset_name = replace(basename(picked_route), r"(\.(imzML|imzml|mzML|mzml))$"i => "")
|
||||
registry = load_registry(registry_path)
|
||||
existing_entry = get(registry, dataset_name, nothing)
|
||||
@ -781,8 +757,8 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
dims = parse.(Int, split(dims_str, " x "))
|
||||
imgWidth, imgHeight = dims[1], dims[2]
|
||||
|
||||
global_msi_data[] = nothing # Ensure data is not held in memory
|
||||
log_memory_usage("Fast Load (global_msi_data[] cleared)", global_msi_data[])
|
||||
msi_data = nothing # Ensure data is not held in memory
|
||||
log_memory_usage("Fast Load (msi_data cleared)", msi_data)
|
||||
btnMetadataDisable = false
|
||||
SpectraEnabled = true
|
||||
selected_folder_main = dataset_name
|
||||
@ -1024,10 +1000,10 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
image_available_folders = deepcopy(img_folders)
|
||||
|
||||
selected_folder_main = dataset_name
|
||||
global_msi_data[] = loaded_data
|
||||
msi_data = loaded_data
|
||||
|
||||
# Determine plot mode from loaded data
|
||||
df = global_msi_data[].spectrum_stats_df
|
||||
df = 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)
|
||||
@ -1037,7 +1013,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
last_plot_mode = "lines" # Default
|
||||
end
|
||||
|
||||
log_memory_usage("Full Load", global_msi_data[])
|
||||
log_memory_usage("Full Load", msi_data)
|
||||
|
||||
eTime = round(time() - sTime, digits=3)
|
||||
msg = "Active file loaded in $(eTime) seconds. Dataset '$(dataset_name)' is ready for analysis."
|
||||
@ -1045,21 +1021,18 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
SpectraEnabled = true
|
||||
|
||||
catch e
|
||||
global_msi_data[] = nothing
|
||||
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
|
||||
end
|
||||
end
|
||||
|
||||
@ -1259,11 +1232,6 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
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
|
||||
@ -1291,9 +1259,9 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
end
|
||||
target_path = entry["source_path"]
|
||||
|
||||
# Ensure global_msi_data[] is for the currently selected file and load if needed
|
||||
# Ensure 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 global_msi_data[] used for plotting/interactive exploration.
|
||||
# with the global msi_data used for plotting/interactive exploration.
|
||||
println("DEBUG: Opening isolated MSIData instance for pipeline stability...")
|
||||
pipeline_msi_data = OpenMSIData(target_path)
|
||||
|
||||
@ -1718,9 +1686,9 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
|
||||
# Determine plot mode for this specific spectrum
|
||||
spectrum_mode_for_plot = "lines" # Default to lines
|
||||
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 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 mode == MSI_src.CENTROID
|
||||
spectrum_mode_for_plot = "stem"
|
||||
end
|
||||
@ -1971,7 +1939,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
@onbutton recalculate_suggestions_btn begin
|
||||
is_processing = true
|
||||
push!(__model__)
|
||||
if global_msi_data[] === nothing
|
||||
if msi_data === nothing
|
||||
msg = "Please load a file first."
|
||||
warning_msg = true
|
||||
return
|
||||
@ -1985,7 +1953,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
for p in reference_peaks_list
|
||||
if tryparse(Float64, string(p["mz"])) !== nothing
|
||||
)
|
||||
recommended_params = main_precalculation(global_msi_data[], reference_peaks=ref_peaks)
|
||||
recommended_params = main_precalculation(msi_data, reference_peaks=ref_peaks)
|
||||
|
||||
|
||||
for (step_name, params) in recommended_params
|
||||
@ -2380,11 +2348,6 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
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..."
|
||||
@ -2609,21 +2572,21 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
return
|
||||
end
|
||||
|
||||
if global_msi_data[] === nothing || full_route != target_path
|
||||
if global_msi_data[] !== nothing
|
||||
close(global_msi_data[])
|
||||
if msi_data === nothing || full_route != target_path
|
||||
if msi_data !== nothing
|
||||
close(msi_data)
|
||||
end
|
||||
msg = "Reloading $(basename(target_path)) for analysis..."
|
||||
full_route = target_path
|
||||
global_msi_data[] = OpenMSIData(target_path)
|
||||
msi_data = OpenMSIData(target_path)
|
||||
if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing
|
||||
raw_min = entry["metadata"]["global_min_mz"]
|
||||
raw_max = entry["metadata"]["global_max_mz"]
|
||||
min_val = isa(raw_min, Dict) ? get(raw_min, "value", raw_min) : raw_min
|
||||
max_val = isa(raw_max, Dict) ? get(raw_max, "value", raw_max) : raw_max
|
||||
set_global_mz_range!(global_msi_data[], convert(Float64, min_val), convert(Float64, max_val))
|
||||
set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val))
|
||||
else
|
||||
precompute_analytics(global_msi_data[])
|
||||
precompute_analytics(msi_data)
|
||||
end
|
||||
end
|
||||
|
||||
@ -2636,7 +2599,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
end
|
||||
end
|
||||
|
||||
plotdata, plotlayout, xSpectraMz, ySpectraMz = meanSpectrumPlot(global_msi_data[], selected_folder_main, mask_path=mask_path_for_plot)
|
||||
plotdata, plotlayout, xSpectraMz, ySpectraMz = meanSpectrumPlot(msi_data, selected_folder_main, mask_path=mask_path_for_plot)
|
||||
plotdata_before = plotdata
|
||||
plotlayout_before = plotlayout
|
||||
last_plot_type = "mean"
|
||||
@ -2644,7 +2607,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
fTime = time()
|
||||
eTime = round(fTime - sTime, digits=3)
|
||||
msg = "Plot loaded in $(eTime) seconds"
|
||||
log_memory_usage("Mean Plot Generated", global_msi_data[])
|
||||
log_memory_usage("Mean Plot Generated", msi_data)
|
||||
catch e
|
||||
msg = "Could not generate mean spectrum plot: $e"
|
||||
warning_msg = true
|
||||
@ -2698,21 +2661,21 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
return
|
||||
end
|
||||
|
||||
if global_msi_data[] === nothing || full_route != target_path
|
||||
if global_msi_data[] !== nothing
|
||||
close(global_msi_data[])
|
||||
if msi_data === nothing || full_route != target_path
|
||||
if msi_data !== nothing
|
||||
close(msi_data)
|
||||
end
|
||||
msg = "Reloading $(basename(target_path)) for analysis..."
|
||||
full_route = target_path
|
||||
global_msi_data[] = OpenMSIData(target_path)
|
||||
msi_data = OpenMSIData(target_path)
|
||||
if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing
|
||||
raw_min = entry["metadata"]["global_min_mz"]
|
||||
raw_max = entry["metadata"]["global_max_mz"]
|
||||
min_val = isa(raw_min, Dict) ? get(raw_min, "value", raw_min) : raw_min
|
||||
max_val = isa(raw_max, Dict) ? get(raw_max, "value", raw_max) : raw_max
|
||||
set_global_mz_range!(global_msi_data[], convert(Float64, min_val), convert(Float64, max_val))
|
||||
set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val))
|
||||
else
|
||||
precompute_analytics(global_msi_data[])
|
||||
precompute_analytics(msi_data)
|
||||
end
|
||||
end
|
||||
local mask_path_for_plot::Union{String, Nothing} = nothing
|
||||
@ -2724,7 +2687,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
end
|
||||
end
|
||||
|
||||
plotdata, plotlayout, xSpectraMz, ySpectraMz = sumSpectrumPlot(global_msi_data[], selected_folder_main, mask_path=mask_path_for_plot)
|
||||
plotdata, plotlayout, xSpectraMz, ySpectraMz = sumSpectrumPlot(msi_data, selected_folder_main, mask_path=mask_path_for_plot)
|
||||
plotdata_before = plotdata
|
||||
plotlayout_before = plotlayout
|
||||
last_plot_type = "sum"
|
||||
@ -2732,7 +2695,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
fTime = time()
|
||||
eTime = round(fTime - sTime, digits=3)
|
||||
msg = "Total plot loaded in $(eTime) seconds"
|
||||
log_memory_usage("Sum Plot Generated", global_msi_data[])
|
||||
log_memory_usage("Sum Plot Generated", msi_data)
|
||||
catch e
|
||||
msg = "Could not generate total spectrum plot: $e"
|
||||
warning_msg = true
|
||||
@ -2789,21 +2752,21 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
return
|
||||
end
|
||||
|
||||
if global_msi_data[] === nothing || full_route != target_path
|
||||
if global_msi_data[] !== nothing
|
||||
close(global_msi_data[])
|
||||
if msi_data === nothing || full_route != target_path
|
||||
if msi_data !== nothing
|
||||
close(msi_data)
|
||||
end
|
||||
msg = "Reloading $(basename(target_path)) for analysis..."
|
||||
full_route = target_path
|
||||
global_msi_data[] = OpenMSIData(target_path)
|
||||
msi_data = OpenMSIData(target_path)
|
||||
if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing
|
||||
raw_min = entry["metadata"]["global_min_mz"]
|
||||
raw_max = entry["metadata"]["global_max_mz"]
|
||||
min_val = isa(raw_min, Dict) ? get(raw_min, "value", raw_min) : raw_min
|
||||
max_val = isa(raw_max, Dict) ? get(raw_max, "value", raw_max) : raw_max
|
||||
set_global_mz_range!(global_msi_data[], convert(Float64, min_val), convert(Float64, max_val))
|
||||
set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val))
|
||||
else
|
||||
precompute_analytics(global_msi_data[])
|
||||
precompute_analytics(msi_data)
|
||||
end
|
||||
end
|
||||
|
||||
@ -2818,7 +2781,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
|
||||
# Convert to positive coordinates for processing
|
||||
y_positive = yCoord < 0 ? abs(yCoord) : yCoord
|
||||
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, plotlayout, xSpectraMz, ySpectraMz, spectrum_id = xySpectrumPlot(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"
|
||||
@ -2859,7 +2822,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
fTime = time()
|
||||
eTime = round(fTime - sTime, digits=3)
|
||||
msg = "Plot loaded in $(eTime) seconds"
|
||||
log_memory_usage("XY Plot Generated", global_msi_data[])
|
||||
log_memory_usage("XY Plot Generated", msi_data)
|
||||
catch e
|
||||
msg = "Could not retrieve spectrum: $e"
|
||||
warning_msg = true
|
||||
@ -2904,21 +2867,21 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
return
|
||||
end
|
||||
|
||||
if global_msi_data[] === nothing || full_route != target_path
|
||||
if global_msi_data[] !== nothing
|
||||
close(global_msi_data[])
|
||||
if msi_data === nothing || full_route != target_path
|
||||
if msi_data !== nothing
|
||||
close(msi_data)
|
||||
end
|
||||
msg = "Reloading $(basename(target_path)) for analysis..."
|
||||
full_route = target_path
|
||||
global_msi_data[] = OpenMSIData(target_path)
|
||||
msi_data = OpenMSIData(target_path)
|
||||
if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing
|
||||
raw_min = entry["metadata"]["global_min_mz"]
|
||||
raw_max = entry["metadata"]["global_max_mz"]
|
||||
min_val = isa(raw_min, Dict) ? get(raw_min, "value", raw_min) : raw_min
|
||||
max_val = isa(raw_max, Dict) ? get(raw_max, "value", raw_max) : raw_max
|
||||
set_global_mz_range!(global_msi_data[], convert(Float64, min_val), convert(Float64, max_val))
|
||||
set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val))
|
||||
else
|
||||
precompute_analytics(global_msi_data[])
|
||||
precompute_analytics(msi_data)
|
||||
end
|
||||
end
|
||||
|
||||
@ -2932,7 +2895,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
end
|
||||
|
||||
# Call the new nSpectrumPlot function
|
||||
plotdata, plotlayout, xSpectraMz, ySpectraMz, spectrum_id = nSpectrumPlot(global_msi_data[], idSpectrum, selected_folder_main, mask_path=mask_path_for_plot)
|
||||
plotdata, plotlayout, xSpectraMz, ySpectraMz, spectrum_id = nSpectrumPlot(msi_data, idSpectrum, selected_folder_main, mask_path=mask_path_for_plot)
|
||||
plotdata_before = plotdata
|
||||
plotlayout_before = plotlayout
|
||||
last_plot_type = "single"
|
||||
@ -2942,7 +2905,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
fTime = time()
|
||||
eTime = round(fTime - sTime, digits=3)
|
||||
msg = "Plot loaded in $(eTime) seconds"
|
||||
log_memory_usage("nSpectrum Plot Generated", global_msi_data[])
|
||||
log_memory_usage("nSpectrum Plot Generated", msi_data)
|
||||
catch e
|
||||
msg = "Could not retrieve spectrum: $e"
|
||||
warning_msg = true
|
||||
@ -3223,7 +3186,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
|
||||
# This handler will now correctly load the first image from the newly selected folder.
|
||||
@onchange selected_folder_main begin
|
||||
# The global_msi_data[] object lifecycle is managed by the btnSearch handler.
|
||||
# The 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)
|
||||
@ -3455,7 +3418,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
fTime = time()
|
||||
eTime = round(fTime - sTime, digits=3)
|
||||
msg = "Plot loaded in $(eTime) seconds"
|
||||
log_memory_usage("Mean Plot Generated", global_msi_data[])
|
||||
log_memory_usage("Mean Plot Generated", msi_data)
|
||||
catch e
|
||||
msg = "Failed to load and process image: $e"
|
||||
warning_msg = true
|
||||
@ -3513,7 +3476,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
fTime = time()
|
||||
eTime = round(fTime - sTime, digits=3)
|
||||
msg = "Plot loaded in $(eTime) seconds"
|
||||
log_memory_usage("Mean Plot Generated", global_msi_data[])
|
||||
log_memory_usage("Mean Plot Generated", msi_data)
|
||||
catch e
|
||||
msg = "Failed to load and process image: $e"
|
||||
warning_msg = true
|
||||
@ -3614,7 +3577,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
# To include a visualization in the spectrum plot indicating where is the selected mass
|
||||
@onchange Nmass begin
|
||||
if !isempty(xSpectraMz)
|
||||
df = global_msi_data[].spectrum_stats_df
|
||||
df = 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)
|
||||
@ -3873,7 +3836,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
|
||||
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", global_msi_data[])
|
||||
log_memory_usage("App Ready", msi_data)
|
||||
end
|
||||
# is_processing = false
|
||||
GC.gc() # Trigger garbage collection
|
||||
|
||||
@ -546,12 +546,7 @@
|
||||
</q-list>
|
||||
|
||||
<!-- Pipeline Controls -->
|
||||
<div class="row q-mt-md justify-between items-center text-caption text-grey">
|
||||
<div>
|
||||
<q-icon name="info" class="q-mr-xs" /> Feature matrices are exported in highly optimized MatrixMarket Coordinate format (<b>.mtx</b>).
|
||||
<br/><i>Python:</i> <code>scipy.io.mmread('filtered_data.mtx')</code> | <i>R:</i> <code>Matrix::readMM('filtered_data.mtx')</code>
|
||||
</div>
|
||||
<div>
|
||||
<div class="row justify-end items-center q-mt-md">
|
||||
<q-btn class="q-ma-sm" icon="get_app" v-on:click="export_params_btn=true" label="Export Params" outline
|
||||
:disable="is_processing" />
|
||||
<q-btn class="q-ma-sm" icon="upload_file" v-on:click="import_params_btn=true" label="Import Params" outline
|
||||
@ -559,7 +554,6 @@
|
||||
<q-btn :loading="is_processing" class="q-ma-sm btn-style" icon="play_arrow"
|
||||
v-on:click="run_full_pipeline=true" padding="lg" label="Run Pipeline" :disable="is_processing" />
|
||||
</div>
|
||||
</div>
|
||||
</q-tab-panel>
|
||||
<q-tab-panel name="generator">
|
||||
<div class="text-h6">imzML & mzML Data Processor</div>
|
||||
|
||||
@ -1,35 +1,27 @@
|
||||
# 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.")
|
||||
println("This process can take 5-10 minutes.")
|
||||
|
||||
# Define the image path with OS-specific extension
|
||||
ext = Sys.iswindows() ? ".dll" : Sys.isapple() ? ".dylib" : ".so"
|
||||
sysimage_name = headless ? "sys_msi_headless" * ext : "sys_msi_gui" * ext
|
||||
sysimage_path = joinpath(@__DIR__, sysimage_name)
|
||||
sysimage_path = "MSI_sysimage" * ext
|
||||
|
||||
packages_to_include = Symbol[
|
||||
:MSI_src, # <--- Bake our core library
|
||||
:DataFrames,
|
||||
:SparseArrays,
|
||||
:Mmap,
|
||||
# List dependencies to include in the sysimage for maximum stability
|
||||
# We have removed graphical heavy-lifters (CairoMakie, PlotlyBase) to ensure
|
||||
# the build completes successfully on systems with standard RAM.
|
||||
packages_to_include = [
|
||||
:Genie,
|
||||
:GenieFramework,
|
||||
:PlotlyBase,
|
||||
:Libz,
|
||||
:DataFrames,
|
||||
:Serialization,
|
||||
:Printf,
|
||||
:JSON
|
||||
]
|
||||
|
||||
if !headless
|
||||
append!(packages_to_include, [
|
||||
:Genie,
|
||||
:GenieFramework,
|
||||
:PlotlyBase
|
||||
])
|
||||
end
|
||||
|
||||
create_sysimage(
|
||||
packages_to_include,
|
||||
sysimage_path = sysimage_path,
|
||||
@ -39,11 +31,5 @@ create_sysimage(
|
||||
)
|
||||
|
||||
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
|
||||
println("To start Julia with this image, use:")
|
||||
println("julia --threads auto --project=. --sysimage $(sysimage_path) start_MSI_GUI.jl")
|
||||
|
||||
@ -19,9 +19,9 @@ function create_warmup_datasets(dir, T_mz::Type, T_int::Type)
|
||||
mz_bytes = n_points * sizeof(T_mz)
|
||||
int_bytes = n_points * sizeof(T_int)
|
||||
|
||||
# Write some dummy binary data (m/z must be sorted for validation)
|
||||
# Write some dummy binary data
|
||||
open(ibd_path, "w") do f
|
||||
write(f, sort(rand(T_mz, n_points)))
|
||||
write(f, rand(T_mz, n_points))
|
||||
write(f, rand(T_int, n_points))
|
||||
end
|
||||
|
||||
@ -83,7 +83,7 @@ function create_warmup_datasets(dir, T_mz::Type, T_int::Type)
|
||||
|
||||
# 2. Create minimal mzML (Base64)
|
||||
mzml_path = joinpath(dir, "warmup_$(T_mz)_$(T_int).mzML")
|
||||
b64_mz = base64encode(sort(rand(T_mz, n_points)))
|
||||
b64_mz = base64encode(rand(T_mz, n_points))
|
||||
b64_int = base64encode(rand(T_int, n_points))
|
||||
|
||||
mzml_content = """<?xml version="1.0" encoding="utf-8"?>
|
||||
@ -133,21 +133,6 @@ try
|
||||
imzml_data = MSI_src.load_imzml_lazy(imzml_path; use_mmap=true)
|
||||
MSI_src.get_multiple_mz_slices(imzml_data, [10.0, 20.0], 0.1)
|
||||
|
||||
# Exercise Sprint 2 Streaming Pipeline
|
||||
config = MSI_src.PipelineConfig(
|
||||
steps=[
|
||||
MSI_src.StreamingStep(:baseline_correction, Dict(:method => :snip, :iterations => 5)),
|
||||
MSI_src.StreamingStep(:normalization, Dict(:method => :tic)),
|
||||
MSI_src.StreamingStep(:peak_picking, Dict(:method => :profile, :snr_threshold => 3.0))
|
||||
],
|
||||
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)
|
||||
@ -166,8 +151,8 @@ try
|
||||
p = PlotlyBase.Plot(PlotlyBase.scatter(x=1:2, y=[1,2]))
|
||||
end
|
||||
catch e
|
||||
if !(e isa InterruptException)
|
||||
@warn "Global Warmup failed: $e"
|
||||
end
|
||||
catch
|
||||
# Silent fail for interrupt
|
||||
end
|
||||
println("Precompilation workload completed.")
|
||||
|
||||
@ -20,7 +20,30 @@ 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
|
||||
|
||||
@ -1,37 +1,77 @@
|
||||
# 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
|
||||
|
||||
# --- Buffer Pooling ---
|
||||
"""
|
||||
posix_madvise(buffer::AbstractArray, advice::Integer)
|
||||
BufferPool
|
||||
|
||||
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.
|
||||
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.
|
||||
"""
|
||||
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
|
||||
mutable struct BufferPool
|
||||
lock::ReentrantLock
|
||||
buffers::Dict{Int, Vector{Vector{UInt8}}}
|
||||
|
||||
function BufferPool()
|
||||
new(ReentrantLock(), Dict{Int, Vector{Vector{UInt8}}}())
|
||||
end
|
||||
end
|
||||
|
||||
# --- Buffer Pooling ---
|
||||
"""
|
||||
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
|
||||
|
||||
|
||||
"""
|
||||
@ -71,16 +111,11 @@ If no suitable buffer is available, a new one is allocated.
|
||||
- A `Vector{UInt8}` buffer.
|
||||
"""
|
||||
function get_buffer!(pool::SimpleBufferPool, size::Int)::Vector{UInt8}
|
||||
buf = lock(pool.lock) do
|
||||
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)
|
||||
|
||||
876
src/MSIData.jl
876
src/MSIData.jl
File diff suppressed because it is too large
Load Diff
@ -20,18 +20,7 @@ export OpenMSIData,
|
||||
MSIData,
|
||||
_iterate_spectra_fast,
|
||||
validate_spectrum,
|
||||
get_mz_slice,
|
||||
REGISTRY_LOCK,
|
||||
SpectrumMetadata,
|
||||
SpectrumAsset,
|
||||
SpectrumMode,
|
||||
CENTROID,
|
||||
PROFILE,
|
||||
MzMLSource,
|
||||
ImzMLSource,
|
||||
SimpleBufferPool,
|
||||
get_buffer!,
|
||||
release_buffer!
|
||||
REGISTRY_LOCK
|
||||
|
||||
# Define shared registry lock
|
||||
const REGISTRY_LOCK = ReentrantLock()
|
||||
@ -71,21 +60,9 @@ 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")
|
||||
@ -95,8 +72,6 @@ include("Preprocessing.jl")
|
||||
include("ImageProcessing.jl")
|
||||
include("Precalculations.jl")
|
||||
include("PreprocessingPipeline.jl")
|
||||
include("StreamingKernels.jl")
|
||||
include("StreamingPipeline.jl")
|
||||
|
||||
using Setfield # For immutable struct updates
|
||||
|
||||
@ -120,17 +95,13 @@ 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(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)
|
||||
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)
|
||||
else
|
||||
error("Unsupported file type: $norm_filepath. Please provide a .mzML or .imzML file.")
|
||||
error("Unsupported file type: $filepath. Please provide a .mzML or .imzML file.")
|
||||
end
|
||||
|
||||
# Apply spectrum type map if provided
|
||||
|
||||
@ -1314,16 +1314,18 @@ 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
|
||||
count = end_idx - start_idx + 1
|
||||
if count < 3
|
||||
if (end_idx - start_idx + 1) < 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 = float(intensity[peak_idx])
|
||||
A_est = intensity[peak_idx]
|
||||
# Mean (μ): m/z at peak intensity
|
||||
mu_est = float(mz[peak_idx])
|
||||
mu_est = 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)
|
||||
@ -1337,27 +1339,19 @@ function _fit_gaussian_and_r2(mz::AbstractVector{<:Real}, intensity::AbstractVec
|
||||
return 0.0
|
||||
end
|
||||
|
||||
# Calculate SS_res, mean_y in a single pass to avoid allocations
|
||||
SS_res = 0.0
|
||||
sum_y = 0.0
|
||||
# Gaussian function
|
||||
gaussian(x, A, mu, sigma) = A * exp.(-(x .- mu).^2 ./ (2 * sigma^2))
|
||||
|
||||
@inbounds for i in start_idx:end_idx
|
||||
x_val = float(mz[i])
|
||||
y_val = float(intensity[i])
|
||||
# Generate estimated Gaussian curve
|
||||
y_est = gaussian(x_data, A_est, mu_est, sigma_est)
|
||||
|
||||
# Gaussian function estimate
|
||||
y_est = A_est * exp(-((x_val - mu_est)^2) / (2 * sigma_est^2))
|
||||
# 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 += (y_val - y_est)^2
|
||||
sum_y += y_val
|
||||
end
|
||||
|
||||
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
|
||||
SS_res = sum((y_data .- y_est).^2)
|
||||
SS_tot = sum((y_data .- mean(y_data)).^2)
|
||||
|
||||
if SS_tot == 0
|
||||
return 1.0 # Perfect fit if all y_data are the same
|
||||
|
||||
@ -12,7 +12,6 @@ generation.
|
||||
# =============================================================================
|
||||
|
||||
using Statistics # For mean, median
|
||||
using SparseArrays
|
||||
using StatsBase # For mad (Median Absolute Deviation)
|
||||
using SavitzkyGolay # For SavitzkyGolay filtering
|
||||
using Dates # For now()
|
||||
@ -41,7 +40,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::AbstractMatrix{Float64}
|
||||
matrix::Array{Float64,2}
|
||||
mz_bins::Vector{Tuple{Float64,Float64}}
|
||||
sample_ids::Vector{Int}
|
||||
end
|
||||
@ -488,23 +487,32 @@ 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 the baseline estimate array once
|
||||
# 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
|
||||
b1 = collect(float.(y))
|
||||
b2 = similar(b1)
|
||||
|
||||
current_b = b1
|
||||
next_b = b2
|
||||
|
||||
for k in 1:iterations
|
||||
prev_val = b1[1]
|
||||
b1[1] = min(b1[1], b1[2])
|
||||
# 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
|
||||
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)
|
||||
next_b[i] = min(current_b[i], 0.5 * (current_b[i-1] + current_b[i+1]))
|
||||
end
|
||||
|
||||
# Return the final baseline estimate
|
||||
return b1
|
||||
# Swap references for the next iteration (no data copy here)
|
||||
current_b, next_b = next_b, current_b
|
||||
end
|
||||
|
||||
# Return the final baseline estimate (which is in current_b after the last swap)
|
||||
return current_b
|
||||
end
|
||||
|
||||
"""
|
||||
@ -712,32 +720,18 @@ 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}}[]
|
||||
|
||||
# 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)
|
||||
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
|
||||
|
||||
candidate_peak_indices = Int[]
|
||||
sizehint!(candidate_peak_indices, div(n, 10)) # Pre-allocate memory capacity
|
||||
|
||||
@inbounds for i in 2:n-1
|
||||
for i in 2:n-1
|
||||
left = max(1, i - half_window)
|
||||
right = min(n, i + half_window)
|
||||
|
||||
# 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
|
||||
# Prominence check
|
||||
prominence = ys[i] - max(minimum(@view ys[left:i]), minimum(@view ys[i: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 &&
|
||||
if ys[i] >= maximum(@view ys[left:right]) &&
|
||||
(ys[i] > snr_threshold * noise_level) &&
|
||||
(prominence > min_peak_prominence * ys[i])
|
||||
push!(candidate_peak_indices, i)
|
||||
@ -774,14 +768,7 @@ function detect_peaks_profile_core(mz::AbstractVector{<:Real}, y::AbstractVector
|
||||
|
||||
left = max(1, p_idx - half_window)
|
||||
right = min(n, p_idx + half_window)
|
||||
|
||||
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)
|
||||
prominence = ys[p_idx] - max(minimum(@view ys[left:p_idx]), minimum(@view ys[p_idx:right]))
|
||||
|
||||
push!(detected_peaks, (mz=peak_mz, intensity=peak_int, fwhm=fwhm_ppm, shape_r2=shape_r2, snr=peak_snr, prominence=prominence))
|
||||
end
|
||||
|
||||
@ -1,79 +0,0 @@
|
||||
# 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
|
||||
@ -1,317 +0,0 @@
|
||||
# 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
|
||||
@ -1,402 +0,0 @@
|
||||
# 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
|
||||
594
src/imzML.jl
594
src/imzML.jl
@ -1,4 +1,5 @@
|
||||
using Images, Statistics, CairoMakie, DataFrames, Printf, ColorSchemes, StatsBase, Mmap
|
||||
# src/imzML.jl
|
||||
using Images, Statistics, CairoMakie, DataFrames, Printf, ColorSchemes, StatsBase
|
||||
|
||||
"""
|
||||
This file provides a library for parsing `.imzML` and `.ibd` files in pure Julia.
|
||||
@ -194,6 +195,86 @@ function axes_config_img(stream::IO)
|
||||
return param_groups
|
||||
end
|
||||
|
||||
"""
|
||||
get_spectrum_tag_offset(stream)
|
||||
|
||||
Calculates the character offset within a `<spectrum>` 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 (.+)")
|
||||
first = 1
|
||||
|
||||
while true
|
||||
value = match(r"[^=]+\"([^\"]+)\"", tag.captures[1][first:end])
|
||||
if value === nothing
|
||||
break
|
||||
end
|
||||
first += value.offsets[1] + length(value.captures[1])
|
||||
offset += length(value.captures[1])
|
||||
end
|
||||
return offset
|
||||
end
|
||||
|
||||
"""
|
||||
get_spectrum_attributes(stream, hIbd)
|
||||
|
||||
Reads metadata to determine the byte offsets and data types for reading spectra.
|
||||
|
||||
# Arguments
|
||||
|
||||
- `stream::IO`: The stream to parse.
|
||||
- `hIbd::IO`: The IBD file handle.
|
||||
|
||||
# Returns
|
||||
|
||||
- `mz_first::Int`: The index of the first array
|
||||
- `array_lengths::Vector{Int}`: The lengths of the arrays.
|
||||
- `offset_mz::Int`: The byte offset of the mz array.
|
||||
- `offset_intensity::Int`: The byte offset of the intensity array.
|
||||
- `data_type_mz::DataType`: The data type of the mz array.
|
||||
- `data_type_intensity::DataType`: The data type of the intensity array.
|
||||
"""
|
||||
function get_spectrum_attributes(stream::IO, hIbd::IO)
|
||||
# Look for position x
|
||||
readuntil(stream, "IMS:1000050")
|
||||
x_skip = 0
|
||||
|
||||
# Look for position y
|
||||
readuntil(stream, "IMS:1000051")
|
||||
y_skip = 0
|
||||
|
||||
# Determine order of mz/intensity arrays
|
||||
pos_before = position(stream)
|
||||
|
||||
# Look for first external offset (could be mz or intensity)
|
||||
readuntil(stream, "external offset")
|
||||
current_line = readline(stream)
|
||||
|
||||
# Check which array comes first by looking at the param group reference
|
||||
mz_first = occursin("mzArray", current_line) ? 3 : 4
|
||||
|
||||
# Find array length - look for the first external array length after coordinates
|
||||
seek(stream, pos_before)
|
||||
readuntil(stream, "IMS:1000103")
|
||||
array_len_skip = 0
|
||||
|
||||
# Find spectrum end
|
||||
readuntil(stream, "</spectrum>")
|
||||
spectrum_end_skip = 0
|
||||
|
||||
return [x_skip, y_skip, mz_first, array_len_skip, spectrum_end_skip]
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
read_spectrum_block(stream::IO)
|
||||
|
||||
@ -374,7 +455,7 @@ end
|
||||
# Arguments
|
||||
|
||||
- `stream::IO`: The stream to parse.
|
||||
- `hIbd::IO`: The IBD file handle.
|
||||
- `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.
|
||||
@ -389,7 +470,7 @@ end
|
||||
|
||||
- `spectra_metadata::Vector{SpectrumMetadata}`: The spectrum metadata.
|
||||
"""
|
||||
function parse_imzml_spectrum_block(stream::IO, hIbd::IO, param_groups::Dict{String, SpecDim},
|
||||
function parse_imzml_spectrum_block(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)
|
||||
@ -406,8 +487,8 @@ function parse_imzml_spectrum_block(stream::IO, hIbd::IO, param_groups::Dict{Str
|
||||
@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, 0.0, 0.0)
|
||||
int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, Int64(0), 0, :intensity, 0.0, 0.0)
|
||||
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)
|
||||
spectra_metadata[j] = SpectrumMetadata(Int32(0), Int32(0), "", :sample, global_mode, mz_asset, int_asset)
|
||||
end
|
||||
break
|
||||
@ -431,8 +512,8 @@ function parse_imzml_spectrum_block(stream::IO, hIbd::IO, param_groups::Dict{Str
|
||||
|
||||
if length(mz_data) != 1 || length(int_data) != 1
|
||||
println("DEBUG: Spectrum $k is empty or invalid - creating placeholder metadata")
|
||||
mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, Int64(0), 0, :mz, 0.0, 0.0)
|
||||
int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, Int64(0), 0, :intensity, 0.0, 0.0)
|
||||
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]
|
||||
@ -446,9 +527,9 @@ function parse_imzml_spectrum_block(stream::IO, hIbd::IO, param_groups::Dict{Str
|
||||
end
|
||||
|
||||
mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, mz_info.offset,
|
||||
mz_is_compressed ? mz_info.encoded_length : mz_info.array_length, :mz, 0.0, 0.0)
|
||||
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, 0.0, 0.0)
|
||||
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)
|
||||
@ -474,7 +555,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, use_mmap::Bool=true)
|
||||
function load_imzml_lazy(file_path::String; cache_size::Int=100)
|
||||
println("DEBUG: Checking for .imzML file at $file_path")
|
||||
if !isfile(file_path)
|
||||
throw(FileFormatError("Provided path is not a file: $(file_path)"))
|
||||
@ -488,118 +569,283 @@ function load_imzml_lazy(file_path::String; cache_size::Int=100, use_mmap::Bool=
|
||||
|
||||
println("DEBUG: Opening file streams for .imzML and .ibd")
|
||||
stream = open(file_path, "r")
|
||||
|
||||
# --- 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
|
||||
ts_hIbd = ThreadSafeFileHandle(ibd_path)
|
||||
|
||||
try
|
||||
@debug "Parsing imzML header..."
|
||||
# --- NEW: Parse all header information in a more efficient single pass ---
|
||||
println("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
|
||||
@debug "Image dimensions: $(width)x$(height), $num_spectra spectra."
|
||||
println("DEBUG: Image dimensions: $(width)x$(height), $num_spectra spectra.")
|
||||
|
||||
# ... (format extraction logic stays the same) ...
|
||||
# [Simplified for brevity in replacement chunk, but keeping the logic]
|
||||
# Extract default formats from the parsed param_groups
|
||||
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; 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
|
||||
if group.Axis == 1
|
||||
mz_group = group
|
||||
elseif group.Axis == 2
|
||||
int_group = group
|
||||
end
|
||||
end
|
||||
|
||||
# Check for compression status once
|
||||
any_comp = mz_is_compressed || int_is_compressed
|
||||
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
|
||||
end
|
||||
|
||||
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,
|
||||
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)
|
||||
|
||||
@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: Metadata parsing complete.")
|
||||
|
||||
# Build coordinate map ...
|
||||
# Build coordinate map for imzML files
|
||||
println("DEBUG: Building 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.")
|
||||
|
||||
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)
|
||||
# --- NEW: Update acquisition mode based on spectrum parsing ---
|
||||
acq_mode_symbol = if global_mode == CENTROID
|
||||
:centroid
|
||||
elseif global_mode == PROFILE
|
||||
:profile
|
||||
else
|
||||
:unknown
|
||||
end
|
||||
|
||||
# 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."
|
||||
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(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)
|
||||
|
||||
# Close the XML stream as it's no longer needed
|
||||
close(stream)
|
||||
|
||||
return msi_data
|
||||
|
||||
catch e
|
||||
close(stream)
|
||||
# Check if handles exist before closing
|
||||
if @isdefined(ibd_handles)
|
||||
for h in ibd_handles
|
||||
isopen(h) && close(h)
|
||||
end
|
||||
end
|
||||
close(ts_hIbd) # Ensure IBD handle is closed on error
|
||||
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("<spectrum ", line)
|
||||
spectrum_start_line = line
|
||||
break
|
||||
end
|
||||
end
|
||||
|
||||
# Parse lines within the spectrum block
|
||||
while !eof(stream)
|
||||
if !occursin("<spectrum ", line) # Avoid re-reading the first line
|
||||
line = readline(stream)
|
||||
end
|
||||
|
||||
if occursin("</spectrum>", 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("<binaryDataArray", line)
|
||||
bda_lines = [line]
|
||||
if !occursin("</binaryDataArray>", line)
|
||||
while !eof(stream)
|
||||
bda_line = readline(stream)
|
||||
push!(bda_lines, bda_line)
|
||||
if occursin("</binaryDataArray>", 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
|
||||
|
||||
|
||||
# =============================================================================
|
||||
#
|
||||
@ -623,7 +869,7 @@ This optimized version uses binary search for efficiency.
|
||||
# Returns
|
||||
- The intensity (`Float64`) of the peak if found, otherwise `0.0`.
|
||||
"""
|
||||
@inline function find_mass(mz_array::AbstractVector{<:Real}, intensity_array::AbstractVector{<:Real},
|
||||
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)
|
||||
@ -684,91 +930,59 @@ 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)
|
||||
|
||||
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}
|
||||
|
||||
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)
|
||||
|
||||
# 1. Find all candidate spectra first for efficient filtering
|
||||
candidate_indices = Set{Int}()
|
||||
indices_to_check = masked_indices === nothing ? (1:length(data.spectra_metadata)) : masked_indices
|
||||
|
||||
for i in indices_to_check
|
||||
# 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)
|
||||
# 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
|
||||
@inbounds for mass_int in min_mass_int:max_mass_int
|
||||
if mass_int in bf
|
||||
for mass_int in min_mass_int:max_mass_int
|
||||
if mass_int in bloom_filters[i]
|
||||
found = true
|
||||
break
|
||||
end
|
||||
end
|
||||
!found && continue
|
||||
|
||||
if !found
|
||||
continue # Definitely not in this spectrum
|
||||
end
|
||||
end
|
||||
|
||||
candidate_count += 1
|
||||
@inbounds candidate_indices[candidate_count] = i
|
||||
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
|
||||
end
|
||||
|
||||
# Use a view of the pre-allocated vector to avoid collect() allocations
|
||||
valid_candidates = view(candidate_indices, 1:candidate_count)
|
||||
println("Found $(length(candidate_indices)) candidate spectra (filtered from $(length(indices_to_check)) initial spectra)")
|
||||
|
||||
# 2. Iterate using the optimized, low-allocation iterator
|
||||
# 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]
|
||||
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
|
||||
@inbounds slice_matrix[meta.y, meta.x] = intensity
|
||||
Base.Threads.atomic_add!(results_count, 1)
|
||||
end
|
||||
slice_matrix[meta.y, meta.x] = intensity
|
||||
results_count += 1
|
||||
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
|
||||
@ -793,14 +1007,13 @@ 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 and allow binary search
|
||||
# Sort masses to improve cache locality during search
|
||||
sorted_masses = sort(masses)
|
||||
n_masses = length(sorted_masses)
|
||||
|
||||
# 1. Initialize a dictionary to hold the output slice matrices (using Float32 for 50% RAM savings)
|
||||
slice_dict = Dict{Real, Matrix{Float32}}()
|
||||
# 1. Initialize a dictionary to hold the output slice matrices
|
||||
slice_dict = Dict{Real, Matrix{Float64}}()
|
||||
for mass in sorted_masses
|
||||
slice_dict[mass] = zeros(Float32, height, width)
|
||||
slice_dict[mass] = zeros(Float64, height, width)
|
||||
end
|
||||
|
||||
local masked_indices::Union{Set{Int}, Nothing} = nothing
|
||||
@ -816,50 +1029,42 @@ function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, t
|
||||
precompute_analytics(data)
|
||||
end
|
||||
|
||||
println("Filtering candidate spectra for $n_masses m/z values...")
|
||||
println("Filtering candidate spectra for $(length(masses)) m/z values...")
|
||||
stats_df = get_spectrum_stats(data)
|
||||
bloom_filters = get_bloom_filters(data)
|
||||
|
||||
# Use a BitSet for faster index tracking
|
||||
candidate_indices = BitSet()
|
||||
candidate_indices = Set{Int}()
|
||||
indices_to_check = masked_indices === nothing ? (1:length(data.spectra_metadata)) : masked_indices
|
||||
|
||||
# 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
|
||||
# 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
|
||||
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 already a candidate, no need to check again
|
||||
if i in candidate_indices
|
||||
continue
|
||||
end
|
||||
# NEW: Bloom filter check with discretization
|
||||
if bloom_filters !== nothing && !is_empty(bloom_filters[i])
|
||||
found_any = false
|
||||
@inbounds for m_idx in m_start_idx:m_end_idx
|
||||
mass = sorted_masses[m_idx]
|
||||
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_any = true
|
||||
found = true
|
||||
break
|
||||
end
|
||||
end
|
||||
found_any && break
|
||||
|
||||
if !found
|
||||
continue # Definitely not in this spectrum
|
||||
end
|
||||
if found_any
|
||||
push!(candidate_indices, i)
|
||||
end
|
||||
else
|
||||
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
|
||||
end
|
||||
@ -868,40 +1073,29 @@ 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]
|
||||
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]
|
||||
# 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
|
||||
# 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)
|
||||
slice_dict[mass][meta.y, meta.x] = intensity
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
# 5. Clean up - replaces NaNs with 0.0 directly in Float32 matrices
|
||||
# 5. Clean up and return
|
||||
for mass in sorted_masses
|
||||
replace!(slice_dict[mass], NaN32 => 0.0f0)
|
||||
replace!(slice_dict[mass], NaN => 0.0)
|
||||
end
|
||||
|
||||
println("Finished generating $n_masses slices in a single pass.")
|
||||
println("Finished generating $(length(masses)) slices in a single pass.")
|
||||
return slice_dict
|
||||
end
|
||||
|
||||
@ -1539,7 +1733,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{Real, Real};
|
||||
bounds::Tuple{Float64, Float64};
|
||||
use_triq::Bool=false, triq_prob::Float64=0.98,
|
||||
mask_path::Union{String, Nothing}=nothing)
|
||||
# Use the provided bounds instead of recalculating
|
||||
|
||||
76
src/mzML.jl
76
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, 0.0, 0.0)
|
||||
return SpectrumAsset(data_format, compression_flag, binary_offset, encoded_length, axis)
|
||||
end
|
||||
|
||||
# This function is updated to return the generic SpectrumMetadata struct
|
||||
@ -394,48 +394,38 @@ 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")
|
||||
|
||||
# --- 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]
|
||||
ts_stream = ThreadSafeFileHandle(file_path, "r")
|
||||
|
||||
try
|
||||
# --- NEW: Parse instrument metadata from header ---
|
||||
println("DEBUG: Parsing instrument metadata from header...")
|
||||
instrument_meta = parse_instrument_metadata_mzml(primary_handle)
|
||||
instrument_meta = parse_instrument_metadata_mzml(ts_stream.handle)
|
||||
|
||||
seekstart(primary_handle) # Reset stream after header 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("------------------------------------")
|
||||
|
||||
seekstart(ts_stream.handle) # Reset stream after header parsing
|
||||
|
||||
println("DEBUG: Finding index offset...")
|
||||
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
|
||||
|
||||
index_offset = find_index_offset(ts_stream.handle)
|
||||
println("DEBUG: Seeking to index list at offset $index_offset.")
|
||||
seek(primary_handle, index_offset)
|
||||
seek(ts_stream.handle, index_offset)
|
||||
|
||||
println("DEBUG: Searching for '<index name=\"spectrum\">'.")
|
||||
if find_tag(primary_handle, r"<index\s+name=\"spectrum\"") === nothing
|
||||
if find_tag(ts_stream.handle, r"<index\s+name=\"spectrum\"") === nothing
|
||||
throw(FileFormatError("Could not find spectrum index."))
|
||||
end
|
||||
println("DEBUG: Found spectrum index tag.")
|
||||
|
||||
println("DEBUG: Parsing spectrum offsets...")
|
||||
spectrum_offsets = parse_offset_list(primary_handle)
|
||||
spectrum_offsets = parse_offset_list(ts_stream.handle)
|
||||
if isempty(spectrum_offsets)
|
||||
throw(FileFormatError("No spectrum offsets found."))
|
||||
end
|
||||
@ -443,25 +433,31 @@ function load_mzml_lazy(file_path::String; cache_size::Int=100)
|
||||
println("DEBUG: Found $num_spectra spectrum offsets.")
|
||||
|
||||
println("DEBUG: Parsing metadata for each spectrum...")
|
||||
# Pre-allocate the metadata vector for better performance
|
||||
spectra_metadata = Vector{SpectrumMetadata}(undef, num_spectra)
|
||||
|
||||
# Use @inbounds for faster indexing in the loop
|
||||
@inbounds for i in 1:num_spectra
|
||||
spectra_metadata[i] = parse_spectrum_metadata(primary_handle, spectrum_offsets[i])
|
||||
spectra_metadata[i] = parse_spectrum_metadata(ts_stream.handle, spectrum_offsets[i])
|
||||
|
||||
# Progress reporting for large files
|
||||
if i % 1000 == 0
|
||||
println("DEBUG: Processed $i/$num_spectra spectra")
|
||||
end
|
||||
end
|
||||
println("DEBUG: Metadata parsing complete for all $num_spectra spectra.")
|
||||
|
||||
# Inferred global formats from first spectrum
|
||||
# Assuming uniform data formats, take from the first spectrum
|
||||
first_meta = spectra_metadata[1]
|
||||
mz_format = first_meta.mz_asset.format
|
||||
intensity_format = first_meta.int_asset.format
|
||||
println("DEBUG: Inferred global m/z format: $mz_format")
|
||||
println("DEBUG: Inferred global intensity format: $intensity_format")
|
||||
|
||||
# Determine overall acquisition mode ...
|
||||
num_centroid = count(m -> m.mode == CENTROID, spectra_metadata)
|
||||
num_profile = count(m -> m.mode == PROFILE, spectra_metadata)
|
||||
# --- NEW: Determine overall acquisition mode ---
|
||||
modes = [meta.mode for meta in spectra_metadata]
|
||||
num_centroid = count(m -> m == CENTROID, modes)
|
||||
num_profile = count(m -> m == PROFILE, modes)
|
||||
|
||||
acq_mode_symbol = if num_centroid > 0 && num_profile == 0
|
||||
:centroid
|
||||
@ -472,28 +468,26 @@ 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,
|
||||
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
|
||||
instrument_meta.vendor_preprocessing_steps # Add this new field
|
||||
)
|
||||
|
||||
source = MzMLSource(mzml_handles, mz_format, intensity_format, mmap_data)
|
||||
source = MzMLSource(ts_stream, mz_format, intensity_format)
|
||||
println("DEBUG: Creating MSIData object.")
|
||||
return MSIData(source, spectra_metadata, final_instrument_meta, (0, 0), nothing, cache_size)
|
||||
|
||||
catch e
|
||||
# Close all handles in the pool if initialization fails
|
||||
for h in mzml_handles
|
||||
isopen(h) && close(h)
|
||||
end
|
||||
close(ts_stream) # Ensure stream is closed on error
|
||||
rethrow(e)
|
||||
end
|
||||
end
|
||||
|
||||
@ -19,52 +19,6 @@ 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()
|
||||
|
||||
|
||||
@ -49,16 +49,12 @@ 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/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_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_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 2GB"), # Mouse Stomach
|
||||
BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/salida_Tims/Stomach_DHB.imzML",804.3,0.1, "Mouse Stomach 4GB"), # Mouse Stomach
|
||||
|
||||
|
||||
BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML",804.3,0.1, "Mouse Stomach"), # Mouse Stomach
|
||||
]
|
||||
|
||||
const NUM_REPETITIONS = 50 # Number of times to generate the image for averaging
|
||||
|
||||
@ -1,215 +0,0 @@
|
||||
# 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
|
||||
@ -1,48 +0,0 @@
|
||||
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
|
||||
@ -16,9 +16,8 @@ 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 = ""
|
||||
|
||||
#=
|
||||
|
||||
@ -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"
|
||||
|
||||
|
||||
@ -31,10 +31,9 @@ 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 ---
|
||||
@ -59,17 +58,16 @@ 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/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML"
|
||||
const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML" # centroid
|
||||
# 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.2
|
||||
const MZ_TOLERANCE = 0.1
|
||||
|
||||
# Coordinates to plot a specific spectrum from imzML
|
||||
const COORDS_TO_PLOT = (50, 50) # Example coordinates (X, Y)
|
||||
|
||||
@ -1,58 +0,0 @@
|
||||
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
|
||||
@ -1,124 +0,0 @@
|
||||
#!/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()
|
||||
Loading…
x
Reference in New Issue
Block a user