diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md new file mode 100644 index 0000000..385ecf9 --- /dev/null +++ b/.github/pull_request_template.md @@ -0,0 +1,4 @@ +## Pull Request (PR) description + +## Issues resolved +n/a \ No newline at end of file diff --git a/.github/workflows/main_publish-release.yaml b/.github/workflows/main_publish-release.yaml new file mode 100644 index 0000000..6af4750 --- /dev/null +++ b/.github/workflows/main_publish-release.yaml @@ -0,0 +1,42 @@ +name: Main Branch Create New Release Workflow + +on: + push: + branches: [main] + +permissions: + contents: write + +jobs: + release: + name: Build and publish release + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v7 + + - name: Install uv + uses: astral-sh/setup-uv@ae62891fec2bb8e7d6c99fc78c9fec3a63790f8d + + - name: Read version + id: version + run: | + VERSION=$(grep -m1 '^version *= *' pyproject.toml | sed -E 's/version *= *"([^"]+)"/\1/') + if [ -z "$VERSION" ]; then + echo "Could not read version from pyproject.toml" >&2 + exit 1 + fi + echo "version=$VERSION" >> "$GITHUB_OUTPUT" + + - name: Build package + run: uv build + + - name: Publish release + uses: softprops/action-gh-release@v3 + with: + tag_name: v${{ steps.version.outputs.version }} + name: v${{ steps.version.outputs.version }} + generate_release_notes: true + draft: false + files: | + dist/*.whl + dist/*.tar.gz diff --git a/.github/workflows/main_version-update-check.yaml b/.github/workflows/main_version-update-check.yaml new file mode 100644 index 0000000..92edc24 --- /dev/null +++ b/.github/workflows/main_version-update-check.yaml @@ -0,0 +1,66 @@ +name: Main Branch Version Check Workflow + +on: + pull_request: + branches: [main] + +permissions: + contents: read + +jobs: + version-check: + name: Check version incremented + runs-on: ubuntu-latest + steps: + - name: Checkout PR branch + uses: actions/checkout@v7 + + - name: Read PR version + id: pr_version + run: | + VERSION=$(grep -m1 '^version *= *' pyproject.toml | sed -E 's/version *= *"([^"]+)"/\1/') + if [ -z "$VERSION" ]; then + echo "Could not read version from pyproject.toml" >&2 + exit 1 + fi + echo "version=$VERSION" >> "$GITHUB_OUTPUT" + + - name: Read main version + id: main_version + run: | + git fetch origin main --depth=1 + VERSION=$(git show origin/main:pyproject.toml | grep -m1 '^version *= *' | sed -E 's/version *= *"([^"]+)"/\1/') + if [ -z "$VERSION" ]; then + echo "Could not read version from main's pyproject.toml" >&2 + exit 1 + fi + echo "version=$VERSION" >> "$GITHUB_OUTPUT" + + - name: Install uv + uses: astral-sh/setup-uv@ae62891fec2bb8e7d6c99fc78c9fec3a63790f8d + + - name: Compare versions + run: | + PR_VERSION="${{ steps.pr_version.outputs.version }}" + MAIN_VERSION="${{ steps.main_version.outputs.version }}" + echo "PR version: $PR_VERSION" + echo "main version: $MAIN_VERSION" + uv run --with packaging python3 - "$PR_VERSION" "$MAIN_VERSION" <<'EOF' + import sys + from packaging.version import Version, InvalidVersion + + pr_raw, main_raw = sys.argv[1], sys.argv[2] + + try: + pr = Version(pr_raw) + main = Version(main_raw) + except InvalidVersion as e: + print(f"Could not parse version as PEP 440: {e}", file=sys.stderr) + sys.exit(1) + + if pr <= main: + print(f'Version in pyproject.toml ({pr_raw}) is not greater than main ({main_raw}). Versions must strictly increase.', file=sys.stderr) + sys.exit(1) + + print(f'OK: {pr_raw} > {main_raw}') + EOF diff --git a/docs/analyse.md b/docs/analyse.md index b31cbe2..46529b5 100644 --- a/docs/analyse.md +++ b/docs/analyse.md @@ -28,6 +28,7 @@ Arguments: Options: -o, --out-dir PATH Path to output directory. Results will be written under '.../evaluator/analyse/'. [default: .] --fill-threshold FLOAT [0-1] Override configuration fill threshold parameter for this run. + --membrane-thickness-nm FLOAT [x>0] Override membrane shell thickness (nm) used for diameter derivation, instead of estimating from file. -h, --help Show this message and exit. Batch Options: @@ -59,6 +60,9 @@ fill_ratio = (filled_volume - original_volume) / filled_volume A component is classified as enclosed (`is_enclosed = True`) if `fill_ratio > fill_threshold`. The default of `0.05` is deliberately permissive, to accommodate incomplete segmentations where a small fraction of the membrane may not have been captured. Decreasing this threshold includes more components as enclosed; increasing it requires a larger enclosed cavity relative to the membrane. +#### `--membrane-thickness-nm` +This option controls the assumed membrane thickness value used during derivation of a vesicle's outer diameter from the meaured membrane volume. The default value (`"auto"`) causes EValuator to estimate the thickness for each segmentation file as twice the median distance-transform value along the ridge of the membrane mask. Using any other value will override this estimation for every file. + #### `--qc-max-sphere-rmse-rel`, `--qc-max-aspect-ratio`, `--qc-min-solidity`, `--qc-min-arc-coverage` These four thresholds underpin the vesicle-vs-debris QC flag included in results (`is_vesicle_like` / `qc_flags`, see [Output CSV columns](#output-csv-columns)). @@ -87,7 +91,8 @@ Column | Description | Units[^units] --|--|-- `tomogram` | Filename of the input MRC file from which the EV was identified | `label` | Unique integer identifier for the EV within its file | -`equiv_diameter_nm`[^equivdiameter] | Diameter of a sphere with the same volume as the membrane component | nm (rounded to 2 d.p.) +`equiv_diameter_nm`[^equivdiameter] | Outer diameter of the vesicle, derived from the membrane shell volume and thickness | nm (rounded to 2 d.p.) +`xy_diameter_nm` | Outer diameter of the widest xy slice (mean of y and x extents), independent of z | nm (rounded to 2 d.p.) `major_axis_diameter`[^diameters] | Length of the longest axis of the best-fit ellipsoid | nm (rounded to 2 d.p.) `minor_axis_diameter`[^diameters] | Length of the shortest axis of the best-fit ellipsoid | nm (rounded to 2 d.p.) `aspect_ratio`[^aspectratio] | Ratio of `major_axis_diameter` to `minor_axis_diameter` | unitless (rounded to 2 d.p.) @@ -101,8 +106,11 @@ Column | Description | Units[^units] `solidity` | Voxel count ÷ convex-hull volume (QuickHull over the component's voxels; ~1 for a solid convex blob, low for sprawling/concave debris) | unitless (0 < s ≲ 1; rounded to 4 d.p.) `arc_coverage` | Fraction of the fitted sphere's surface (18×36 grid) occupied by the component | unitless (0 ≤ c ≤ 1; rounded to 4 d.p.; `NaN` if < 4 voxels) `bbox_extent` | Bounding-box fill ratio (`skimage` `regionprops.extent`) | unitless (0 < e ≤ 1; rounded to 4 d.p.) +`touches_z_boundary` | Whether the component's bounding box reaches the first or last z slice (and component therefore may be clipped by reconstruction box) | boolean (True/False) `is_vesicle_like` | Composite QC flag — `True` unless a vesicle-vs-debris check fails (see [`--qc-*` options](#--qc-max-sphere-rmse-rel---qc-max-aspect-ratio---qc-min-solidity---qc-min-arc-coverage)) | boolean (True/False) `qc_flags` | Comma-joined names of failed QC checks (`sphere_rmse`, `aspect_ratio`, `solidity`, `open_shell`, `too_small`); empty when clean | string +`membrane_thickness_nm` | Membrane shell thickness used to derive `equiv_diameter_nm` for this file | nm (rounded to 2 d.p.) +`membrane_thickness_source` | How membrane thickness was set: `estimated` (per-file estimate), `override` (config number or CLI) or `fallback` (estimate failed, `2.6` nm used) | text `voxel_size_nm` | Voxel size in nanometres as read from MRC file header, or None if not present | nm (rounded to 4 d.p.) `measurement_units` | Units used for measurements | nm if voxel size was available, otherwise vox @@ -118,7 +126,7 @@ Pipeline run summary - Segmentation files with EVs: 9 (90.0%) - EVs processed: 87 - Number of enclosed EVs: 71 (81.6%) -- Equivalent diameters: 112.4 ± 48.3 nm (mean ± SD) +- Outer diameters: 112.4 ± 48.3 nm (mean ± SD) - Vesicle-like components: 63 / 87 (72.4%) Results saved to: .../evaluator/analyse/evaluator-analyse_results.csv ``` @@ -132,14 +140,14 @@ Results saved to: .../evaluator/analyse/evaluator-analyse_results.csv [^units]: Units given here assume the voxel size in nanometres was read from the MRC file header. If this is not the case, units will not be scaled to physical units and will be in voxels/voxels³. -[^equivdiameter]: `equiv_diameter_nm` is not the measured diameter of an EV. Given the volume of all membrane voxels (Vm), the equivalent diameter is calculated as `(6Vm/π)**(1/3)`, which assumes the shape of the component is a perfect sphere. This is likely not a valid assumption for biological EVs, and `equiv_diameter_nm` is therefore only used as a rough proxy during size filtering and for subsequent axis scaling. +[^equivdiameter]: `equiv_diameter_nm` is not the measured diameter of an EV but a derived outer diameter estimate. Given the volume of all membrane voxels (Vm), the equivalent diameter is derived from solving `4π/3·(R^3 - (R-t)**3) = Vm` for the outer radius `R` (which provides equivalent diameter via `equiv_diameter_nm = 2R`) and membrane thickness `t` (estimated or provided via `--membrane-thickness-nm`). This estimate assumes a fixed membrane thickness per-file and roughly spherical shell, so is unlikely to be a valid model for biological EVs, and is therefore only used as a rough proxy during size filtering and for subsequent axis scaling. -[^diameters]: `major/minor_axis_diameter` are the more accurate measurements of EV size, calculated from the best-fit ellipsoid that matches the EV morphology. Both measurements are derived from the eigenvalues of the component's inertia tensor, which are scaled using `equiv_diameter_nm`. +[^diameters]: `major/minor_axis_diameter` approximate EV diameters using the best-fit ellipsoid implied by the components inertia tensor eigenvalues after scaling using `equiv_diameter_nm`. They inherit the equivalent diameter's assumptions, but the ratio between diameters reflects the true shape. [^aspectratio]: A perfect sphere would have an `aspect_ratio` of 1.0. An EV with `aspect_ratio` greater than 1.0 can be described as a prolate ellipsoid (i.e. elongated), whereas an EV with `aspect_ratio` less than 1.0 can be described as an oblate ellipsoid (i.e. flattened). [^eccentricity]: Eccentricity is calculated as `sqrt(1 - (c/a)^2)` where `a` and `c` are the semi-major and semi-minor axes respectively. A perfect sphere would have `eccentricity` approaching 0, whereas an infinitely elongated ellipsoid would have `eccentricity` approaching 1. -[^lumenvol]: This volume is calculated by filling the membrane mask and subtracting all membrane voxels. Non-enclosed components will have `lumen_volume` of 0. +[^lumenvol]: This volume is calculated by filling the gap-bridged membrane mask (used for is_enclosed check) and subtracting all membrane voxels. Non-enclosed components will have `lumen_volume` of 0. [^surfacearea]: Surface area is computed using the marching cubes algorithm as implemented by `skimage.measure.marching_cubes`. If the marching cubes algorithm fails (e.g. for very small or degenerate components), `NaN` will be returned. \ No newline at end of file diff --git a/docs/config.md b/docs/config.md index 1b312e8..b1e3d2d 100644 --- a/docs/config.md +++ b/docs/config.md @@ -100,11 +100,13 @@ format = "json" # Label command default configuration parameters [label] min_arc_coverage = 0.30 -merge_centre_tol_factor = 1.5 +merge_centre_tol_factor = 0.573 merge_radius_tol_pct = 0.30 +merge_concentric_radius_pct = 0.10 +merge_concentric_centre_frac = 0.25 maximum_diameter_nm = 500.0 minimum_diameter_nm = 20.0 -membrane_thickness_nm = 7 +membrane_thickness_nm = 2.6 max_workers = 0 # Model command default configuration parameters @@ -117,6 +119,7 @@ max_workers = 0 # Analyse command default configuration parameters [analyse] fill_threshold = 0.05 +membrane_thickness_nm = "auto" max_workers = 0 qc_max_sphere_rmse_rel = 0.25 qc_max_aspect_ratio = 2.5 diff --git a/docs/label.md b/docs/label.md index 46448c2..6a911ff 100644 --- a/docs/label.md +++ b/docs/label.md @@ -53,7 +53,7 @@ After merging, each component's coverage of the surface of its best-fit sphere i These options filter components by their equivalent sphere diameter, using the configured membrane thickness to convert diameter to an expected voxel count range. The default range (`20`-`500` nm, set in `config.toml`) is appropriate for typical EV preparations, which include small exosomes (~30-150 nm) through to larger MVB-derived vesicles (~200-500 nm). If no voxel size is present in the input MRC header, this filter is skipped and all components pass. #### Component merge configuration settings -The merge step (centroid-proximity and radius-agreement tolerances) is controlled separately via `merge_centre_tol_factor` and `merge_radius_tol_pct` in `config.toml` under `[label]`; there is no CLI override for these. +The merge step (centroid-proximity and radius-agreement tolerances) is controlled separately via `merge_centre_tol_factor`, `merge_radius_tol_pct`, `merge_concentric_radius_pct` and `merge_concentric_centre_frac` in `config.toml` under `[label]`; there is no CLI override for these. Pairs whose centres are within `merge_concentric_centre_frac` (default a quarter) of the smaller radius but have radii that differ by more than `merge_concentric_radius_pct` are treated as concentric membranes (MVBs, multilamellar vesicles) and are not merged. Concentric pairs with radii differing by less than this fraction are indistinguishable from split fragments, and thus merged to a single component. ## Output diff --git a/docs/plot.md b/docs/plot.md index 5f09571..b03b8ba 100644 --- a/docs/plot.md +++ b/docs/plot.md @@ -3,13 +3,11 @@ ## Overview The `plot` command generates plots and summary tables from [`analyse`](analyse.md) output and/or [`model`](model.md) outputs. Each run produces one or more independent sections (depending on provided results), each writing SVG plots and `.xlsx` summary statistics under `evaluator/plot/
`. -`plot` requires R to be installed, and either available on `PATH` or passed using the `--rscript` option. - ## Pipeline description For a given run, the pipeline: 1. Resolves `--analyse`/`--model` into one or more runs (each option accepts either a direct result file or a sample sheet). 2. Determines which sections are runnable for the given inputs and which were requested on command line. -3. Dispatches each section to its own R script via `Rscript`, writing outputs to `evaluator/plot/
/`. A section directory that already exists is skipped unless `--overwrite` is given. A section that fails (non-zero `Rscript` exit) is logged as a warning, but should not stop the remaining sections from running. +3. Runs each section, writing outputs to `evaluator/plot/
/`. A section directory that already exists is skipped unless `--overwrite` is given. A section that fails is logged as a warning, but should not stop the remaining sections from running. 4. Writes `evaluator/plot/index.md`, linking every section that completed successfully and `evaluator/plot/params.toml`, the resolved `[plot]` parameters used for the run. ## Usage @@ -17,14 +15,15 @@ For a given run, the pipeline: Usage: evaluator plot [OPTIONS] Options: - --analyse PATH Path to an analyse results CSV, or a sample sheet referencing multiple. + --analyse PATH Path to an analyse results CSV, or a sample sheet referencing multiple. --model PATH Path to a model results file (JSON or CSV), or a sample sheet referencing multiple. -o, --out-dir PATH Path to output directory, outputs written under '.../evaluator/plot/'. [default: .] --section TEXT Section(s) to run: distributions, qc, scatter, concordance, compare (repeatable). --all Run every applicable section. --overwrite Overwrite existing section outputs instead of skipping them. - --rscript PATH Path to the Rscript binary (default: resolved from PATH). - -h, --help Show this message and exit. + --include-flagged Include components flagged as non-vesicle-like instead of excluding them. + --agreement TEXT,TEXT Bland-Altman plot for a pair of measurements (repeatable). + -h, --help Show this message and exit. ``` At least one of `--analyse`/`--model` must be provided. Global verbosity options (`-v` / `-vv`) are set on the `evaluator` command itself, and should be included before the `plot` subcommand: @@ -37,27 +36,41 @@ evaluator -vv plot --analyse evaluator/results/analyse/evaluator-analyse_results ### `--section` and `--all` By default (no `--section`, no `--all`), the sections listed under `[plot].default_sections` in the resolved configuration file are run. `--section` may be repeated to run a specific subset, or `--all` can be used to run every runnable section. +### QC filtering +Any vesicles where `analyse`'s `is_vesicle_like` column is `False` are excluded before computing any distribution, statistic or plot. The number of dropped rows and a breakdown of the reasons for exclusion are logged automatically. To include these vesicles in outputs, use the `--include-flagged` option. + +### Caveats +- `distributions` reports each summary statistic for three populations (`all`, `vesicle_like`, `vesicle_like_no_z_clip`), so the effect of the QC filter and of z clipping is visible. Conclusions which hold in all three should be independent of QC filters/z clipping. +- `scatter` and `correlations` figures and `stats`/`pearson`/`spearman` sheets pool all vesicles, treating each as an independent observation so should be considered descriptive only. Each section also writes a `*_tomogram` sheet computed on per-tomogram medians, which is a cluster-aware estimate. +- `compare` aggregates to one median per tomogram before plotting, summarising and testing (see below). Tomograms within a replicate are still non-independent, so `compare` p-values are optimistic if *nreplicate* < *ntomogram*. + ## Sections Section | Requires | Outputs --|--|-- -`distributions` | `--analyse` | A density + rug plot per morphology feature (`_density.svg`), and `summary_stats.xlsx` (*n*, mean, median, sd, p5, p95 per feature) -`qc` | `--analyse` | `count_per_tomogram.svg` (EV count per tomogram), `closure_pass_rate.svg` (fraction `is_enclosed` per tomogram), `closure_fill_ratio_distribution.svg` (with the configured flag threshold marked), and `qc_summary.xlsx` (per-tomogram EV count, pass rate, median fill ratio, count flagged below threshold) -`scatter` | `--analyse` | A scatter plot with fitted trend line per feature pair (`equiv_diameter_nm` vs `lumen_volume`, `surface_area` vs `lumen_volume`, `aspect_ratio` vs `equiv_diameter_nm`), and `scatter_stats.xlsx` (Pearson/Spearman correlation, slope, slope SE per pair) +`distributions` | `--analyse` | A histogram per morphology feature (`_histogram.svg`), and `summary_stats.xlsx` (*n*, mean, median, sd, p5, p95 per feature, for each population) +`qc` | `--analyse` | `count_per_tomogram.svg`, `closure_pass_rate.svg`, `closure_fill_ratio_distribution.svg` (configured flag threshold marked), and `qc_summary.xlsx` +`scatter` | `--analyse` | A scatter plot with fitted trend line and R^2 annotation per feature pair (`_vs_.svg`), and `scatter_stats.xlsx` (Pearson/Spearman correlation, slope, slope SE, R^2, *n* per pair; `stats` is pooled over vesicles, `stats_tomogram` over per-tomogram medians) +`correlations` | `--analyse` | `correlation_spearman.svg` and `correlation_pearson.svg` heatmaps across the morphology features, and `correlations.xlsx` `concordance` | `--analyse` **and** `--model` | See [Concordance](#concordance) below. -`compare` | Sample sheet with more than one distinct `group` (multi-sample mode) | Per-feature violin+box plots by group (`violin_.svg`), a PCA plot across standardised features (`pca.svg`, only if ≥3 complete-case rows and ≥2 features are available), and `group_summary.xlsx` (per-group n, mean, median per feature) +`reliability` | `--analyse` **and** `--model` | `reliability.svg` (fit RMSE against closure fill ratio, coloured by `reliability.is_reliable`) and `reliability_summary.xlsx` +`agreement` | `--analyse` **and** `--model`; runs when `--agreement` is given | See [Agreement](#agreement) below. +`compare` | Sample sheet with more than one distinct `group` | Per-feature violin+box plots by group of per-tomogram medians (`violin_.svg`), `pca.svg` (only if at least 3 complete-case tomograms and 2 features), `group_summary.xlsx` (`n_tomograms`, `n_samples`, `n_evs`, and `n_replicates` when the sheet has a `replicate` column, plus mean and median of per-tomogram medians) and `group_tests.xlsx` (Mann-Whitney U for two groups, Kruskal-Wallis for more, on per-tomogram medians) ### Concordance `concordance` provides a comparison of `model`'s least-squares-fitted geometry and the raw voxel measurements that `analyse` produces. -Rows are joined on `analyse`'s `(tomogram, label)` against `model`'s `(source_file, -label_id)`, matching `source_file`'s basename against `tomogram`, therefore both must have been run against the same labelled MRC(s) for rows to match. Only matched rows are plotted; `concordance_stats.xlsx` reports `n_matched`, `n_analyse_only`, and `n_model_only` so unmatched counts are visible. +Rows are joined on `analyse`'s `(tomogram, label)` against `model`'s `(source_file, label_id)`, matching `source_file`'s basename against `tomogram`, therefore both must have been run against the same labelled MRC(s) for rows to match. Only matched rows are plotted; `concordance_stats.xlsx` reports `n_matched`, `n_analyse_only`, and `n_model_only` so unmatched counts are visible. -Produces, per matching vesicle set: -- `model_vs_equiv_diameter_nm.svg` and `model_vs_major_axis_diameter.svg` — `model`'s fitted diameter (`2 × radius`) plotted against each of `analyse`'s two independent diameter estimates, with a `y = x` reference line. -- `reliability_vs_closure.svg` — whether `model`'s `reliability.is_reliable` gate tracks `analyse`'s independently-measured `closure_fill_ratio`. -- `fit_rmse_vs_enclosed.svg` — whether `model`'s fit RMSE tracks `analyse`'s `is_enclosed`. -- `concordance_stats.xlsx` — match counts (above) plus Pearson correlation between `model`'s fitted diameter and each of `analyse`'s two diameter estimates, and the reliable-fit pass rate among matched rows. +Produces `model_vs_equiv_diameter_nm.svg` and `model_vs_major_axis_diameter.svg` (`model`'s fitted diameter against each `analyse` diameter, with a `y = x` line, trend line and bias/RMSE annotation) and `concordance_stats.xlsx` (bias, RMSE and Pearson correlation per diameter, plus `n_matched`, `n_analyse_only`, `n_model_only`). With the default QC filter, `model` rows for components `analyse` flagged are counted in `n_model_only`. + +### Agreement +`--agreement A,B` creates a Bland-Altman plot (difference `A - B` against the mean, with bias and 95% limits of agreement) for each measurement pair, written to `agreement/_vs_.svg`, with `agreement_stats.xlsx`. The option can be repeated to produce several Bland-Altman plots: +```sh +evaluator plot --analyse analyse.csv --model model.json --agreement equiv_diameter_nm,model_diameter_nm --agreement major_axis_diameter,model_diameter_nm +``` + +`A` and `B` are column names from the joined analyse and model results. Pairs naming a missing column are skipped with a warning. `--section agreement` or `--all` without `--agreement` plots `equiv_diameter_nm,model_diameter_nm`. ## Sample sheets (multi-sample mode) `--analyse` and/or `--model` may point to a sample sheet instead of a single result file, which is a `.tsv` or `.txt` file with a header row containing `sample_id`. @@ -67,32 +80,26 @@ Column | Required | Description `sample_id` | Yes | Unique identifier for the sample/run; joins rows across the `--analyse` and `--model` sheets when both are sheets `path` | No | Path to that sample's result file. If omitted and the *other* option's value is a single file (not a sheet), that single file is treated as applying to every sample on this side `group` | No | Group label, used by the `compare` section. `compare` only becomes runnable when more than one distinct `group` value is present across the sheet. -`replicate` | No | Optional integer replicate number; parsed but not currently used by any section +`replicate` | No | Optional integer replicate number; `compare` reports the number of replicates per group If both `--analyse` and `--model` are sheets, rows are matched by `sample_id`; a `sample_id` present in only one sheet gets `None` for the missing side. If only one side is a sheet, every sample on that side is paired with the other side's single (non-sheet) file, if one was given. -## R prerequisites -`plot` requires `Rscript` to render plots and write spreadsheets, as well as requires the following R packages to be installed: `tidyverse`, `jsonlite`, `openxlsx`, -`viridisLite`. - -By default, `plot` resolves `Rscript` from `PATH`. If R isn't on `PATH` (e.g. an `renv`/`conda` R install), pass its location explicitly with `--rscript`: -```sh -evaluator plot --analyse evaluator/results/analyse/evaluator-analyse_results.csv --all --rscript /path/to/Rscript -``` - ## Output Results are written under the output directory (default: current working directory) to `evaluator/plot/`: -- `
/`: one directory per completed section, containing that section's SVG plots and `.xlsx` summary table(s) +- `
/`: one directory per completed section, containing that section's SVG plots and `.xlsx` summary table(s); in multi-sample mode, per-sample sections are split into `
//` - `index.md`: links to every section that completed successfully in this run - `params.toml`: the resolved `[plot]` parameters used for this run ### Configuration `plot` reads the following default values from the configuration file's `[plot]` section: -- `default_sections` (default `["distributions", "qc", "scatter"]`): sections run when - neither `--section` nor `--all` is given. -- `fill_ratio_flag_threshold` (default `0.05`): `closure_fill_ratio` values below this are marked in the `qc` section's fill-ratio histogram and counted in `qc_summary.xlsx`. +- `default_sections` (default `["distributions", "qc", "scatter", "correlations"]`): sections run when neither `--section` nor `--all` is given +- `fill_ratio_flag_threshold` (default `0.05`): `closure_fill_ratio` values below this are marked in the `qc` fill-ratio histogram and counted in `qc_summary.xlsx` +- `theme` (default `"Okabe-Ito"`): figure colours; one of the viewer's named themes +- `bin_width` (default `25.0`): histogram bin width in feature units, anchored at 0; `0` for automatic bins +- `[plot.bin_width_overrides]`: per-feature bin widths replacing `bin_width` +- `scatter_pairs`: the `[x, y]` feature pairs drawn by `scatter`
diff --git a/docs/tools.md b/docs/tools.md index 34c28d0..cca3b47 100644 --- a/docs/tools.md +++ b/docs/tools.md @@ -36,7 +36,7 @@ Benchmarks missing-wedge correction/fitting accuracy across a range of synthetic Defaults are defined in `evaluator.commands.tools.benchmark_missing_wedge.defaults`. ### `evaluator tools benchmark geometric-proxy` -Benchmarks the `estimateCentroidRadius` geometric proxy (naive vs isotropy-aware bounding box) across partial spherical caps and equatorial bands. +Benchmarks the `estimateCentroidRadius` geometric proxy (naive vs isotropy-aware bounding box) and an independent `diagonal/2 baseline` method across partial spherical caps and equatorial bands. | Option | Description | |---|---| diff --git a/pyproject.toml b/pyproject.toml index d6b87c4..1179e4b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "EValuator" -version = "1.0.0" +version = "2.0.0" description = "A command line tool for automated morphological analysis and visualisation of extracellular vesicles (EVs) from cryo-electron tomography (cryo-ET) data." authors = [ { name = "Samuel Holland", email = "hol_sam@icloud.com" } @@ -36,6 +36,11 @@ viewer = [ "plotly>=6.0", "streamlit>=1.57", ] +plot = [ + "kaleido>=1.0", + "openpyxl>=3.1.5", + "plotly>=6.0", +] [build-system] requires = ["uv_build<0.12"] @@ -49,7 +54,6 @@ include = [ "LICENSE", # Include the LICENSE file "README.md", # Include README or other files "src/evaluator/config.toml", - "src/evaluator/commands/plot/r/**", "src/evaluator/commands/viewer/app/**", "docs/*", ] diff --git a/src/evaluator/commands/analyse/analyse.py b/src/evaluator/commands/analyse/analyse.py index 32bc55a..b2ec41d 100644 --- a/src/evaluator/commands/analyse/analyse.py +++ b/src/evaluator/commands/analyse/analyse.py @@ -70,6 +70,7 @@ def analyse( worker = partial( processSegmentation, fill_threshold=params.fill_threshold, + membrane_thickness_nm=params.membrane_thickness_nm, qc_params=qc_params, ) per_file_results = batchutil.run_batch(seg_files, worker=worker, desc="Segmentation files processed", max_workers=params.max_workers) @@ -88,12 +89,18 @@ def analyse( # ========================= # DEFINE FUNCTION: processSegmentation # ========================= -def processSegmentation(seg_path: Path, fill_threshold, qc_params): +def processSegmentation(seg_path: Path, fill_threshold, membrane_thickness_nm, qc_params): ''' Process a given labelled segmentation file by calling the component processing function for each component. ''' lg.info(f"analyse | {seg_path.name} | Started processing segmentation file.") + # A31: a fitted MRC (solid geometric primitives written by `model`) is easy to feed + # in by accident — both live as sibling .mrc files under evaluator/ — and would be + # measured as if it were a real segmentation with no marker distinguishing the two. + if mrcutil.FITTED_MRC_LABEL in mrcutil.readMRCLabels(seg_path): + lg.warning(f"analyse | {seg_path.name} | Input file has been generated by `evaluator model` - skipping.") + return [] lg.debug(f"analyse | {seg_path.name} | Reading file...") data, voxel_size_nm = mrcutil.readMRCFile(seg_path) # Support both binary segmentations and pre-labelled MRC files from `label` @@ -111,6 +118,19 @@ def processSegmentation(seg_path: Path, fill_threshold, qc_params): lg.warning(f"analyse | {seg_path.name} | No components identified - skipping file.") return [] lg.info(f"analyse | {seg_path.name} | {n_components} components identified for analysis.") + # resolve thickness once per file: a number (config or CLI) is used as given, 'auto' is estimated from this file's mask + if membrane_thickness_nm == 'auto': + lg.debug(f"analyse | {seg_path.name} | Estimating membrane thickness from segmentation...") + thickness_nm = measurement.estimateMembraneThickness(components > 0, voxel_size_nm) + if thickness_nm is None: + thickness_nm, thickness_source = measurement.DEFAULT_MEMBRANE_THICKNESS_NM, 'fallback' + lg.warning(f"analyse | {seg_path.name} | Could not estimate membrane thickness — using fallback of {thickness_nm} nm.") + else: + thickness_source = 'estimated' + lg.info(f"analyse | {seg_path.name} | Membrane thickness ({thickness_source}): {thickness_nm:.2f}") + else: + thickness_nm, thickness_source = membrane_thickness_nm, 'override' + lg.debug(f"analyse | {seg_path.name} | Membrane thickness ({thickness_source}): {thickness_nm:.2f}") lg.debug(f"analyse | {seg_path.name} | Measuring component properties...") component_list = measure.regionprops(components) file_results = [] @@ -118,7 +138,7 @@ def processSegmentation(seg_path: Path, fill_threshold, qc_params): with logging_redirect_tqdm(): for component in tqdm(component_list, desc="Components processed"): lg.debug(f"analyse | {seg_path.name} | Component {component.label} | Measuring component features...") - component_data = processComponent(component.label, components, component, voxel_size_nm, seg_path.name, fill_threshold, qc_params) + component_data = processComponent(component.label, components, component, voxel_size_nm, seg_path.name, fill_threshold, thickness_nm, thickness_source, qc_params) if component_data is None: lg.warning(f"analyse | {seg_path.name} | Component {component.label} | Component processing failed — skipping.") continue @@ -131,22 +151,23 @@ def processSegmentation(seg_path: Path, fill_threshold, qc_params): # ========================= # DEFINE FUNCTION: processComponent # ========================= -def processComponent(component_label, labelled_volumes, component_properties, voxel_size_nm, filename, fill_threshold, qc_params): +def processComponent(component_label, labelled_volumes, component_properties, voxel_size_nm, filename, fill_threshold, membrane_thickness_nm, membrane_thickness_source, qc_params): ''' For a given component, make all defined measurements and return as a dictionary. ''' - lg.debug(f"analyse | {filename} | Component {component_label} | Setting scale...") scale = voxel_size_nm if voxel_size_nm is not None else 1.0 scale_label = "nm" if voxel_size_nm is not None else "vox" + lg.debug(f"analyse | {filename} | Component {component_label} | Scale set to {scale}{scale_label}") lg.debug(f"analyse | {filename} | Component {component_label} | Creating component mask...") component_mask = measurement.createComponentMask(component=component_properties, labelled_vol=labelled_volumes, label_val=component_label) lg.debug(f"analyse | {filename} | Component {component_label} | Measuring membrane volume and equivalent diameter...") - membrane_vol_nm3, equiv_diameter_nm = measurement.measureMembraneVolumeDiameter(component=component_properties, scale=scale) + membrane_vol_nm3, equiv_diameter_nm = measurement.measureMembraneVolumeDiameter(component=component_properties, scale=scale, membrane_thickness_nm=membrane_thickness_nm) + xy_diameter_nm = measurement.measureXYDiameter(component_mask, scale) component_mask_dilated = filtering.morphologicalDilation(component_mask) lg.debug(f"analyse | {filename} | Component {component_label} | Checking if component is enclosed...") enclosed, fill_ratio = filtering.checkEnclosed(component_mask=component_mask_dilated, threshold=fill_threshold) lg.debug(f"analyse | {filename} | Component {component_label} | Measuring lumen volume...") - lumen_vol_nm3 = measurement.measureLumenVolume(component_mask=component_mask, scale=scale) + lumen_vol_nm3 = measurement.measureLumenVolume(component_mask=component_mask_dilated, scale=scale) lg.debug(f"analyse | {filename} | Component {component_label} | Measuring surface area...") surface_area = measurement.computeSurfaceArea(component_mask, voxel_size_nm) lg.debug(f"analyse | {filename} | Component {component_label} | Measuring major/minor axes diameters...") @@ -161,11 +182,13 @@ def processComponent(component_label, labelled_volumes, component_properties, vo arc_coverage = geometry.arcCoverage(coords, _max_fit_points) solidity = geometry.solidity(coords, n_voxels, _max_fit_points) bbox_extent = float(component_properties.extent) + touches_z_boundary = measurement.touchesZBoundary(component_properties, labelled_volumes.shape) is_vesicle_like, qc_flags = measurement.classifyVesicle(n_voxels, sphere_rmse_rel, aspect_ratio, solidity, arc_coverage, enclosed, qc_params) return { "tomogram": filename, "label": component_label, "equiv_diameter_nm": round(equiv_diameter_nm, 2), + "xy_diameter_nm": round(xy_diameter_nm, 2), "major_axis_diameter": round(major_axis_diameter, 2), "minor_axis_diameter": round(minor_axis_diameter, 2), "aspect_ratio": round(aspect_ratio, 2), @@ -179,8 +202,11 @@ def processComponent(component_label, labelled_volumes, component_properties, vo "solidity": _r(solidity, 4), "arc_coverage": _r(arc_coverage, 4), "bbox_extent": _r(bbox_extent, 4), + "touches_z_boundary": bool(touches_z_boundary), "is_vesicle_like": is_vesicle_like, "qc_flags": qc_flags, + "membrane_thickness_nm": round(membrane_thickness_nm, 2), + "membrane_thickness_source": membrane_thickness_source, "voxel_size_nm": round(scale, 4) if voxel_size_nm is not None else None, "measurement_units": scale_label, } diff --git a/src/evaluator/commands/analyse/cli.py b/src/evaluator/commands/analyse/cli.py index cf2900f..cffc5cc 100644 --- a/src/evaluator/commands/analyse/cli.py +++ b/src/evaluator/commands/analyse/cli.py @@ -51,6 +51,10 @@ def analyse( Optional[float], typer.Option('--fill-threshold', help='Override configuration fill threshold parameter for this run', min=0, max=1) ] = None, + membrane_thickness_nm: Annotated[ + Optional[float], + typer.Option('--membrane-thickness-nm', help='Override membrane shell thickness (nm) used for diameter derivation, instead of estimating from file', min=0) + ] = None, qc_max_sphere_rmse_rel: Annotated[ Optional[float], typer.Option('--qc-max-sphere-rmse-rel', help='Max best-fit-sphere relative RMSE for a vesicle-like component', min=0, rich_help_panel='Vesicle-Like Check Overrides') @@ -83,6 +87,7 @@ def analyse( input, output, fill_threshold=fill_threshold, + membrane_thickness_nm=membrane_thickness_nm, qc_max_sphere_rmse_rel=qc_max_sphere_rmse_rel, qc_max_aspect_ratio=qc_max_aspect_ratio, qc_min_solidity=qc_min_solidity, diff --git a/src/evaluator/commands/analyse/utils/filtering.py b/src/evaluator/commands/analyse/utils/filtering.py index 35d3872..b6b8f7f 100644 --- a/src/evaluator/commands/analyse/utils/filtering.py +++ b/src/evaluator/commands/analyse/utils/filtering.py @@ -22,10 +22,14 @@ # ========================= def morphologicalDilation(binary_vol: numpy.ndarray): ''' - Applies morphological dilation to bridge small gaps in thin membrane shells. + Applies morphological dilation to bridge small gaps in thin membrane shells. Use over morphologicalClosure to as erosion can remove gap-filling voxels added by dilation. + To avoid outward asymmetric clipping from bounding-box input mask, pad volume by structuring element's radius then crop back to original shape. ''' - return ndimage.binary_dilation(binary_vol, structure=ball(2)) + pad = 2 + padded = numpy.pad(binary_vol, pad_width=pad, mode='constant', constant_values=False) + dilated = ndimage.binary_dilation(padded, structure=ball(2)) + return dilated[pad:-pad, pad:-pad, pad:-pad] # ========================= # DEFINE FUNCTION: checkEnclosed diff --git a/src/evaluator/commands/analyse/utils/geometry.py b/src/evaluator/commands/analyse/utils/geometry.py index aed9683..354a6d9 100644 --- a/src/evaluator/commands/analyse/utils/geometry.py +++ b/src/evaluator/commands/analyse/utils/geometry.py @@ -74,26 +74,32 @@ def arcCoverage(coords: numpy.ndarray, max_fit_points: int = _QC_MAX_FIT_POINTS) # ========================= def deriveAxes(intertia_tensor, voxel_size_nm=None): ''' - Derive semi-axes (a ≥ b ≥ c) of the best-fit ellipsoid from inertia tensor eigenvalues. eigvalsh returns eigenvalues in ascending order (I_a ≤ I_b ≤ I_c). + Derive semi-axes (a ≥ b ≥ c) of the best-fit ellipsoid from inertia tensor eigenvalues. + eigvalsh returns eigenvalues in ascending order (I_a ≤ I_b ≤ I_c). + For ellipsoid with principal moments I_a=(m/5)(b^2+c^2), use inversion a^2 ∝ (I_b + I_c - I_a) instead of 1/sqrt(I_a). ''' eigvals = numpy.linalg.eigvalsh(intertia_tensor) eigvals = numpy.clip(eigvals, 0, None) - with numpy.errstate(divide="ignore", invalid="ignore"): - inv_sqrt = numpy.where(eigvals > 0, 1.0 / numpy.sqrt(eigvals), 0.0) - return inv_sqrt + ia, ib, ic = eigvals + # index 0 pairs with the smallest eigenvalue (ia) -> longest axis; index 2 pairs with ic -> shortest axis + sq_lengths = numpy.clip([ib + ic - ia, ia + ic - ib, ia + ib - ic], 0, None) + return numpy.sqrt(sq_lengths) # ========================= # DEFINE FUNCTION: measureAxes # ========================= def measureAxes(component, equiv_diameter_nm): ''' - Measure major and minor axes by approximating the EV as an ellipsoid. Semi-axes are derived from inertia tensor eigenvalues and scaled to real-world size using the equivalent diameter. + Measure major and minor axes by approximating the EV as an ellipsoid. + Relative semi-axes are derived from inertia tensor eigenvalues and rescaled so geometric mean matches the calculated equivalent diameter. + Degenerate components collapse to 0 and return NaN. ''' - inv_sqrt_axes = deriveAxes(intertia_tensor=component.inertia_tensor) - geomean_inv_sqrt = (inv_sqrt_axes[0] * inv_sqrt_axes[1] * inv_sqrt_axes[2]) ** (1 / 3) - if geomean_inv_sqrt > 0: - axis_scale = (equiv_diameter_nm / 2.0) / geomean_inv_sqrt - principal_semiaxis_a, principal_semiaxis_b, principal_semiaxis_c = inv_sqrt_axes * axis_scale + axis_lengths = deriveAxes(intertia_tensor=component.inertia_tensor) + geomean = (axis_lengths[0] * axis_lengths[1] * axis_lengths[2]) ** (1 / 3) + if not (geomean > 0) or not numpy.isfinite(geomean): + return numpy.nan, numpy.nan + axis_scale = (equiv_diameter_nm / 2.0) / geomean + principal_semiaxis_a, principal_semiaxis_b, principal_semiaxis_c = axis_lengths * axis_scale major_axis = 2 * principal_semiaxis_a minor_axis = 2 * principal_semiaxis_c return major_axis, minor_axis diff --git a/src/evaluator/commands/analyse/utils/io.py b/src/evaluator/commands/analyse/utils/io.py index a66ddab..32b54cc 100644 --- a/src/evaluator/commands/analyse/utils/io.py +++ b/src/evaluator/commands/analyse/utils/io.py @@ -40,5 +40,5 @@ def printSummaryMessage(results, nfiles: int, startt: datetime.datetime, endt: d print(f"- Number of enclosed EVs: {results['is_enclosed'].sum()} ({100 * results['is_enclosed'].mean():.1f}%)") if 'is_vesicle_like' in results: print(f"- Vesicle-like components: {results['is_vesicle_like'].sum()} / {len(results)} ({100 * results['is_vesicle_like'].mean():.1f}%)") - print(f"- Equivalent diameters: {results['equiv_diameter_nm'].mean():.1f} ± {results['equiv_diameter_nm'].std():.1f} nm (mean ± SD)") + print(f"- Outer diameters: {results['equiv_diameter_nm'].mean():.1f} ± {results['equiv_diameter_nm'].std():.1f} nm (mean ± SD)") print(f"Results saved to: {out_path}\n") \ No newline at end of file diff --git a/src/evaluator/commands/analyse/utils/measurement.py b/src/evaluator/commands/analyse/utils/measurement.py index 0f0b312..cd64166 100644 --- a/src/evaluator/commands/analyse/utils/measurement.py +++ b/src/evaluator/commands/analyse/utils/measurement.py @@ -11,6 +11,9 @@ from scipy import ndimage from skimage import measure +# DEFAULT_MEMBRANE_THICKNESS_NM: measured median membrane segmentation thickness produced by MemBrain-seg +DEFAULT_MEMBRANE_THICKNESS_NM = 2.6 + # ========================= # DEFINE FUNCTION: computeSurfaceArea # ========================= @@ -31,18 +34,68 @@ def computeSurfaceArea(component_mask: numpy.ndarray, voxel_size_nm: float): return sa_vox * (voxel_size_nm ** 2) return sa_vox +# ========================= +# DEFINE FUNCTION: estimateMembraneThickness +# ========================= +def estimateMembraneThickness(mask: numpy.ndarray, voxel_size_nm: float | None, *, min_ridge_voxels: int = 100) -> float | None: + ''' + Estimate segmented membrane shell thickness (twice the median Euclidean distance-transform value along local maxima). + ''' + if not mask.any(): + return None + # crop to the mask's bounding box to keep the distance transform small + mask = mask[ndimage.find_objects(mask.astype(numpy.uint8))[0]] + dist = ndimage.distance_transform_edt(mask, sampling=voxel_size_nm if voxel_size_nm is not None else 1.0) + # plateau-tolerant ridge: voxels at least as far from the background as all 26 neighbours + ridge = mask & (dist >= ndimage.maximum_filter(dist, size=3)) + if ridge.sum() < min_ridge_voxels: + return None + return float(2 * numpy.median(dist[ridge])) + + # ========================= # DEFINE FUNCTION: measureMembraneVolumeDiameter # ========================= -def measureMembraneVolumeDiameter(component, scale): +def measureMembraneVolumeDiameter(component, scale, membrane_thickness_nm: float = 2.6): ''' - Measures the volume of membrane components in voxels and converts to nm^3. Calculates the equivalent spherical diameter in nm. + Measures the volume of membrane components in voxels and converts to nm^3. Derives the equivalent spherical diameter by treating membrane volume as a spherical shell of known thickness. ''' vol_vox = component.area vol_nm3 = vol_vox * (scale ** 3) - equiv_diameter_nm = (6 * vol_nm3 / numpy.pi) ** (1 / 3) + t = membrane_thickness_nm + if t is not None and t > 0: + # use the positive root from solving 4/3 pi (R^3 - (R-t)^3) = vol_nm3 + discriminant = max(vol_nm3 / (numpy.pi * t) - (t ** 2) / 3, 0.0) + outer_radius_nm = (t / 2) + (numpy.sqrt(discriminant) / 2) + equiv_diameter_nm = 2 * outer_radius_nm + else: + # no thickness available: fall back to a biased solid-sphere-equivalent diameter + equiv_diameter_nm = (6 * vol_nm3 / numpy.pi) ** (1 / 3) return vol_nm3, equiv_diameter_nm +# ========================= +# DEFINE FUNCTION: measureXYDiameter +# ========================= +def measureXYDiameter(component_mask, scale): + ''' + Diameter from the widest xy slice (mean of the y and x extents), in nm if scale is a voxel size. + Independent of z, so unaffected by clipping at the z faces. + ''' + per_slice = (_sliceExtent(component_mask.any(axis=2)) + _sliceExtent(component_mask.any(axis=1))) / 2 + return float(per_slice.max() * scale) + +# ========================= +# DEFINE HELPER FUNCTION: _sliceExtent +# ========================= +def _sliceExtent(slice_mask: numpy.ndarray) -> numpy.ndarray: + ''' + Per-row extent (voxels) between the first and last True along axis 1; 0 for empty rows. + ''' + idx = numpy.arange(slice_mask.shape[1]) + hi = numpy.where(slice_mask, idx, -1).max(axis=1) + lo = numpy.where(slice_mask, idx, slice_mask.shape[1]).min(axis=1) + return numpy.where(slice_mask.any(axis=1), hi - lo + 1, 0) + # ========================= # DEFINE FUNCTION: createComponentMask # ========================= @@ -88,7 +141,8 @@ def measureEccentricityAspectRatio(major_axis_diameter, minor_axis_diameter): # ========================= def classifyVesicle(n_voxels, sphere_rmse_rel, aspect_ratio, solidity, arc_coverage, is_enclosed, qc_params): ''' - Classify component as `is_vesicle_like` unless any enabled check fails''' + Classify component as `is_vesicle_like` unless any enabled check fails. + ''' if n_voxels < 4: return False, 'too_small' flags = [] @@ -101,3 +155,13 @@ def classifyVesicle(n_voxels, sphere_rmse_rel, aspect_ratio, solidity, arc_cover if not (is_enclosed or arc_coverage >= qc_params['qc_min_arc_coverage']): flags.append('open_shell') return len(flags) == 0, ','.join(flags) + +# ========================= +# DEFINE FUNCTION: touchesZBoundary +# ========================= +def touchesZBoundary(component_properties, volume_shape) -> bool: + ''' + If a component's bounding box reaches the Z boundary, it may be truncated by the reconstruction box. + ''' + z_min, z_max = component_properties.bbox[0], component_properties.bbox[3] + return z_min == 0 or z_max == volume_shape[0] diff --git a/src/evaluator/commands/label/label.py b/src/evaluator/commands/label/label.py index b8da316..bb2e271 100644 --- a/src/evaluator/commands/label/label.py +++ b/src/evaluator/commands/label/label.py @@ -24,7 +24,7 @@ # Import EValuator label utilities # ==================== from evaluator.commands.label.utils.filters import filterComponentsBySize -from evaluator.commands.label.utils.geometric_proxies import estimateCentroidRadius, estimateArcCoverage +from evaluator.commands.label.utils.geometric_proxies import estimateCentroidRadius, estimateArcCoverage, ARC_COVERAGE_N_BINS from evaluator.commands.label.utils.merge import findMergeGroups, applyMerges # ==================== @@ -64,6 +64,7 @@ def label_components( # Label components lg.info(f"label | Labelling connected components...") labelled, n_components = mrcutil.labelComponents(seg_data) + n_initial = n_components lg.info(f"label | {n_components} components identified.") # Step 1: merge split components (geometric proxies only) @@ -73,20 +74,34 @@ def label_components( component_points, centre_tol_factor=params.merge_centre_tol_factor, radius_tol_pct=params.merge_radius_tol_pct, + concentric_radius_pct=params.merge_concentric_radius_pct, + concentric_centre_frac=params.merge_concentric_centre_frac, ) labelled = applyMerges(labelled, merge_groups) + n_merge_groups_merged = sum(1 for group in merge_groups if len(group) > 1) + n_after_merge = len(numpy.unique(labelled)) - 1 + # Step 2: arc-coverage filter (post-merge) lg.debug(f"label | Applying arc-coverage filter...") + n_excluded_arc_coverage = 0 + n_excluded_arc_coverage_size_floor = 0 for component in measure.regionprops(labelled): centroid, radius_estimate = estimateCentroidRadius(component.coords) coverage = estimateArcCoverage(component.coords, centroid, radius_estimate) if coverage < params.min_arc_coverage: - lg.info(f"label | Component {component.label} | Arc coverage {coverage:.2f} below threshold ({params.min_arc_coverage}) — excluding.") + max_possible_coverage = min(component.area, ARC_COVERAGE_N_BINS) / ARC_COVERAGE_N_BINS + size_floor = max_possible_coverage < params.min_arc_coverage + reason = "size floor (too few voxels to reach threshold)" if size_floor else "incomplete coverage" + lg.info(f"label | Component {component.label} | Arc coverage {coverage:.2f} below threshold ({params.min_arc_coverage}) — excluding ({reason}).") + n_excluded_arc_coverage += 1 + if size_floor: + n_excluded_arc_coverage_size_floor += 1 labelled[labelled == component.label] = 0 # Step 3: size/extent filter lg.debug(f"label | Applying size filter...") + n_before_size_filter = len(numpy.unique(labelled)) - 1 labelled = filterComponentsBySize( labelled, voxel_size_nm, @@ -94,9 +109,14 @@ def label_components( maximum_diameter_nm=params.maximum_diameter_nm, membrane_thickness_nm=params.membrane_thickness_nm, ) + n_excluded_size = n_before_size_filter - (len(numpy.unique(labelled)) - 1) - # Relabel sequentially so downstream commands see a 1..N label range - labelled, n_components = mrcutil.labelComponents(labelled.astype(bool)) + surviving_labels = sorted(int(l) for l in numpy.unique(labelled) if l != 0) + remapped = numpy.zeros_like(labelled) + for new_id, old_id in enumerate(surviving_labels, start=1): + remapped[labelled == old_id] = new_id + labelled = remapped + n_components = len(surviving_labels) lg.info(f"label | {n_components} components retained after merge/filter.") # Build output path @@ -116,4 +136,22 @@ def label_components( # Write labelled MRC lg.debug(f"label | Writing labelled MRC to {out_file.name}...") mrcutil.writeMRCFile(labelled.astype(numpy.float32), voxel_size_nm, out_file) - lg.info(f"label | Finished labelling.") \ No newline at end of file + lg.info(f"label | Finished labelling.") + + # Write per-file exclusion tally (visible at warning for non-interactive batch runs) + tally = ( + f"components identified: {n_initial}\n" + f"merge groups merged (>1 component): {n_merge_groups_merged}\n" + f"components after merge: {n_after_merge}\n" + f"excluded by arc-coverage filter: {n_excluded_arc_coverage}\n" + f"excluded by size/extent filter: {n_excluded_size}\n" + f"components retained: {n_components}\n" + ) + tally_file = out_file.with_name(f"{out_file.stem}_exclusions.txt") + tally_file.write_text(tally) + lg.warning( + f"label | {segmentation.name}: {n_initial} identified -> {n_components} retained; " + f"{n_excluded_arc_coverage} excluded by arc-coverage filter ({n_excluded_arc_coverage_size_floor} too small to reach threshold), " + f"{n_excluded_size} excluded by size filter; " + f"tally written to {tally_file.name}." + ) diff --git a/src/evaluator/commands/label/utils/geometric_proxies.py b/src/evaluator/commands/label/utils/geometric_proxies.py index 2c0bfd4..02d3114 100644 --- a/src/evaluator/commands/label/utils/geometric_proxies.py +++ b/src/evaluator/commands/label/utils/geometric_proxies.py @@ -9,6 +9,13 @@ # ==================== import numpy +# ========================= +# DEFINE CONSTANT: arc-coverage angular grid +# ========================= +ARC_COVERAGE_N_LAT_BINS, ARC_COVERAGE_N_LON_BINS = 18, 36 +ARC_COVERAGE_N_BINS = ARC_COVERAGE_N_LAT_BINS * ARC_COVERAGE_N_LON_BINS + + # ========================= # DEFINE FUNCTION: estimateCentroidRadius # ========================= @@ -69,9 +76,9 @@ def estimateArcCoverage(points: numpy.ndarray, centroid: numpy.ndarray, radius_e norms[norms == 0] = 1e-9 unit_vecs = vecs / norms[:, None] # Coarse lat/long grid (10-degree bins each), equal-area: bin sin(lat) - n_lat_bins, n_lon_bins = 18, 36 + n_lat_bins, n_lon_bins = ARC_COVERAGE_N_LAT_BINS, ARC_COVERAGE_N_LON_BINS lon = numpy.degrees(numpy.arctan2(unit_vecs[:, 1], unit_vecs[:, 0])) lat_idx = numpy.clip((((unit_vecs[:, 2] + 1) / 2) * n_lat_bins).astype(int), 0, n_lat_bins - 1) lon_idx = numpy.clip(((lon + 180) / 360 * n_lon_bins).astype(int), 0, n_lon_bins - 1) occupied = set(zip(lat_idx.tolist(), lon_idx.tolist())) - return len(occupied) / (n_lat_bins * n_lon_bins) + return len(occupied) / ARC_COVERAGE_N_BINS diff --git a/src/evaluator/commands/label/utils/merge.py b/src/evaluator/commands/label/utils/merge.py index 4a06c1d..16feb73 100644 --- a/src/evaluator/commands/label/utils/merge.py +++ b/src/evaluator/commands/label/utils/merge.py @@ -21,6 +21,8 @@ def findMergeGroups( component_points: dict[int, numpy.ndarray], centre_tol_factor: float, radius_tol_pct: float, + concentric_radius_pct: float = 0.10, + concentric_centre_frac: float = 0.25, ) -> list[list[int]]: ''' Returns groups of label_ids that should be merged, based on centroid proximity and radius-estimate consistency using union-find @@ -44,7 +46,7 @@ def union(a: int, b: int) -> None: centroid_a, radius_a = proxies[id_a] for id_b in label_ids[i + 1:]: centroid_b, radius_b = proxies[id_b] - if _shouldMerge(centroid_a, radius_a, centroid_b, radius_b, centre_tol_factor, radius_tol_pct): + if _shouldMerge(centroid_a, radius_a, centroid_b, radius_b, centre_tol_factor, radius_tol_pct, concentric_radius_pct, concentric_centre_frac): union(id_a, id_b) groups: dict[int, list[int]] = {} @@ -55,15 +57,18 @@ def union(a: int, b: int) -> None: # ========================= # DEFINE FUNCTION: _shouldMerge # ========================= -def _shouldMerge(centroid_a, radius_a, centroid_b, radius_b, centre_tol_factor, radius_tol_pct) -> bool: - '''Mutually-consistent check: centroid proximity and radius agreement.''' +def _shouldMerge(centroid_a, radius_a, centroid_b, radius_b, centre_tol_factor, radius_tol_pct, concentric_radius_pct=0.10, concentric_centre_frac=0.25) -> bool: + '''Mutually-consistent check: centroid proximity and radius agreement, unless the pair is concentric.''' centre_distance = float(numpy.linalg.norm(centroid_a - centroid_b)) centre_ok = centre_distance <= centre_tol_factor * (radius_a + radius_b) if max(radius_a, radius_b) == 0: - radius_ok = False - else: - radius_ok = abs(radius_a - radius_b) / max(radius_a, radius_b) <= radius_tol_pct - return centre_ok and radius_ok + return False + radius_diff = abs(radius_a - radius_b) / max(radius_a, radius_b) + if not (centre_ok and radius_diff <= radius_tol_pct): + return False + # concentric signature: near-identical centre but clearly different radius; split fragments of one sphere share both + concentric = centre_distance <= concentric_centre_frac * min(radius_a, radius_b) and radius_diff > concentric_radius_pct + return not concentric # ========================= # DEFINE FUNCTION: applyMerges diff --git a/src/evaluator/commands/model/model.py b/src/evaluator/commands/model/model.py index 7fd9ec7..17f4336 100644 --- a/src/evaluator/commands/model/model.py +++ b/src/evaluator/commands/model/model.py @@ -35,7 +35,7 @@ def model_batch(input_path, output_dir, max_workers=None, **overrides) -> None: ''' mrc_files = batchutil.resolve_mrc_inputs(input_path) config, _ = confutil.load_config(output_dir) - max_workers = max_workers if max_workers is not None else config.label.max_workers + max_workers = max_workers if max_workers is not None else config.model.max_workers worker = partial(model_evs, output_dir=output_dir, **overrides) batchutil.run_batch(mrc_files, worker=worker, desc="Labelled files processed", max_workers=max_workers) @@ -56,6 +56,9 @@ def model_evs( labelled_data, voxel_size_nm = mrcutil.readMRCFile(input_file) labels = np.unique(labelled_data) labels = labels[labels != 0] + # Check if only one non-zero label found and warn if so + if len(labels) == 1 and len(np.unique(labelled_data)) == 2: + lg.warning(f"{input_file.name}: only one non-zero label found. Ensure this is a labelled tomogram (i.e. run through evaluator label) otherwise all segmentations will be fitted as a single object.") records: list[dict] = [] # Fit each labelled vesicle lg.info("model | Fitting {} labelled components...", len(labels)) @@ -101,5 +104,5 @@ def model_evs( voxel_size_nm=voxel_size_nm, ) mrc_out_file = out_dir / f"{input_file.stem}_model_fitted.mrc" - mrcutil.writeMRCFile(fitted_volume.astype(np.uint16), voxel_size_nm, mrc_out_file) + mrcutil.writeMRCFile(fitted_volume.astype(np.uint16), voxel_size_nm, mrc_out_file, label=mrcutil.FITTED_MRC_LABEL) lg.debug(f"model | Wrote fitted MRC file") diff --git a/src/evaluator/commands/model/utils/least_squares_fit.py b/src/evaluator/commands/model/utils/least_squares_fit.py index 98f1d1f..ee8fbe5 100644 --- a/src/evaluator/commands/model/utils/least_squares_fit.py +++ b/src/evaluator/commands/model/utils/least_squares_fit.py @@ -106,6 +106,8 @@ def _ellipsoid_geometric_residuals( return np.abs(norm_dist - 1.0) * np.mean(radii) # _bic: returns a float for the calculated Bayesian Information Criterion for an array of residuals +# N.B. sphere residuals are exact orthogonal distances, but ellipsoid residuals are normalised algebraic distance rescaled by the mean radius. +# As the two residual sets use different scales, bic_ellipsoid < bic_sphere is not a like-for-like selection test but more of a heuristic. def _bic(residuals: np.ndarray, k: int) -> float: n = len(residuals) rss = np.sum(residuals ** 2) @@ -129,9 +131,11 @@ def _assess_reliability( vecs = points - centre norms = np.linalg.norm(vecs, axis=1) with np.errstate(invalid='ignore'): - latitudes = np.degrees(np.arcsin(np.clip(vecs[:, 2] / norms, -1, 1))) + latitudes = np.degrees(np.arcsin(np.clip(vecs[:, 0] / norms, -1, 1))) lat_span = float(latitudes.max() - latitudes.min()) if n > 0 else 0.0 span_ok = lat_span >= min_latitude_span_deg + # size_limited: rmse_ok is the only failing criterion + size_limited = (not rmse_ok) and count_ok and span_ok return { 'is_reliable': rmse_ok and count_ok and span_ok, 'relative_rmse': float(relative_rmse), @@ -139,6 +143,7 @@ def _assess_reliability( 'count_ok': count_ok, 'lat_span_deg': lat_span, 'span_ok': span_ok, + 'size_limited': size_limited, } # _beam_axis_flag: identifies if an ellipsoid's major axis is close to the z-axis, returning a dictionary containing this check and the angles @@ -146,7 +151,7 @@ def _beam_axis_flag(evecs: np.ndarray, radii: np.ndarray, tol_deg: float) -> dic major_idx = np.argmax(radii) major_axis = evecs[:, major_idx] # Angle between major axis and z - cos_angle = np.clip(np.abs(major_axis[2]), 0, 1) + cos_angle = np.clip(np.abs(major_axis[0]), 0, 1) angle_from_z = float(np.degrees(np.arccos(cos_angle))) flagged = angle_from_z < tol_deg return { diff --git a/src/evaluator/commands/plot/cli.py b/src/evaluator/commands/plot/cli.py index db071e8..7e7dba6 100644 --- a/src/evaluator/commands/plot/cli.py +++ b/src/evaluator/commands/plot/cli.py @@ -15,6 +15,7 @@ # Import EValuator utilities # ==================== from evaluator.commands.plot import plot as plotFuncs +from evaluator.commands.plot.utils.input import parse_agreement # ==================== # Initialise typer as evaluatorPlot @@ -23,7 +24,7 @@ add_completion=False, ) -VALID_SECTIONS = ["distributions", "qc", "scatter", "concordance", "compare"] +VALID_SECTIONS = ["distributions", "qc", "scatter", "correlations", "concordance", "reliability", "agreement", "compare"] @evaluatorPlot.command(help='Generate plots from evaluator analyse and/or model output', rich_help_panel='Component Analysis') def plot( @@ -51,9 +52,13 @@ def plot( bool, typer.Option('--overwrite', help='Overwrite existing section outputs instead of skipping them') ] = False, - rscript: Annotated[ - Optional[Path], - typer.Option('--rscript', help='Path to the Rscript binary (default: resolved from PATH)') + include_flagged: Annotated[ + bool, + typer.Option('--include-flagged', help='Include components flagged as non-vesicle-like instead of excluding them') + ] = False, + agreement: Annotated[ + Optional[list[str]], + typer.Option('--agreement', help='Bland-Altman plot for a pair of measurements (repeatable)') ] = None, ): ''' @@ -61,6 +66,10 @@ def plot( ''' if analyse_input is None and model_input is None: raise typer.BadParameter('At least one of --analyse or --model must be provided.') + try: + agreement_pairs = parse_agreement(agreement or []) + except ValueError as e: + raise typer.BadParameter(str(e)) plotFuncs.run_plot( analyse_input, model_input, @@ -68,5 +77,6 @@ def plot( sections=section, all_sections=all_sections, overwrite=overwrite, - rscript=rscript, + include_flagged=include_flagged, + agreement=agreement_pairs, ) diff --git a/src/evaluator/commands/plot/plot.py b/src/evaluator/commands/plot/plot.py index 7b7b097..ba63875 100644 --- a/src/evaluator/commands/plot/plot.py +++ b/src/evaluator/commands/plot/plot.py @@ -14,30 +14,24 @@ # ==================== from evaluator.utils import config as confutil from evaluator.utils import paths as pathutil +from evaluator.utils.figures import plots as plotutil +from evaluator.utils.figures import theme as themeutil from evaluator.utils.settings import lg +from evaluator.commands.plot import sections as sectionutil +from evaluator.commands.plot.utils import export as exportutil from evaluator.commands.plot.utils.input import resolve_plot_inputs, available_sections, PlotRun -from evaluator.commands.plot.utils.dispatch import resolve_rscript, dispatch, RscriptError - -# ==================== -# Define R scripts for each section -# ==================== -SECTION_SCRIPTS = { - 'distributions': 'distributions.R', - 'qc': 'qc.R', - 'scatter': 'scatter.R', - 'concordance': 'concordance.R', - 'compare': 'compare.R', -} # ==================== # run_plot: orchestration logic for plot command # ==================== -def run_plot(analyse_input, model_input, output, sections, all_sections, overwrite, rscript): +def run_plot(analyse_input, model_input, output, sections, all_sections, overwrite, include_flagged=False, agreement=None): config, evaluator_dir = confutil.load_config(output) params = config.plot runs, multi_run, sheet_path = resolve_plot_inputs(analyse_input, model_input) runnable = available_sections(runs, multi_run) requested = runnable if all_sections else (sections or params.default_sections) + if agreement and 'agreement' not in sections: + requested.append('agreement') to_run = [s for s in requested if s in runnable] skipped = set(requested) - set(to_run) if skipped: @@ -45,7 +39,12 @@ def run_plot(analyse_input, model_input, output, sections, all_sections, overwri if not to_run: lg.error('plot | No runnable sections for the given input.') return - rscript_bin = resolve_rscript(rscript) + try: + exportutil.check_kaleido() + except RuntimeError as e: + lg.error(str(e)) + return + plotutil.use_theme(themeutil.resolve(params.theme)) out_dir = pathutil.generate_command_output_dir(evaluator_dir, 'plot') confutil.write_params(params, out_dir) completed = [] @@ -56,30 +55,25 @@ def run_plot(analyse_input, model_input, output, sections, all_sections, overwri continue section_dir.mkdir(parents=True, exist_ok=True) try: - _run_section(section, runs, multi_run, section_dir, rscript_bin, params, sheet_path) + _run_section(section, runs, multi_run, section_dir, params, include_flagged, agreement) completed.append(section) - except RscriptError as e: - lg.warning('plot | Section {} failed: {}', section, e) + except Exception as e: + lg.warning(f'plot | Section {section} failed: {e}') _write_index(out_dir, completed) lg.info(f'plot | Finished. Sections completed: {completed}') # ==================== # _run_section: orchestration logic for each section # ==================== -def _run_section(section: str, runs, multi_run: bool, section_dir: Path, rscript_bin: Path, params, sheet_path: Path | None = None) -> None: - mode = 'multi' if multi_run else 'single' - script = SECTION_SCRIPTS[section] - if section == 'concordance': - run = runs[0] - dispatch(rscript_bin, script, [section_dir, run.analyse_path, run.model_path, mode]) - elif section == 'qc': - run = runs[0] - dispatch(rscript_bin, script, [section_dir, run.analyse_path, mode, params.fill_ratio_flag_threshold]) - elif section in ('distributions', 'scatter'): - run = runs[0] - dispatch(rscript_bin, script, [section_dir, run.analyse_path, mode]) - elif section == 'compare': - dispatch(rscript_bin, script, [section_dir, sheet_path, mode]) +def _run_section(section: str, runs: list[PlotRun], multi_run: bool, section_dir: Path, params, include_flagged: bool, agreement) -> None: + if section == 'compare': + sectionutil.compare(section_dir, runs, params, include_flagged=include_flagged) + return + # per-sample sections: in multi-run mode every sample is plotted into its own subdirectory + for run in (runs if multi_run else runs[:1]): + run_dir = section_dir / run.sample_id if multi_run else section_dir + run_dir.mkdir(parents=True, exist_ok=True) + sectionutil.SECTIONS[section](run_dir, sectionutil.load(run, include_flagged=include_flagged, agreement=agreement), params) # ==================== # _write_index: write the markdown index of all generated files @@ -88,4 +82,4 @@ def _write_index(out_dir: Path, sections: list[str]) -> None: lines = ['# EValuator plot results', ''] for s in sections: lines.append(f'- [`{s}/`]({s}/)') - (out_dir / 'index.md').write_text('\n'.join(lines) + '\n') \ No newline at end of file + (out_dir / 'index.md').write_text('\n'.join(lines) + '\n') diff --git a/src/evaluator/commands/plot/r/compare.R b/src/evaluator/commands/plot/r/compare.R deleted file mode 100644 index d28bccb..0000000 --- a/src/evaluator/commands/plot/r/compare.R +++ /dev/null @@ -1,55 +0,0 @@ -# compare.R: plot group comparisons - -# Parse positional arguments -args <- commandArgs(trailingOnly = TRUE) -output_dir <- args[1]; sheet_path <- args[2]; mode <- args[3] - -# Get script directory for path traversal -script_dir <- local({ - args <- commandArgs(trailingOnly = FALSE) - script <- grep("^--file=", args, value = TRUE) - dirname(normalizePath(sub("^--file=", "", script))) -}) - -# Import internal utilities -source(file.path(script_dir, "utils", "import.R")) -source(file.path(script_dir, "utils", "io.R")) -source(file.path(script_dir, "utils", "theme.R")) - -# Load data sheets -sheet <- readr::read_tsv(sheet_path, show_col_types = FALSE) -loaded <- lapply(seq_len(nrow(sheet)), function(i) { - row <- sheet[i, ] - df <- loadAnalyseCSV(row$path) - df$sample_id <- row$sample_id - df$group <- row$group - df -}) -combined <- dplyr::bind_rows(loaded) - -# Define features to plot and plot violin plot for each group -features <- c("equiv_diameter_nm", "aspect_ratio", "closure_fill_ratio", "lumen_volume") -pal <- groupPalette(combined$group) -for (feat in intersect(features, names(combined))) { - p <- ggplot2::ggplot(combined, ggplot2::aes(x = group, y = .data[[feat]], fill = group)) + - ggplot2::geom_violin(alpha = 0.6) + ggplot2::geom_boxplot(width = 0.1, outlier.shape = NA) + - ggplot2::scale_fill_manual(values = pal) + evTheme() + ggplot2::theme(legend.position = "none") - saveSVG(p, file.path(output_dir, paste0("violin_", feat, ".svg"))) -} - -# Create matrix for features and save PCA plot -feature_matrix <- combined[, intersect(features, names(combined))] -feature_matrix <- feature_matrix[stats::complete.cases(feature_matrix), ] -if (nrow(feature_matrix) >= 3 && ncol(feature_matrix) >= 2) { - pca <- stats::prcomp(scale(feature_matrix)) - pca_df <- data.frame(pca$x[, 1:2], group = combined$group[stats::complete.cases(combined[, features])]) - p_pca <- ggplot2::ggplot(pca_df, ggplot2::aes(x = PC1, y = PC2, colour = group)) + - ggplot2::geom_point(alpha = 0.7) + ggplot2::scale_colour_manual(values = pal) + evTheme() - saveSVG(p_pca, file.path(output_dir, "pca.svg"), width = 6, height = 6) -} - -# Write results to spreadsheet -writeXLSX(dplyr::summarise(dplyr::group_by(combined, group), - n_evs = dplyr::n(), - across(all_of(intersect(features, names(combined))), list(mean = ~mean(.x, na.rm=TRUE), median = ~median(.x, na.rm=TRUE)))), - file.path(output_dir, "group_summary.xlsx")) diff --git a/src/evaluator/commands/plot/r/concordance.R b/src/evaluator/commands/plot/r/concordance.R deleted file mode 100644 index b1c46b8..0000000 --- a/src/evaluator/commands/plot/r/concordance.R +++ /dev/null @@ -1,58 +0,0 @@ -# concordance.R: plot concordance between model and analyse - -# Parse positional arguments -args <- commandArgs(trailingOnly = TRUE) -output_dir <- args[1]; analyse_path <- args[2]; model_path <- args[3]; mode <- args[4] - -# Get script directory for path traversal -script_dir <- local({ - args <- commandArgs(trailingOnly = FALSE) - script <- grep("^--file=", args, value = TRUE) - dirname(normalizePath(sub("^--file=", "", script))) -}) - -# Import internal utilities -source(file.path(script_dir, "utils", "import.R")) -source(file.path(script_dir, "utils", "io.R")) -source(file.path(script_dir, "utils", "theme.R")) - -# Import data -analyse_df <- loadAnalyseCSV(analyse_path) -model_df <- loadModelResults(model_path) -joined <- joinAnalyseModel(analyse_df, model_df) -if (nrow(joined) == 0) { - stop("No matching (tomogram, label) rows between analyse and model output — check both were run on the same labelled MRC(s).") -} -joined$model_diameter_nm <- joined$radius * 2 - -# Diameter concordance: model's fitted diameter vs analyse's two independent diameter estimates, with y = x reference line. -for (analyse_diam in c("equiv_diameter_nm", "major_axis_diameter")) { - p <- ggplot2::ggplot(joined, ggplot2::aes(x = .data[[analyse_diam]], y = model_diameter_nm, colour = tomogram)) + - ggplot2::geom_point(alpha = 0.6) + - ggplot2::geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = "grey40") + - evTheme() + ggplot2::theme(legend.position = "none") + - ggplot2::labs(x = analyse_diam, y = "model_diameter_nm (2 x fitted radius)") - saveSVG(p, file.path(output_dir, paste0("model_vs_", analyse_diam, ".svg")), width = 5, height = 5) -} - -# See if model reliability tracks analyse's independent closure signal -p_reliab <- ggplot2::ggplot(joined, ggplot2::aes(x = reliability.is_reliable, y = closure_fill_ratio)) + - ggplot2::geom_boxplot(fill = "#56B4E9") + evTheme() + - ggplot2::labs(x = "model reliability.is_reliable", y = "analyse closure_fill_ratio") -saveSVG(p_reliab, file.path(output_dir, "reliability_vs_closure.svg")) - -# RMSE of the fit vs whether analyse independently judged the EV enclosed -p_rmse <- ggplot2::ggplot(joined, ggplot2::aes(x = is_enclosed, y = rmse_nm)) + - ggplot2::geom_boxplot(fill = "#D55E00") + evTheme() + - ggplot2::labs(x = "analyse is_enclosed", y = "model rmse_nm") -saveSVG(p_rmse, file.path(output_dir, "fit_rmse_vs_enclosed.svg")) - -conc_stats <- data.frame( - n_matched = nrow(joined), - n_analyse_only = nrow(analyse_df) - nrow(joined), - n_model_only = nrow(model_df) - nrow(joined), - pearson_equiv = cor(joined$equiv_diameter_nm, joined$model_diameter_nm, use = "complete.obs"), - pearson_major = cor(joined$major_axis_diameter, joined$model_diameter_nm, use = "complete.obs"), - reliable_pass_rate = mean(joined$reliability.is_reliable, na.rm = TRUE) -) -writeXLSX(conc_stats, file.path(output_dir, "concordance_stats.xlsx")) diff --git a/src/evaluator/commands/plot/r/distributions.R b/src/evaluator/commands/plot/r/distributions.R deleted file mode 100644 index a350e8e..0000000 --- a/src/evaluator/commands/plot/r/distributions.R +++ /dev/null @@ -1,40 +0,0 @@ -# distributions.R: plot distribution of EV morphometrics - -# Parse positional arguments -args <- commandArgs(trailingOnly = TRUE) -output_dir <- args[1]; analyse_path <- args[2]; mode <- args[3] - -# Get script directory for path traversal -script_dir <- local({ - args <- commandArgs(trailingOnly = FALSE) - script <- grep("^--file=", args, value = TRUE) - dirname(normalizePath(sub("^--file=", "", script))) -}) - -# Import internal utilities -source(file.path(script_dir, "utils", "import.R")) -source(file.path(script_dir, "utils", "io.R")) -source(file.path(script_dir, "utils", "theme.R")) - -# Import data -df <- loadAnalyseCSV(analyse_path) - -# Define morphology features to plot -features <- c("equiv_diameter_nm", "major_axis_diameter", "minor_axis_diameter", "aspect_ratio", "eccentricity", "membrane_volume", "lumen_volume", "surface_area", "closure_fill_ratio") - -# For each feature, create a plot and spreadsheet showing distribution -stats <- list() -for (feat in intersect(features, names(df))) { - p <- ggplot2::ggplot(df, ggplot2::aes(x = .data[[feat]])) + - ggplot2::geom_density(fill = "#56B4E9", alpha = 0.4) + - ggplot2::geom_rug(alpha = 0.3) + - evTheme() + ggplot2::labs(title = feat) - saveSVG(p, file.path(output_dir, paste0(feat, "_density.svg"))) - v <- df[[feat]] - stats[[feat]] <- data.frame( - feature = feat, n = sum(!is.na(v)), mean = mean(v, na.rm = TRUE), - median = median(v, na.rm = TRUE), sd = sd(v, na.rm = TRUE), - p5 = quantile(v, 0.05, na.rm = TRUE), p95 = quantile(v, 0.95, na.rm = TRUE) - ) -} -writeXLSX(do.call(rbind, stats), file.path(output_dir, "summary_stats.xlsx")) diff --git a/src/evaluator/commands/plot/r/qc.R b/src/evaluator/commands/plot/r/qc.R deleted file mode 100644 index 0c3b3dc..0000000 --- a/src/evaluator/commands/plot/r/qc.R +++ /dev/null @@ -1,47 +0,0 @@ -# qc.R: plot qc metrics - -# Parse positional arguments -args <- commandArgs(trailingOnly = TRUE) -output_dir <- args[1]; analyse_path <- args[2]; mode <- args[3]; fill_ratio_flag_threshold <- as.numeric(args[4]) - -# Get script directory for path traversal -script_dir <- local({ - args <- commandArgs(trailingOnly = FALSE) - script <- grep("^--file=", args, value = TRUE) - dirname(normalizePath(sub("^--file=", "", script))) -}) - -# Import internal utilities -source(file.path(script_dir, "utils", "import.R")) -source(file.path(script_dir, "utils", "io.R")) -source(file.path(script_dir, "utils", "theme.R")) - -# Import data -df <- loadAnalyseCSV(analyse_path) - -# Calculate number of EVs per tomogram and plot -count_per_tomo <- dplyr::count(df, tomogram, name = "n_evs") -saveSVG(ggplot2::ggplot(count_per_tomo, ggplot2::aes(x = tomogram, y = n_evs)) + - ggplot2::geom_col(fill = "#0072B2") + evTheme() + - ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)), - file.path(output_dir, "count_per_tomogram.svg")) - -# Calculate number of enclosed EVs per tomogram and plot -pass_rate <- dplyr::summarise(dplyr::group_by(df, tomogram), - pass_rate = mean(is_enclosed, na.rm = TRUE), .groups = "drop") -saveSVG(ggplot2::ggplot(pass_rate, ggplot2::aes(x = tomogram, y = pass_rate)) + - ggplot2::geom_col(fill = "#009E73") + ggplot2::ylim(0, 1) + evTheme() + - ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)), - file.path(output_dir, "closure_pass_rate.svg")) - -saveSVG(ggplot2::ggplot(df, ggplot2::aes(x = closure_fill_ratio)) + - ggplot2::geom_histogram(bins = 30, fill = "#D55E00") + evTheme() + - ggplot2::geom_vline(xintercept = fill_ratio_flag_threshold, linetype = "dashed", colour = "black"), - file.path(output_dir, "closure_fill_ratio_distribution.svg")) - -writeXLSX(dplyr::summarise(dplyr::group_by(df, tomogram), - n_evs = dplyr::n(), pass_rate = mean(is_enclosed, na.rm = TRUE), - median_fill_ratio = median(closure_fill_ratio, na.rm = TRUE), - n_flagged_low_fill = sum(closure_fill_ratio < fill_ratio_flag_threshold, na.rm = TRUE), - .groups = "drop"), - file.path(output_dir, "qc_summary.xlsx")) diff --git a/src/evaluator/commands/plot/r/scatter.R b/src/evaluator/commands/plot/r/scatter.R deleted file mode 100644 index 519e14a..0000000 --- a/src/evaluator/commands/plot/r/scatter.R +++ /dev/null @@ -1,46 +0,0 @@ -# scatter.R: plot scatter graph of morphometrics - -# Parse positional arguments -args <- commandArgs(trailingOnly = TRUE) -output_dir <- args[1]; analyse_path <- args[2]; mode <- args[3] - -# Get script directory for path traversal -script_dir <- local({ - args <- commandArgs(trailingOnly = FALSE) - script <- grep("^--file=", args, value = TRUE) - dirname(normalizePath(sub("^--file=", "", script))) -}) - -# Import internal utilities -source(file.path(script_dir, "utils", "import.R")) -source(file.path(script_dir, "utils", "io.R")) -source(file.path(script_dir, "utils", "theme.R")) - -# Import data -df <- loadAnalyseCSV(analyse_path) - -# Define list of pairs to analyse -pairs <- list( - c("equiv_diameter_nm", "lumen_volume"), - c("surface_area", "lumen_volume"), - c("aspect_ratio", "equiv_diameter_nm") -) -stats_rows <- list() - -# For each pair, plot the scatter and write a spreadsheet with stats -for (pair in pairs) { - x <- pair[1]; y <- pair[2] - p <- ggplot2::ggplot(df, ggplot2::aes(x = .data[[x]], y = .data[[y]], colour = tomogram)) + - ggplot2::geom_point(alpha = 0.6) + - ggplot2::geom_smooth(method = "lm", se = FALSE, linewidth = 0.5) + - evTheme() + ggplot2::theme(legend.position = "none") - saveSVG(p, file.path(output_dir, paste0(x, "_vs_", y, ".svg"))) - fit <- lm(df[[y]] ~ df[[x]]) - stats_rows[[length(stats_rows) + 1]] <- data.frame( - x = x, y = y, - pearson_r = cor(df[[x]], df[[y]], use = "complete.obs"), - spearman_r = cor(df[[x]], df[[y]], use = "complete.obs", method = "spearman"), - slope = coef(fit)[2], slope_se = summary(fit)$coefficients[2, 2] - ) -} -writeXLSX(do.call(rbind, stats_rows), file.path(output_dir, "scatter_stats.xlsx")) diff --git a/src/evaluator/commands/plot/r/utils/import.R b/src/evaluator/commands/plot/r/utils/import.R deleted file mode 100644 index b89b472..0000000 --- a/src/evaluator/commands/plot/r/utils/import.R +++ /dev/null @@ -1,45 +0,0 @@ -# Data ingestion for plot sections - -loadAnalyseCSV <- function(path) { - df <- readr::read_csv(path, show_col_types = FALSE) - required <- c("tomogram", "label", "equiv_diameter_nm", "major_axis_diameter", "minor_axis_diameter", "aspect_ratio", "eccentricity", "membrane_volume", "lumen_volume", "surface_area", "is_enclosed", "closure_fill_ratio") - validateColumns(df, required) - df -} - -# model output may be JSON (default) or the flattened CSV form (nested fields JSON-encoded per-cell) -loadModelResults <- function(path) { - if (grepl("\\.json$", path)) { - payload <- jsonlite::fromJSON(path, simplifyDataFrame = TRUE, flatten = TRUE) - df <- payload$results - } else { - df <- readr::read_csv(path, show_col_types = FALSE) - for (col in intersect(c("reliability.is_reliable", "reliability"), names(df))) { - # nested reliability is JSON-encoded per cell in CSV mode; unpack is_reliable only - if (col == "reliability") { - parsed <- lapply(df$reliability, function(x) jsonlite::fromJSON(x)) - df$reliability.is_reliable <- vapply(parsed, function(x) isTRUE(x$is_reliable), logical(1)) - } - } - } - required <- c("label_id", "source_file", "chosen_model", "radius", "rmse_nm") - validateColumns(df, required) - df -} - -validateColumns <- function(df, required) { - missing <- setdiff(required, names(df)) - if (length(missing) > 0) { - stop(sprintf("Input missing required column(s): %s", paste(missing, collapse = ", "))) - } -} - -# Join analyse and model rows for the same vesicle -joinAnalyseModel <- function(analyse_df, model_df) { - model_df$tomogram <- basename(model_df$source_file) - dplyr::inner_join( - analyse_df, model_df, - by = c("tomogram" = "tomogram", "label" = "label_id"), - suffix = c(".analyse", ".model") - ) -} diff --git a/src/evaluator/commands/plot/r/utils/io.R b/src/evaluator/commands/plot/r/utils/io.R deleted file mode 100644 index 322b1cc..0000000 --- a/src/evaluator/commands/plot/r/utils/io.R +++ /dev/null @@ -1,40 +0,0 @@ -# Shared utilities for input/output operations - -# Save a plot to SVG file -saveSVG <- function(plot, path, width = 6, height = 4) { - ggplot2::ggsave(path, plot = plot, width = width, height = height, device = "svg") -} - -# Save a dataframe to a XLSX spreadsheet -writeXLSX <- function(df, path, sheet = "Sheet1") { - openxlsx::write.xlsx(df, path, sheetName = sheet, overwrite = TRUE) - stripDanglingDrawingRefs(path) -} - -# stripDanglingDrawingRefs: fix a bug in openxlsx (>= 4.2.8) where worksheet entries aren't included in zip file causing issue with some readers -stripDanglingDrawingRefs <- function(path) { - tmp <- tempfile() - dir.create(tmp) - on.exit(unlink(tmp, recursive = TRUE), add = TRUE) - zip::unzip(path, exdir = tmp) - - relsPath <- file.path(tmp, "xl", "worksheets", "_rels", "sheet1.xml.rels") - if (file.exists(relsPath)) { - s <- readLines(relsPath, warn = FALSE) - s <- gsub(']*Type="[^"]*(drawing|vmlDrawing)[^"]*"[^>]*/>', "", s) - writeLines(s, relsPath) - } - - ctPath <- file.path(tmp, "[Content_Types].xml") - if (file.exists(ctPath)) { - s <- readLines(ctPath, warn = FALSE) - s <- gsub(']*PartName="/xl/drawings/drawing1.xml"[^>]*/>', "", s) - writeLines(s, ctPath) - } - - files <- list.files(tmp, recursive = TRUE, all.files = TRUE) - unlink(path) - old <- setwd(tmp) - on.exit(setwd(old), add = TRUE) - zip::zip(path, files) -} \ No newline at end of file diff --git a/src/evaluator/commands/plot/r/utils/theme.R b/src/evaluator/commands/plot/r/utils/theme.R deleted file mode 100644 index 2a1135e..0000000 --- a/src/evaluator/commands/plot/r/utils/theme.R +++ /dev/null @@ -1,17 +0,0 @@ -# Shared theme for plotting - -evTheme <- function() { - ggplot2::theme_minimal(base_size = 11) + - ggplot2::theme( - panel.grid.minor = ggplot2::element_blank(), - legend.position = "right", - plot.title = ggplot2::element_text(face = "bold") - ) -} - -# Okabe-Ito for <=8 groups (colour-blind safe), Viridis beyond. -groupPalette <- function(groups) { - okabe_ito <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#000000") - n <- length(unique(groups)) - if (n <= length(okabe_ito)) okabe_ito[seq_len(n)] else viridisLite::viridis(n) -} diff --git a/src/evaluator/commands/plot/sections.py b/src/evaluator/commands/plot/sections.py new file mode 100644 index 0000000..07cc7c8 --- /dev/null +++ b/src/evaluator/commands/plot/sections.py @@ -0,0 +1,270 @@ +''' +======================================= +EValuator: PLOT SECTIONS +======================================= +''' + +# ==================== +# Import external dependencies +# ==================== +import numpy as np, pandas as pd, scipy.stats as stats +from dataclasses import dataclass +from pathlib import Path + +# ==================== +# Import internal EValuator utilities +# ==================== +from evaluator.utils.figures import plots +from evaluator.utils.settings import lg +from evaluator.commands.plot.utils import data as datautil +from evaluator.commands.plot.utils.export import write_figure, write_xlsx +from evaluator.commands.plot.utils.input import PlotRun + +# ==================== +# Define constant variables +# ==================== +# Morphology features reported by distributions and correlations +FEATURES = ['equiv_diameter_nm', 'xy_diameter_nm', datautil.SIZE_COL, 'major_axis_diameter', 'minor_axis_diameter', 'aspect_ratio', 'eccentricity', 'membrane_volume', 'lumen_volume', 'surface_area', 'closure_fill_ratio', 'sphere_rmse_rel',] +COMPARE_FEATURES = ['equiv_diameter_nm', 'aspect_ratio', 'closure_fill_ratio', 'lumen_volume'] +DIAMETERS = ['equiv_diameter_nm', 'major_axis_diameter'] +DEFAULT_AGREEMENT = [('equiv_diameter_nm', datautil.SIZE_COL)] + +# ==================== +# Define per-run section input +# ==================== +@dataclass(frozen=True) +class SectionData: + analyse: pd.DataFrame | None + model: pd.DataFrame | None + joined: pd.DataFrame | None + n_analyse_only: int = 0 + n_model_only: int = 0 + agreement: tuple[tuple[str, str], ...] = () + analyse_all: pd.DataFrame | None = None + +def load(run: PlotRun, *, include_flagged: bool, agreement=None) -> SectionData: + pairs = tuple(agreement or ()) + analyse = datautil.load_analyse(run.analyse_path, include_flagged=include_flagged) if run.analyse_path else None + analyse_all = datautil.load_analyse(run.analyse_path, include_flagged=True) if run.analyse_path else None # unfiltered, for report-both-ways stats + model = datautil.load_model(run.model_path) if run.model_path else None + if analyse is None or model is None: + return SectionData(analyse, model, None, agreement=pairs, analyse_all=analyse_all) + joined, n_analyse_only, n_model_only = datautil.join_analyse_model(analyse, model) # before the attach below, so the joined frame keeps one model_diameter_nm + analyse = datautil.attach_fitted_diameter(analyse, model) + analyse_all = datautil.attach_fitted_diameter(analyse_all, model) + return SectionData(analyse, model, joined, n_analyse_only, n_model_only, pairs, analyse_all) + +# ==================== +# Define sections +# ==================== +def distributions(out: Path, data: SectionData, params) -> None: + df = data.analyse + rows = [] + for feat in FEATURES: + if feat not in df.columns: + continue + width = params.bin_width_overrides.get(feat, params.bin_width) + write_figure(plots.distribution(df, feat, set(), bin_size=width or None), out / f'{feat}_histogram.svg') + for population, subset in datautil.populations(data.analyse_all).items(): + v = pd.to_numeric(subset[feat], errors='coerce') + rows.append({ + 'population': population, 'feature': feat, 'n': int(v.notna().sum()), 'mean': v.mean(), 'median': v.median(), + 'sd': v.std(), 'p5': v.quantile(0.05), 'p95': v.quantile(0.95), + }) + notes = pd.DataFrame({'note': [ + 'population: all = every component; vesicle_like = passes the QC filters; vesicle_like_no_z_clip = also excludes components touching a z face', + 'a conclusion that holds in all populations does not depend on the QC filter or on z clipping', + 'xy_diameter_nm is independent of z clipping; aspect_ratio, eccentricity and sphere_rmse_rel are biased low for clipped components, preferring vesicle_like_no_z_clip may be useful with caveat that excluding clipped EVs biases size small', + ]}) + write_xlsx(out / 'summary_stats.xlsx', {'summary': pd.DataFrame(rows), 'notes': notes}) + +def scatter(out: Path, data: SectionData, params) -> None: + df = data.analyse + _log_pooling(df) + pairs = [(_size(df, x), _size(df, y)) for x, y in params.scatter_pairs] + tomo = datautil.aggregate_tomograms(df, list(dict.fromkeys(c for pair in pairs for c in pair if c in df.columns)), keys=('tomogram',)) + rows, tomo_rows = [], [] + for x, y in pairs: + if x not in df.columns or y not in df.columns: + lg.warning(f'plot | scatter: skipping {x!r} vs {y!r}, column missing') + continue + write_figure(plots.feature_scatter(df, x, y, set()), out / f'{x}_vs_{y}.svg') + rows.append(_pair_stats(df, x, y)) + tomo_rows.append(_pair_stats(tomo, x, y)) + notes = pd.concat([_notes(df), pd.DataFrame({'note': [f'stats_tomogram: regression and correlation on per-tomogram medians (n = {len(tomo)} tomograms), the cluster-aware estimate; stats is pooled over EVs and descriptive only']})], ignore_index=True) + write_xlsx(out / 'scatter_stats.xlsx', {'stats': pd.DataFrame(rows), 'stats_tomogram': pd.DataFrame(tomo_rows), 'notes': notes}) + +def correlations(out: Path, data: SectionData, params) -> None: + df = data.analyse + cols = [c for c in FEATURES if c in df.columns and pd.to_numeric(df[c], errors='coerce').notna().any()] + if len(cols) < 2: + raise ValueError('correlations needs at least two numeric feature columns') + _log_pooling(df) + tomo = datautil.aggregate_tomograms(df, cols, keys=('tomogram',)) + sheets = {} + for method in ('spearman', 'pearson'): + write_figure(plots.correlation_matrix(df, cols, method=method), out / f'correlation_{method}.svg', 'heatmap') + sheets[method] = df[cols].apply(pd.to_numeric, errors='coerce').corr(method=method).reset_index(names='feature') + sheets[f'{method}_tomogram'] = tomo[cols].corr(method=method).reset_index(names='feature') + sheets['notes'] = pd.concat([_notes(df), pd.DataFrame({'note': [f'*_tomogram sheets: correlations of per-tomogram medians (n = {len(tomo)} tomograms), the cluster-aware estimate']})], ignore_index=True) + write_xlsx(out / 'correlations.xlsx', sheets) + +def qc(out: Path, data: SectionData, params) -> None: + df = data.analyse + thr = params.fill_ratio_flag_threshold + enclosed = df['is_enclosed'].astype(float) + low = (df['closure_fill_ratio'] < thr).groupby(df['tomogram']).sum().rename('n_flagged_low_fill') + per = ( + df.assign(_enclosed=enclosed).groupby('tomogram') + .agg(n_evs=('label', 'size'), pass_rate=('_enclosed', 'mean'), median_fill_ratio=('closure_fill_ratio', 'median')) + .join(low).reset_index() + ) + write_figure(plots.bar(per, 'tomogram', 'n_evs', title='EV count per tomogram', y_title='EVs'), out / 'count_per_tomogram.svg') + write_figure(plots.bar(per, 'tomogram', 'pass_rate', title='Closure pass rate per tomogram', y_title='Fraction enclosed', y_range=(0, 1)), out / 'closure_pass_rate.svg') + fig = plots.distribution(df, 'closure_fill_ratio', set(), bin_size=0.05) + fig.add_vline(x=thr, line_dash='dash', line_color=plots.ACTIVE['font']) + write_figure(fig, out / 'closure_fill_ratio_distribution.svg') + write_xlsx(out / 'qc_summary.xlsx', {'summary': per}) + +def concordance(out: Path, data: SectionData, params) -> None: + j = data.joined + if j is None or j.empty: + raise ValueError('no matching (tomogram, label) rows between analyse and model output, check both were run on the same labelled MRC(s)') + rows = [] + for col in DIAMETERS: + fig = plots.concordance(j, set(), analyse_col=col) + if fig is None: + continue + write_figure(fig, out / f'model_vs_{col}.svg', 'square') + a, b = _paired(j[col], j[datautil.SIZE_COL]) + diff = b - a + rows.append({ + 'analyse_measure': col, 'n': len(a), 'bias': diff.mean(), 'rmse': float(np.sqrt((diff ** 2).mean())), + 'pearson_r': a.corr(b), + }) + match = pd.DataFrame([{'n_matched': len(j), 'n_analyse_only': data.n_analyse_only, 'n_model_only': data.n_model_only}]) + write_xlsx(out / 'concordance_stats.xlsx', {'agreement': pd.DataFrame(rows), 'matching': match}) + +def reliability(out: Path, data: SectionData, params) -> None: + j = data.joined + fig = plots.reliability(j, set()) if j is not None else None + if fig is None: + raise ValueError('reliability needs model rmse and analyse closure columns') + write_figure(fig, out / 'reliability.svg') + row = {'n': len(j)} + rel_col = plots.find_col(j, 'is_reliable') + if rel_col is not None: + rel = j[rel_col].fillna(False).astype(bool) + row |= {'n_reliable': int(rel.sum()), 'reliable_pass_rate': float(rel.mean())} + if 'is_vesicle_like' in j.columns: + row['n_vesicle_like'] = int(j['is_vesicle_like'].fillna(False).astype(bool).sum()) + write_xlsx(out / 'reliability_summary.xlsx', {'summary': pd.DataFrame([row])}) + +def agreement(out: Path, data: SectionData, params) -> None: + j = data.joined + if j is None or j.empty: + raise ValueError('agreement needs matching analyse and model rows') + rows = [] + for a, b in (data.agreement or DEFAULT_AGREEMENT): + fig = plots.bland_altman(j, a, b, set()) + if fig is None: + lg.warning(f'plot | agreement: skipping {a!r} vs {b!r}, column missing or fewer than two paired values') + continue + write_figure(fig, out / f'{a}_vs_{b}.svg') + x, y = _paired(j[a], j[b]) + diff, sd = x - y, (x - y).std(ddof=1) + rows.append({ + 'a': a, 'b': b, 'n': len(x), 'bias': diff.mean(), + 'loa_lower': diff.mean() - 1.96 * sd, 'loa_upper': diff.mean() + 1.96 * sd, + }) + if not rows: + raise ValueError('agreement: no plottable measurement pair') + write_xlsx(out / 'agreement_stats.xlsx', {'stats': pd.DataFrame(rows)}) + +def compare(out: Path, runs: list[PlotRun], params, *, include_flagged: bool) -> None: + frames = [] + for r in runs: + if not (r.analyse_path and r.group): + continue + df = datautil.load_analyse(r.analyse_path, include_flagged=include_flagged) + if r.model_path: + df = datautil.attach_fitted_diameter(df, datautil.load_model(r.model_path)) + frames.append(df.assign(sample_id=r.sample_id, group=r.group, replicate=r.replicate)) + combined = pd.concat(frames, ignore_index=True) + feats = [_size(combined, f) for f in COMPARE_FEATURES if f in combined.columns] + _log_pooling(combined, n_samples=combined['sample_id'].nunique()) + tomo = datautil.aggregate_tomograms(combined, feats) + lg.info(f'plot | caveat: {datautil.tomogram_note(tomo)}') + for feat in feats: + write_figure(plots.group_distribution(tomo, feat), out / f'violin_{feat}.svg') + fig = plots.pca_scatter(tomo, feats) + if fig is not None: + write_figure(fig, out / 'pca.svg', 'square') + summary = tomo.groupby('group')[feats].agg(['mean', 'median']) + summary.columns = [f'{f}_{s}' for f, s in summary.columns] + counts = tomo.groupby('group').agg(n_tomograms=('tomogram', 'size'), n_samples=('sample_id', 'nunique'), n_evs=('n_evs', 'sum')) + replicates = combined.groupby('group')['replicate'].nunique() + if replicates.any(): # only when the sample sheet carries a replicate column + counts['n_replicates'] = replicates + notes = pd.DataFrame({'note': [datautil.tomogram_note(tomo), 'summary columns are the mean and median across per-tomogram medians']}) + write_xlsx(out / 'group_summary.xlsx', {'summary': counts.join(summary).reset_index(), 'notes': notes}) + write_xlsx(out / 'group_tests.xlsx', {'tests': _group_tests(tomo, feats), 'notes': notes}) + +# SECTIONS: per-run section functions, keyed by section name (compare pools every run so is called separately) +SECTIONS = { + 'distributions': distributions, + 'scatter': scatter, + 'correlations': correlations, + 'qc': qc, + 'concordance': concordance, + 'reliability': reliability, + 'agreement': agreement, +} + +# ==================== +# Define private helpers +# ==================== +# _group_tests: tomogram-level test per feature (Mann-Whitney U for two groups, Kruskal-Wallis for more) +def _group_tests(tomo: pd.DataFrame, feats: list[str]) -> pd.DataFrame: + rows = [] + for feat in feats: + samples = {g: pd.to_numeric(s[feat], errors='coerce').dropna() for g, s in tomo.groupby('group')} + samples = {g: v for g, v in samples.items() if len(v)} + if len(samples) < 2: + continue + name, test = ('mann-whitney', stats.mannwhitneyu) if len(samples) == 2 else ('kruskal-wallis', stats.kruskal) + try: + res = test(*samples.values()) + statistic, p_value = res.statistic, res.pvalue + except ValueError: # e.g. every value identical + statistic = p_value = np.nan + rows.append({'feature': feat, 'test': name, 'statistic': statistic, 'p_value': p_value, **{f'n_tomograms_{g}': len(v) for g, v in samples.items()}}) + return pd.DataFrame(rows) + +# _size: swap the shell-volume equivalent diameter for the fitted diameter when model output is available +def _size(df: pd.DataFrame, col: str) -> str: + if col == 'equiv_diameter_nm' and datautil.SIZE_COL in df.columns and df[datautil.SIZE_COL].notna().any(): + return datautil.SIZE_COL + return col + +def _paired(a, b) -> tuple[pd.Series, pd.Series]: + a, b = pd.to_numeric(a, errors='coerce'), pd.to_numeric(b, errors='coerce') + m = a.notna() & b.notna() + return a[m], b[m] + +def _pair_stats(df: pd.DataFrame, x: str, y: str) -> dict: + xs, ys = _paired(df[x], df[y]) + row = {'x': x, 'y': y, 'n': len(xs), 'pearson_r': np.nan, 'spearman_r': np.nan, 'slope': np.nan, 'slope_se': np.nan, 'r2': np.nan} + if len(xs) >= 3 and xs.nunique() > 1: + fit = stats.linregress(xs, ys) + row |= { + 'pearson_r': fit.rvalue, 'spearman_r': stats.spearmanr(xs, ys).statistic, + 'slope': fit.slope, 'slope_se': fit.stderr, 'r2': fit.rvalue ** 2, + } + return row + +def _notes(df: pd.DataFrame, n_samples: int = 1) -> pd.DataFrame: + return pd.DataFrame({'note': [datautil.pooling_note(df, n_samples)]}) + +def _log_pooling(df: pd.DataFrame, n_samples: int = 1) -> None: + lg.info(f'plot | caveat: {datautil.pooling_note(df, n_samples)}') diff --git a/src/evaluator/commands/plot/utils/data.py b/src/evaluator/commands/plot/utils/data.py new file mode 100644 index 0000000..4fe4a0a --- /dev/null +++ b/src/evaluator/commands/plot/utils/data.py @@ -0,0 +1,111 @@ +''' +======================================= +EValuator: PLOT DATA PREPROCESSING +======================================= +''' + +# ==================== +# Import external dependencies +# ==================== +import pandas as pd +from pathlib import Path + +# ==================== +# Import internal EValuator utilities +# ==================== +from evaluator.utils import io as ioutil +from evaluator.utils.settings import lg + +# ==================== +# Define constant variables +# ==================== +ANALYSE_REQUIRED = ['tomogram', 'label', 'equiv_diameter_nm', 'major_axis_diameter', 'minor_axis_diameter', 'aspect_ratio', 'eccentricity', 'membrane_volume', 'lumen_volume', 'surface_area', 'is_enclosed', 'closure_fill_ratio'] +MODEL_REQUIRED = ['label_id', 'source_file', 'chosen_model', 'radius', 'rmse_nm'] +# SIZE_COL: fitted-sphere diameter from `model`, used as the size axis whenever model output is supplied +SIZE_COL = 'model_diameter_nm' + +# ==================== +# Define loaders +# ==================== +def load_analyse(path: Path, *, include_flagged: bool = False) -> pd.DataFrame: + df = pd.read_csv(path) + _validate(df, ANALYSE_REQUIRED, path.name) + if include_flagged or 'is_vesicle_like' not in df.columns: # older csvs predate the column: nothing to filter + return df + keep = df['is_vesicle_like'].fillna(False).astype(bool) + dropped = df[~keep] + if len(dropped): + lg.info(f'plot | {path.name}: excluded {len(dropped)}/{len(df)} component(s) flagged as not vesicle-like (is_vesicle_like == False)') + _log_drop_reasons(dropped) + return df[keep].reset_index(drop=True) + +def load_model(path: Path) -> pd.DataFrame: + records, _ = ioutil.read_results(path) + df = pd.json_normalize(records) # flattens nested reliability / fit dicts to 'reliability.is_reliable' etc. + _validate(df, MODEL_REQUIRED, path.name) + df[SIZE_COL] = pd.to_numeric(df['radius'], errors='coerce') * 2 + return df + +# join_analyse_model: inner join on (tomogram basename, label); also returns the unmatched counts on each side +def join_analyse_model(analyse: pd.DataFrame, model: pd.DataFrame) -> tuple[pd.DataFrame, int, int]: + model = model.assign(tomogram=model['source_file'].map(lambda s: Path(s).name)) + joined = analyse.merge(model, left_on=['tomogram', 'label'], right_on=['tomogram', 'label_id'], how='inner', suffixes=('', '_model')) + return joined, len(analyse) - len(joined), len(model) - len(joined) + +# aggregate_tomograms: one row per tomogram holding the median of each feature and the EV count, so tomograms (not EVs) are the unit of analysis +def aggregate_tomograms(df: pd.DataFrame, feats: list[str], keys: tuple[str, ...] = ('group', 'sample_id', 'tomogram')) -> pd.DataFrame: + keys = [k for k in keys if k in df.columns] + num = df[keys].join(df[feats].apply(pd.to_numeric, errors='coerce')) + grouped = num.groupby(keys, dropna=False) + out = grouped[feats].median() + out.insert(0, 'n_evs', grouped.size()) + return out.reset_index() + +# tomogram_note: caveat text for statistics computed on per-tomogram medians +def tomogram_note(tomo: pd.DataFrame) -> str: + n_samples = tomo['sample_id'].nunique() if 'sample_id' in tomo.columns else 1 + return (f'n = {len(tomo)} tomogram(s) from {n_samples} sample(s); each tomogram contributes its median, so EVs are not treated as independent. ') + +# attach_fitted_diameter: left-join model's fitted diameter onto analyse, masking unreliable fits, so figures can use it as the size axis +def attach_fitted_diameter(analyse: pd.DataFrame, model: pd.DataFrame) -> pd.DataFrame: + fitted = pd.DataFrame({ + 'tomogram': model['source_file'].map(lambda s: Path(s).name), + 'label': model['label_id'], + SIZE_COL: model[SIZE_COL], + }) + if 'reliability.is_reliable' in model.columns: + fitted[SIZE_COL] = fitted[SIZE_COL].where(model['reliability.is_reliable'].fillna(False).astype(bool)) + return analyse.merge(fitted.drop_duplicates(['tomogram', 'label']), on=['tomogram', 'label'], how='left') + +# populations: named subsets of the unfiltered analyse frame, for reporting the QC filter and z clipping both ways +def populations(df: pd.DataFrame) -> dict[str, pd.DataFrame]: + out = {'all': df} + if 'is_vesicle_like' not in df.columns: + return out + keep = df['is_vesicle_like'].fillna(False).astype(bool) + out['vesicle_like'] = df[keep] + if 'touches_z_boundary' in df.columns: + clipped = df['touches_z_boundary'].fillna(False).astype(bool) + out['vesicle_like_no_z_clip'] = df[keep & ~clipped] + return out + +# pooling_note: caveat text for any statistic that treats pooled EVs as independent observations +def pooling_note(df: pd.DataFrame, n_samples: int = 1) -> str: + return (f'n = {len(df)} EVs pooled across {df["tomogram"].nunique()} tomogram(s) and {n_samples} sample(s); EVs are treated as independent and no clustering-adjusted test is applied, so treat statistics as descriptive only') + +# ==================== +# Define private helpers +# ==================== +def _validate(df: pd.DataFrame, required: list[str], name: str) -> None: + missing = [c for c in required if c not in df.columns] + if missing: + raise ValueError(f'{name}: input missing required column(s): {", ".join(missing)}') + +def _log_drop_reasons(dropped: pd.DataFrame) -> None: + if 'qc_flags' not in dropped.columns: + return + # qc_flags is comma-separated, so a component with two flags counts once under each + flags = dropped['qc_flags'].dropna().astype(str).str.split(',').explode().str.strip() + counts = flags[flags != ''].value_counts() + if not counts.empty: + lg.info('plot | by reason: ' + ', '.join(f'{reason}={n}' for reason, n in counts.items())) diff --git a/src/evaluator/commands/plot/utils/dispatch.py b/src/evaluator/commands/plot/utils/dispatch.py deleted file mode 100644 index 13a71d2..0000000 --- a/src/evaluator/commands/plot/utils/dispatch.py +++ /dev/null @@ -1,57 +0,0 @@ -''' -======================================= -EValuator: PLOT COMMAND R SCRIPT DISPATCH -======================================= -''' - -# ==================== -# Import external dependencies -# ==================== -import shutil, subprocess -from importlib.resources import files as pkg_files -from pathlib import Path - -# ==================== -# Import internal EValuator dependencies -# ==================== -from evaluator.utils.settings import lg - -# ==================== -# Define custom error classes -# ==================== -class RscriptNotFoundError(RuntimeError): - '''Raised when no usable Rscript binary can be resolved''' - -class RscriptError(RuntimeError): - '''Raised when an R script exits non-zero''' - -# ==================== -# Define Rscript dispatch functions -# ==================== -def resolve_rscript(rscript: Path | None) -> Path: - if rscript is not None: - if not rscript.exists(): - raise RscriptNotFoundError(f'Rscript not found at {rscript}') - return rscript - found = shutil.which("Rscript") - if found is None: - raise RscriptNotFoundError('Rscript not found on PATH. Install R to PATH, or pass --rscript.') - return Path(found) - -def _script_path(name: str) -> Path: - return Path(str(pkg_files('evaluator.commands.plot') / 'r' / name)) - -def dispatch(rscript_bin: Path, script_name: str, args: list[str]) -> None: - ''' - Run r/ with positional args, streaming stderr as it's produced - ''' - script = _script_path(script_name) - cmd = [str(rscript_bin), str(script), *[str(a) for a in args]] - lg.debug("plot | Running: {}", ' '.join(cmd)) - proc = subprocess.run(cmd, capture_output=True, text=True) - if proc.stdout: - lg.debug("plot | {} stdout:\n{}", script_name, proc.stdout) - if proc.returncode != 0: - raise RscriptError(f"{script_name} failed (exit {proc.returncode}):\n{proc.stderr}") - if proc.stderr: - lg.debug("plot | {} stderr:\n{}", script_name, proc.stderr) diff --git a/src/evaluator/commands/plot/utils/export.py b/src/evaluator/commands/plot/utils/export.py new file mode 100644 index 0000000..7a88df4 --- /dev/null +++ b/src/evaluator/commands/plot/utils/export.py @@ -0,0 +1,58 @@ +''' +======================================= +EValuator: PLOT AND TABLE EXPORT +======================================= +''' + +# ==================== +# Import external dependencies +# ==================== +import pandas as pd, plotly.graph_objects as go +from pathlib import Path + +# ==================== +# Define constants +# ==================== +# Export size (px) per figure kind +SIZES: dict[str, tuple[int, int]] = { + 'default': (700, 480), + 'square': (560, 560), + 'wide': (900, 480), + 'heatmap': (720, 640), +} + +# ==================== +# Define export functions +# ==================== +# check_kaleido: check if static rendering is available +def check_kaleido() -> None: + try: + go.Figure().to_image(format='svg') + except Exception as e: + raise RuntimeError(f'plot | static export needs the plot extra (kaleido) and a Chrome/Chromium install (try `plotly_get_chrome`): {e}') from e + +def write_figure(fig: go.Figure, path: Path, kind: str = 'default') -> None: + width, height = SIZES[kind] + _vectorise(fig).write_image(path, format='svg', width=width, height=height) + +def write_xlsx(path: Path, sheets: dict[str, pd.DataFrame]) -> None: + with pd.ExcelWriter(path, engine='openpyxl') as xw: + for name, df in sheets.items(): # first sheet is the active one, so put the main table first + df.to_excel(xw, sheet_name=name[:31], index=False) + +# ==================== +# Define private helpers +# ==================== +# _vectorise: scattergl traces are rasterised in SVG output, so swap them for plain scatter traces +def _vectorise(fig: go.Figure) -> go.Figure: + out = go.Figure(fig) + traces = [] + for tr in out.data: + spec = tr.to_plotly_json() + if spec.get('type') == 'scattergl': + spec.pop('type') + tr = go.Scatter(**spec) + traces.append(tr) + out.data = () + out.add_traces(traces) + return out diff --git a/src/evaluator/commands/plot/utils/input.py b/src/evaluator/commands/plot/utils/input.py index 5f0494a..76e5c4f 100644 --- a/src/evaluator/commands/plot/utils/input.py +++ b/src/evaluator/commands/plot/utils/input.py @@ -81,9 +81,19 @@ def available_sections(runs: list[PlotRun], multi_run: bool) -> list[str]: has_model = any(r.model_path for r in runs) sections = [] if has_analyse: - sections += ['distributions', 'qc', 'scatter'] + sections += ['distributions', 'qc', 'scatter', 'correlations'] if has_analyse and has_model: - sections.append("concordance") + sections += ['concordance', 'reliability', 'agreement'] if multi_run and len({r.group for r in runs if r.group} ) > 1: sections.append('compare') return sections + +# parse_agreement: split --agreement values into (A, B) measurement pairs; each value may hold several space-separated pairs +def parse_agreement(values: list[str]) -> list[tuple[str, str]]: + pairs = [] + for token in (t for v in values for t in v.split()): + parts = token.split(',') + if len(parts) != 2 or not all(parts): + raise ValueError(f'invalid --agreement pair {token!r}, expected A,B') + pairs.append((parts[0], parts[1])) + return pairs diff --git a/src/evaluator/commands/tools/benchmark_geometric_proxy/benchmark_geometric_proxy.py b/src/evaluator/commands/tools/benchmark_geometric_proxy/benchmark_geometric_proxy.py index ef43667..b87a0f3 100644 --- a/src/evaluator/commands/tools/benchmark_geometric_proxy/benchmark_geometric_proxy.py +++ b/src/evaluator/commands/tools/benchmark_geometric_proxy/benchmark_geometric_proxy.py @@ -63,6 +63,9 @@ def calculate_benchmark( t0 = time.perf_counter() new_centre, new_radius = estimateCentroidRadius(points) new_time = time.perf_counter() - t0 + t0 = time.perf_counter() + diag_centre, diag_radius = methods.naive_diagonal_radius(points) + diag_time = time.perf_counter() - t0 rows.append({ 'cloud_type': cloud_type, 'completeness_deg': completeness_deg, @@ -74,6 +77,9 @@ def calculate_benchmark( 'new_radius_nm': new_radius, 'new_centre_error_nm': float(np.linalg.norm(new_centre - true_centre)), 'new_time_s': new_time, + 'diag_radius_nm': diag_radius, + 'diag_centre_error_nm': float(np.linalg.norm(diag_centre - true_centre)), + 'diag_time_s': diag_time, }) point_index += 1 console.progress_bar(point_index, total_points) diff --git a/src/evaluator/commands/tools/benchmark_geometric_proxy/defaults.py b/src/evaluator/commands/tools/benchmark_geometric_proxy/defaults.py index 3c8c736..358879c 100644 --- a/src/evaluator/commands/tools/benchmark_geometric_proxy/defaults.py +++ b/src/evaluator/commands/tools/benchmark_geometric_proxy/defaults.py @@ -20,4 +20,5 @@ METHODS = { 'naive bounding box': ('old_radius_nm', 'old_time_s'), 'isotropy-aware bounding box': ('new_radius_nm', 'new_time_s'), + 'diagonal/2 baseline': ('diag_radius_nm', 'diag_time_s'), } diff --git a/src/evaluator/commands/tools/benchmark_geometric_proxy/methods.py b/src/evaluator/commands/tools/benchmark_geometric_proxy/methods.py index 6bb8e5d..73df4bc 100644 --- a/src/evaluator/commands/tools/benchmark_geometric_proxy/methods.py +++ b/src/evaluator/commands/tools/benchmark_geometric_proxy/methods.py @@ -21,3 +21,13 @@ def bbox_centroid_radius(points: np.ndarray) -> tuple[np.ndarray, float]: extent = points.max(axis=0) - points.min(axis=0) radius_estimate = float(np.mean(extent) / 2.0) return centroid, radius_estimate + +# -- naive_diagonal_radius: returns tuple of ndarray, float +def naive_diagonal_radius(points: np.ndarray) -> tuple[np.ndarray, float]: + ''' + Estimate centre/radius as point mean/half the bounding-box diagonal (independent naive baseline) + ''' + centroid = points.mean(axis=0) + extent = points.max(axis=0) - points.min(axis=0) + radius_estimate = float(np.linalg.norm(extent) / 2.0) + return centroid, radius_estimate diff --git a/src/evaluator/commands/tools/benchmark_geometric_proxy/summaries.py b/src/evaluator/commands/tools/benchmark_geometric_proxy/summaries.py index 482a065..42ce6c9 100644 --- a/src/evaluator/commands/tools/benchmark_geometric_proxy/summaries.py +++ b/src/evaluator/commands/tools/benchmark_geometric_proxy/summaries.py @@ -18,7 +18,7 @@ # Define constants # ==================== CLOUD_TYPE_LABELS = {'polar_cap': 'Polar cap', 'equatorial_band': 'Equatorial band'} -METHOD_COLOURS = {'naive bounding box': 'tab:red', 'isotropy-aware bounding box': 'tab:blue'} +METHOD_COLOURS = {'naive bounding box': 'tab:red', 'isotropy-aware bounding box': 'tab:blue', 'diagonal/2 baseline': 'tab:green'} # ==================== # Define summary functions diff --git a/src/evaluator/commands/tools/benchmark_missing_wedge/benchmark_missing_wedge.py b/src/evaluator/commands/tools/benchmark_missing_wedge/benchmark_missing_wedge.py index aa7c8df..4df5f77 100644 --- a/src/evaluator/commands/tools/benchmark_missing_wedge/benchmark_missing_wedge.py +++ b/src/evaluator/commands/tools/benchmark_missing_wedge/benchmark_missing_wedge.py @@ -50,14 +50,14 @@ def calculate_benchmark( offset = tuple(rng.uniform(-0.5, 0.5, size=3)) axes = tuple(rng.normal(1.0, shape_jitter, size=3)) # Generate synthetic EVs across specified diameter range - shell, truth = synthetic.generate_ev_shell(diameter_nm=d, voxel_size_nm=voxel_size_nm, centre_offset_voxels=offset, axis_ratios=axes) + shell, truth = synthetic.generate_ev_shell(diameter_nm=d, thickness_nm=defaults.THICKNESS, voxel_size_nm=voxel_size_nm, centre_offset_voxels=offset, axis_ratios=axes) true_diameter = truth['diameter_nm'] # Apply fourier missing wedge degradation degraded = synthetic.apply_fourier_missing_wedge(shell, tilt_range_deg=tilt_range_deg) # Baseline (no mitigation) t0 = time.perf_counter() - baseline_d = methods._shell_voxel_diameter(degraded, voxel_size_nm) + baseline_d = methods._shell_voxel_diameter(degraded, voxel_size_nm, defaults.THICKNESS) baseline_time = time.perf_counter() - t0 # Apply mitigation 1: anisotropic closing → shell-volume diameter @@ -79,7 +79,7 @@ def calculate_benchmark( # Apply mitigation 3: convex hull outer diameter t0 = time.perf_counter() - hull_d = methods.hull_outer_diameter(shell, voxel_size_nm) + hull_d = methods.hull_outer_diameter(degraded, voxel_size_nm) hull_time = time.perf_counter() - t0 # Apply mitigation 4: XY/Z diagnostic diff --git a/src/evaluator/commands/tools/benchmark_missing_wedge/defaults.py b/src/evaluator/commands/tools/benchmark_missing_wedge/defaults.py index f913cc1..271d206 100644 --- a/src/evaluator/commands/tools/benchmark_missing_wedge/defaults.py +++ b/src/evaluator/commands/tools/benchmark_missing_wedge/defaults.py @@ -9,6 +9,7 @@ DIAMETER_STEP = 40 REPLICATES = 10 VOXEL_SIZE = 2.0 +THICKNESS = 5.0 # matches synthetic.generate_ev_shell's default membrane thickness TILT_RANGE = 60.0 DIAMETER_JITTER = 6.0 SHAPE_JITTER = 0.1 diff --git a/src/evaluator/commands/tools/benchmark_missing_wedge/methods.py b/src/evaluator/commands/tools/benchmark_missing_wedge/methods.py index 3cbef6f..b6449e4 100644 --- a/src/evaluator/commands/tools/benchmark_missing_wedge/methods.py +++ b/src/evaluator/commands/tools/benchmark_missing_wedge/methods.py @@ -15,15 +15,27 @@ # ==================== # Define method functions # ==================== -# -- shell_voxel_diameter: returns float corresponding to raw shell voxel count -def _shell_voxel_diameter(binary: np.ndarray, voxel_size_nm: float) -> float: +# -- _shell_voxel_diameter: returns float corresponding to raw shell voxel count +def _shell_voxel_diameter(binary: np.ndarray, voxel_size_nm: float, thickness_nm: float) -> float: ''' - Calculate equivalent diameter from raw shell voxel count + Calculate outer diameter from raw shell voxel count, given a known membrane thickness. ''' n = int(binary.sum()) + volume_nm3 = n * voxel_size_nm ** 3 + t = thickness_nm + # V = (4/3)pi*(3*r^2*t - 3*r*t^2 + t^3) => 3t*r^2 - 3t^2*r + (t^3 - 3V/(4pi)) = 0 + a = 3 * t + b = -3 * t ** 2 + c = t ** 3 - (3 * volume_nm3) / (4 * np.pi) + discriminant = b ** 2 - 4 * a * c + if discriminant < 0: + # shell is thinner/sparser than a full continuous shell of thickness t would allow at this + # voxel count; fall back to the solid-sphere-equivalent as a conservative estimate + return 2 * (3 * volume_nm3 / (4 * np.pi)) ** (1 / 3) + r_out = (-b + np.sqrt(discriminant)) / (2 * a) + return 2 * r_out if n == 0: return float('nan') - return 2 * (3 * n * voxel_size_nm ** 3 / (4 * np.pi)) ** (1 / 3) # -- anisotropic_closing_per_component: returns ndarray with closed binary segmentation def anisotropic_closing_per_component( diff --git a/src/evaluator/commands/tools/utils/synthetic.py b/src/evaluator/commands/tools/utils/synthetic.py index fc2bc5a..8928e5c 100644 --- a/src/evaluator/commands/tools/utils/synthetic.py +++ b/src/evaluator/commands/tools/utils/synthetic.py @@ -34,9 +34,10 @@ def generate_ev_shell( dx = (x - centre[2]) / axis_ratios[2] dist = np.sqrt(dz ** 2 + dy ** 2 + dx ** 2) shell = np.abs(dist - radius_voxels) <= half_thickness_voxels - + volume_equivalent_diameter_nm = diameter_nm * (float(np.prod(axis_ratios))) ** (1 / 3) truth = { - "diameter_nm": diameter_nm, + "diameter_nm": volume_equivalent_diameter_nm, + "nominal_diameter_nm": diameter_nm, "thickness_nm": thickness_nm, "centre_voxels": centre, "voxel_size_nm": voxel_size_nm, diff --git a/src/evaluator/commands/viewer/app/pages/tomogram.py b/src/evaluator/commands/viewer/app/pages/tomogram.py index 898b2df..a9a2a22 100644 --- a/src/evaluator/commands/viewer/app/pages/tomogram.py +++ b/src/evaluator/commands/viewer/app/pages/tomogram.py @@ -15,13 +15,13 @@ # Import EValuator utilities from evaluator.utils import io as ioutil from evaluator.utils import mrc as mrcutil +from evaluator.utils.figures import plots as plotutil +from evaluator.utils.figures.format import pretty_column # Import EValuator viewer utilities from evaluator.commands.viewer.app.components.plotly_camera_sync import plotly_view -from evaluator.commands.viewer.utils import plots as plotutil from evaluator.commands.viewer.utils import theme as themeutil from evaluator.commands.viewer.utils.export import export_filtered_csv -from evaluator.commands.viewer.utils.format import pretty_column from evaluator.commands.viewer.utils.join import join_analyse_model from evaluator.commands.viewer.utils.mesh import build_label_mesh_traces, build_point_cloud_trace, dim_trace diff --git a/src/evaluator/commands/viewer/utils/theme.py b/src/evaluator/commands/viewer/utils/theme.py index 48edb26..3fbfda5 100644 --- a/src/evaluator/commands/viewer/utils/theme.py +++ b/src/evaluator/commands/viewer/utils/theme.py @@ -10,39 +10,13 @@ import streamlit as st # ==================== -# Define constants +# Import shared theme constants # ==================== -# Named themes -THEMES: dict[str, dict] = { - 'Okabe-Ito': { - 'palette': ['#E69F00', '#56B4E9', '#009E73', '#F0E442', '#0072B2', '#D55E00', '#CC79A7', '#999999'], - 'highlight': '#FFD400', 'base': '#56B4E9', 'reliable': '#009E73', 'unreliable': '#D55E00', 'points': '#888888', - }, - 'Brewer Set2': { - 'palette': ['#66C2A5', '#FC8D62', '#8DA0CB', '#E78AC3', '#A6D854', '#FFD92F', '#E5C494', '#B3B3B3'], - 'highlight': '#FFD400', 'base': '#8DA0CB', 'reliable': '#66C2A5', 'unreliable': '#FC8D62', 'points': '#888888', - }, - 'Viridis': { - 'palette': ['#440154', '#472D7B', '#3B528B', '#2C728E', '#21918C', '#5EC962', '#8FD744', '#FDE725'], - 'highlight': '#FF4B4B', 'base': '#21918C', 'reliable': '#5EC962', 'unreliable': '#FDE725', 'points': '#909090', - }, - 'Grayscale': { - 'palette': ['#111111', '#333333', '#555555', '#777777', '#999999', '#BBBBBB', '#DDDDDD', '#EEEEEE'], - 'highlight': '#FF4B4B', 'base': '#777777', 'reliable': '#111111', 'unreliable': '#BBBBBB', 'points': '#909090', - }, - 'Neon': { - 'palette': ['#FF00A0', '#00E5FF', '#7CFF00', '#FFD000', '#B000FF', '#FF5C00', '#00FF9C', '#3D5AFE'], - 'highlight': '#FFFFFF', 'base': '#00E5FF', 'reliable': '#7CFF00', 'unreliable': '#FF5C00', 'points': '#00FF9C', - }, -} - -# Set default theme to Okabe-Ito -DEFAULT_THEME = 'Okabe-Ito' - -# Single-colour keys to use -ROLE_KEYS = ('highlight', 'base', 'reliable', 'unreliable', 'points') -UI_KEYS = ('scene_bg', 'paper_bg', 'grid', 'font') +from evaluator.utils.figures.theme import THEMES, DEFAULT_THEME, ROLE_KEYS, UI_KEYS +# ==================== +# Define constants +# ==================== # Session constants _SESSION_NAME = 'viewer_theme_name' _SESSION_OVERRIDES = 'viewer_theme_overrides' diff --git a/src/evaluator/config.toml b/src/evaluator/config.toml index ccdcda7..fa5118f 100644 --- a/src/evaluator/config.toml +++ b/src/evaluator/config.toml @@ -10,11 +10,13 @@ format = "json" # Label command default configuration parameters [label] min_arc_coverage = 0.30 -merge_centre_tol_factor = 1.5 +merge_centre_tol_factor = 0.573 merge_radius_tol_pct = 0.30 +merge_concentric_radius_pct = 0.10 +merge_concentric_centre_frac = 0.25 maximum_diameter_nm = 500.0 minimum_diameter_nm = 20.0 -membrane_thickness_nm = 7 +membrane_thickness_nm = 2.6 max_workers = 0 # Model command default configuration parameters @@ -27,6 +29,7 @@ max_workers = 0 # Analyse command default configuration parameters [analyse] fill_threshold = 0.05 +membrane_thickness_nm = "auto" max_workers = 0 qc_max_sphere_rmse_rel = 0.25 qc_max_aspect_ratio = 2.5 @@ -36,8 +39,20 @@ qc_max_fit_points = 4000 # Plot command default configuration parameters [plot] -default_sections = ["distributions", "qc", "scatter"] +default_sections = ["distributions", "qc", "scatter", "correlations"] fill_ratio_flag_threshold = 0.05 +theme = "Okabe-Ito" # figure colours (valid options: Okabe-Ito, Brewer Set2, Viridis, Grayscale, Neon) +bin_width = 25.0 # histogram bin width in feature units (0 for automatic bins) +scatter_pairs = [["equiv_diameter_nm", "lumen_volume"], ["surface_area", "lumen_volume"], ["aspect_ratio", "equiv_diameter_nm"]] + +[plot.bin_width_overrides] +aspect_ratio = 0.1 +eccentricity = 0.05 +closure_fill_ratio = 0.05 +sphere_rmse_rel = 0.02 +membrane_volume = 0.0 +lumen_volume = 0.0 +surface_area = 0.0 # Visualise command default configuration parameters [visualise] diff --git a/src/evaluator/utils/config.py b/src/evaluator/utils/config.py index 0f3a53d..7aa7961 100644 --- a/src/evaluator/utils/config.py +++ b/src/evaluator/utils/config.py @@ -35,11 +35,13 @@ class OutputConfig(_Section): class LabelConfig(_Section): '''Configuration parameters for `evaluator label` command''' min_arc_coverage: float = Field(0.30, description='Minimum fraction of expected spherical surface a component must cover to be retained', gt=0, le=1) - merge_centre_tol_factor: float = Field(1.5, description='Centroids within centre_tol_factor * (r_a + r_b) are merge candidates') + merge_centre_tol_factor: float = Field(0.573, description='Centroids within centre_tol_factor * (r_a + r_b) are merge candidates') merge_radius_tol_pct: float = Field(0.30, description='Maximum relative difference between two components\' estimated radii to be merge candidates', gt=0) + merge_concentric_radius_pct: float = Field(0.10, description='Radii difference threshold to prevent merging concentric pairs', ge=0) + merge_concentric_centre_frac: float = Field(0.25, description='Threshold to consider components concentric if centroids are closer than this fraction of the smaller radius', ge=0) minimum_diameter_nm: float = Field(20, description='Minimum component diameter (nm)') maximum_diameter_nm: float = Field(500, description='Maximum component diameter (nm)') - membrane_thickness_nm: float = Field(7, description='Membrane thickness (nm)') + membrane_thickness_nm: float = Field(2.6, description='Membrane thickness (nm)') max_workers: int | None = Field(None, description='Maximum parallel worker processes for batch input (default: CPU count)') @@ -53,6 +55,7 @@ class ModelConfig(_Section): class AnalyseConfig(_Section): '''Configuration parameters for `evaluator analyse` command''' fill_threshold: float = Field(..., description='Fill ratio threshold for labelling') + membrane_thickness_nm: float | Literal['auto'] = Field('auto', description='Assumed membrane shell thickness (nm), used to derive outer diameter from membrane volume') max_workers: int | None = Field(None, description='Maximum parallel worker processes for batch input (default: CPU count)') # Vesicle-vs-debris QC thresholds qc_max_sphere_rmse_rel: float = Field(0.25, description='Max best-fit-sphere relative RMSE for a vesicle-like component', gt=0) @@ -63,10 +66,12 @@ class AnalyseConfig(_Section): class PlotConfig(_Section): '''Configuration parameters for `evaluator plot` command''' - default_sections: list[Literal['distributions', 'qc', 'scatter', 'concordance', 'compare']] = Field( - default_factory=lambda: ['distributions', 'qc', 'scatter'], description='Sections run when --section/--all are not given', - ) + default_sections: list[Literal['distributions', 'qc', 'scatter', 'correlations', 'concordance', 'reliability', 'agreement', 'compare']] = Field(default_factory=lambda: ['distributions', 'qc', 'scatter', 'correlations'], description='Sections run when --section/--all are not given') fill_ratio_flag_threshold: float = Field(0.05, description='closure_fill_ratio below this is flagged in qc plots') + theme: Literal['Okabe-Ito', 'Brewer Set2', 'Viridis', 'Grayscale', 'Neon'] = Field('Okabe-Ito', description='Colour theme for figures (same themes as the viewer)') + bin_width: float = Field(25.0, ge=0, description='Histogram bin width in feature units, anchored at 0 (0 for automatic bins)') + bin_width_overrides: dict[str, float] = Field(default_factory=lambda: {'aspect_ratio': 0.1, 'eccentricity': 0.05, 'closure_fill_ratio': 0.05, 'sphere_rmse_rel': 0.02, 'membrane_volume': 0.0, 'lumen_volume': 0.0, 'surface_area': 0.0}, description='Per-feature histogram bin widths replacing bin_width (0 for automatic bins)') + scatter_pairs: list[tuple[str, str]] = Field(default_factory=lambda: [('equiv_diameter_nm', 'lumen_volume'), ('surface_area', 'lumen_volume'), ('aspect_ratio', 'equiv_diameter_nm')], description='(x, y) feature pairs drawn by the scatter section') class VisualiseConfig(_Section): '''Configuration parameters for `evaluator visualise` command''' diff --git a/src/evaluator/commands/viewer/utils/format.py b/src/evaluator/utils/figures/format.py similarity index 100% rename from src/evaluator/commands/viewer/utils/format.py rename to src/evaluator/utils/figures/format.py diff --git a/src/evaluator/commands/viewer/utils/plots.py b/src/evaluator/utils/figures/plots.py similarity index 86% rename from src/evaluator/commands/viewer/utils/plots.py rename to src/evaluator/utils/figures/plots.py index 3401062..173af9d 100644 --- a/src/evaluator/commands/viewer/utils/plots.py +++ b/src/evaluator/utils/figures/plots.py @@ -10,10 +10,10 @@ import numpy as np, pandas as pd, plotly.graph_objects as go, scipy.stats as _stats # ==================== -# Import internal viewer utilities +# Import internal figure utilities # ==================== -from evaluator.commands.viewer.utils import theme as themeutil -from evaluator.commands.viewer.utils.format import pretty_column +from evaluator.utils.figures import theme as themeutil +from evaluator.utils.figures.format import pretty_column # ==================== # Theme helpers @@ -331,7 +331,7 @@ def bland_altman(df: pd.DataFrame, col_a: str, col_b: str, selected: set[int]) - for yv, dash in ((bias, 'solid'), (lo, 'dash'), (hi, 'dash')): fig.add_hline(y=yv, line_color=ACTIVE['highlight'], line_dash=dash, line_width=1) _annotate(fig, f'bias {bias:+.3g}
95% limits [{lo:.3g}, {hi:.3g}]
n {int(m.sum())}') - _style(fig, f'Bland-Altman: {pretty_column(col_a)} vs {pretty_column(col_b)}', 'mean of {pretty_column(col_a)} & {pretty_column(col_b)}', 'difference (a - b)') + _style(fig, f'Bland-Altman: {pretty_column(col_a)} vs {pretty_column(col_b)}', f'mean of {pretty_column(col_a)} & {pretty_column(col_b)}', 'difference (a - b)') return fig def correlation_matrix(df: pd.DataFrame, cols: list[str], method: str = 'spearman') -> go.Figure: @@ -361,6 +361,52 @@ def correlation_matrix(df: pd.DataFrame, cols: list[str], method: str = 'spearma fig.update_yaxes(autorange='reversed') return fig +# ==================== +# Define builders for plot-specific figures +# ==================== +def bar(df: pd.DataFrame, x: str, y: str, *, title: str, y_title: str, y_range: tuple[float, float] | None = None) -> go.Figure: + fig = go.Figure(go.Bar(x=df[x], y=df[y], marker_color=ACTIVE['base'])) + _style(fig, title, pretty_column(x), y_title) + fig.update_xaxes(tickangle=45) + if y_range is not None: + fig.update_yaxes(range=list(y_range)) + return fig + +def group_distribution(df: pd.DataFrame, feature: str, group: str = 'group') -> go.Figure: + pal = ACTIVE['palette'] + fig = go.Figure() + for i, g in enumerate(sorted(df[group].dropna().unique())): + fig.add_trace(go.Violin( + y=_numeric(df.loc[df[group] == g, feature]), name=str(g), box_visible=True, meanline_visible=True, + points=False, fillcolor=pal[i % len(pal)], line_color=ACTIVE['font'], opacity=0.7, showlegend=False, + )) + _style(fig, f'{pretty_column(feature)} by {group}', group.title(), pretty_column(feature)) + return fig + +def pca_scatter(df: pd.DataFrame, features: list[str], group: str = 'group') -> go.Figure | None: + '''PC1 vs PC2 of standardised features, or None with < 3 rows or < 2 usable features''' + num = df[features].apply(_numeric) + ok = num.notna().all(axis=1) + x = num[ok] + x = x.loc[:, x.std(ddof=1) > 0] # constant columns cannot be standardised + if len(x) < 3 or x.shape[1] < 2: + return None + z = ((x - x.mean()) / x.std(ddof=1)).to_numpy() + u, s, _ = np.linalg.svd(z, full_matrices=False) + scores, var = u[:, :2] * s[:2], s ** 2 / np.sum(s ** 2) + groups = df.loc[ok, group].to_numpy() + pal = ACTIVE['palette'] + fig = go.Figure() + for i, g in enumerate(sorted(set(groups))): + m = groups == g + fig.add_trace(go.Scattergl( + x=scores[m, 0], y=scores[m, 1], mode='markers', name=str(g), + marker=dict(size=9, color=pal[i % len(pal)], line=dict(width=0)), + )) + _style(fig, 'PCA of standardised features', f'PC1 ({var[0]:.0%})', f'PC2 ({var[1]:.0%})') + fig.update_layout(showlegend=True) + return fig + # ==================== # Define function to extract vesicle labels from selection # ==================== diff --git a/src/evaluator/utils/figures/theme.py b/src/evaluator/utils/figures/theme.py new file mode 100644 index 0000000..6b4e123 --- /dev/null +++ b/src/evaluator/utils/figures/theme.py @@ -0,0 +1,49 @@ +''' +======================================= +EValuator: SHARED FIGURE COLOUR THEMES +======================================= +''' + +# ==================== +# Define constants +# ==================== +# Named themes +THEMES: dict[str, dict] = { + 'Okabe-Ito': { + 'palette': ['#E69F00', '#56B4E9', '#009E73', '#F0E442', '#0072B2', '#D55E00', '#CC79A7', '#999999'], + 'highlight': '#FFD400', 'base': '#56B4E9', 'reliable': '#009E73', 'unreliable': '#D55E00', 'points': '#888888', + }, + 'Brewer Set2': { + 'palette': ['#66C2A5', '#FC8D62', '#8DA0CB', '#E78AC3', '#A6D854', '#FFD92F', '#E5C494', '#B3B3B3'], + 'highlight': '#FFD400', 'base': '#8DA0CB', 'reliable': '#66C2A5', 'unreliable': '#FC8D62', 'points': '#888888', + }, + 'Viridis': { + 'palette': ['#440154', '#472D7B', '#3B528B', '#2C728E', '#21918C', '#5EC962', '#8FD744', '#FDE725'], + 'highlight': '#FF4B4B', 'base': '#21918C', 'reliable': '#5EC962', 'unreliable': '#FDE725', 'points': '#909090', + }, + 'Grayscale': { + 'palette': ['#111111', '#333333', '#555555', '#777777', '#999999', '#BBBBBB', '#DDDDDD', '#EEEEEE'], + 'highlight': '#FF4B4B', 'base': '#777777', 'reliable': '#111111', 'unreliable': '#BBBBBB', 'points': '#909090', + }, + 'Neon': { + 'palette': ['#FF00A0', '#00E5FF', '#7CFF00', '#FFD000', '#B000FF', '#FF5C00', '#00FF9C', '#3D5AFE'], + 'highlight': '#FFFFFF', 'base': '#00E5FF', 'reliable': '#7CFF00', 'unreliable': '#FF5C00', 'points': '#00FF9C', + }, +} + +# Set default theme to Okabe-Ito +DEFAULT_THEME = 'Okabe-Ito' + +# Single-colour keys to use +ROLE_KEYS = ('highlight', 'base', 'reliable', 'unreliable', 'points') +UI_KEYS = ('scene_bg', 'paper_bg', 'grid', 'font') + +# Light-mode plot chrome used outside Streamlit +STATIC_UI = {'scene_bg': '#FFFFFF', 'paper_bg': '#FFFFFF', 'grid': '#E5E5E5', 'font': '#31333F'} + +# ==================== +# Define resolver for use outside Streamlit +# ==================== +def resolve(name: str = DEFAULT_THEME) -> dict: + src = {**THEMES[DEFAULT_THEME], **THEMES[name]} + return {**STATIC_UI, **src, 'palette': list(src['palette'])} diff --git a/src/evaluator/utils/mrc.py b/src/evaluator/utils/mrc.py index 3667a72..f2d73f1 100644 --- a/src/evaluator/utils/mrc.py +++ b/src/evaluator/utils/mrc.py @@ -15,7 +15,13 @@ # ==================== # Import internal dependencies # ==================== -from .settings import lg +from evaluator.utils.settings import lg + +# ==================== +# Define constant: FITTED_MRC_LABEL +# ==================== +FITTED_MRC_LABEL = 'EValuator: reconstructed model output' + # ==================== # Define function: validateMRCFile @@ -42,12 +48,17 @@ def readMRCFile(path: Path): ''' with mrcfile.open(str(path), mode='r', permissive=True) as file: data = file.data.copy() - vox_a = float(file.voxel_size.x) - if vox_a == 0.0: + vox_x, vox_y, vox_z = float(file.voxel_size.x), float(file.voxel_size.y), float(file.voxel_size.z) + if any(v == 0 for v in [vox_x, vox_y, vox_z]): lg.warning(f"{path.name}: voxel size not found in MRC header. Physical measurement units will be voxels.") voxel_size_nm = None else: - voxel_size_nm = vox_a / 10.0 + # Check for anisotropy + if vox_y and abs(vox_y - vox_x) / vox_x > 0.01: + lg.warning(f"{path.name}: anisotropic voxel sizes found (x={vox_x}, y={vox_y} Å). Measurements assume isotropic voxels so will be biased along y.") + if vox_z and abs(vox_z - vox_x) / vox_x > 0.01: + lg.warning(f"{path.name}: anisotropic voxel sizes found (x={vox_x}, z={vox_z} Å). Measurements assume isotropic voxels so will be biased along z.") + voxel_size_nm = vox_x / 10.0 return data, voxel_size_nm # ==================== @@ -63,6 +74,35 @@ def writeMRCFile(data: numpy.ndarray, voxel_size_nm: float | None, path: Path): if voxel_size_nm is not None: vox_a = voxel_size_nm * 10.0 mrc.voxel_size = vox_a +# ==================== +# Define function: writeMRCFile +# ==================== +def writeMRCFile(data: numpy.ndarray, voxel_size_nm: float | None, path: Path, *, label: str | None = None): + ''' + Write a numpy array to an MRC file. + If voxel_size_nm is provided, encodes it in the header (converting nm back to Angstroms). + If label is provided, it is written as the MRC header's first text label (max 80 chars), readable via readMRCLabels. + ''' + with mrcfile.new(str(path), overwrite=True) as mrc: + mrc.set_data(data) + if voxel_size_nm is not None: + vox_a = voxel_size_nm * 10.0 + mrc.voxel_size = vox_a + if label is not None: + mrc.header.label[0] = label.encode('ascii', errors='replace')[:80] + mrc.header.nlabl = 1 + +# ==================== +# Define function: readMRCLabels +# ==================== +def readMRCLabels(path: Path) -> list[str]: + ''' + Return non-empty text labels stored in an MRC file's header. + ''' + with mrcfile.open(str(path), mode='r', permissive=True) as file: + n = int(file.header.nlabl) + raw_labels = file.header.label[:n] + return [bytes(l).decode('ascii', errors='replace').strip('\x00').strip() for l in raw_labels] # ==================== # Define function: labelComponents diff --git a/src/evaluator/utils/r/deps/check_deps.R b/src/evaluator/utils/r/deps/check_deps.R deleted file mode 100644 index 9c38838..0000000 --- a/src/evaluator/utils/r/deps/check_deps.R +++ /dev/null @@ -1,21 +0,0 @@ -# check_deps.R: report which R dependencies are installed (as JSON to stdout) - -script_dir <- local({ - args <- commandArgs(trailingOnly = FALSE) - script <- grep("^--file=", args, value = TRUE) - dirname(normalizePath(sub("^--file=", "", script))) -}) -source(file.path(script_dir, "dependencies.R")) - -all_packages <- c(R_DEPENDENCIES$cran) -installed <- character(0) -missing <- character(0) -for (pkg in all_packages) { - ok <- tryCatch(requireNamespace(pkg, quietly = TRUE), error = function(e) FALSE) - if (ok) installed <- c(installed, pkg) else missing <- c(missing, pkg) -} -cat(sprintf( - '{"installed":[%s],"missing":[%s]}', - paste(sprintf('"%s"', installed), collapse = ","), - paste(sprintf('"%s"', missing), collapse = ",") -)) diff --git a/src/evaluator/utils/r/deps/dependencies.R b/src/evaluator/utils/r/deps/dependencies.R deleted file mode 100644 index 26fa332..0000000 --- a/src/evaluator/utils/r/deps/dependencies.R +++ /dev/null @@ -1,5 +0,0 @@ -# dependencies.R: list of all R packages required by evaluator plot command - -R_DEPENDENCIES <- list( - cran = c("tidyverse", "openxlsx", "viridisLite", "jsonlite") -) diff --git a/src/evaluator/utils/r/deps/install_deps.R b/src/evaluator/utils/r/deps/install_deps.R deleted file mode 100644 index 95b639b..0000000 --- a/src/evaluator/utils/r/deps/install_deps.R +++ /dev/null @@ -1,30 +0,0 @@ -# install_deps.R: install R dependencies for evaluator plot command (only installing missing packages) - -script_dir <- local({ - args <- commandArgs(trailingOnly = FALSE) - script <- grep("^--file=", args, value = TRUE) - dirname(normalizePath(sub("^--file=", "", script))) -}) -source(file.path(script_dir, "dependencies.R")) - -is_missing <- function(pkg) !requireNamespace(pkg, quietly = TRUE) - -missing_cran <- Filter(is_missing, R_DEPENDENCIES$cran) - -if (length(missing_cran) == 0) { - message("All R report dependencies are already installed.") - quit(status = 0) -} - -if (length(missing_cran) > 0) { - message(sprintf("Installing %d CRAN package(s): %s", length(missing_cran), paste(missing_cran, collapse = ", "))) - install.packages(missing_cran, repos = "https://cloud.r-project.org") -} - -still_missing <- Filter(is_missing, c(R_DEPENDENCIES$cran, R_DEPENDENCIES$bioc)) -if (length(still_missing) > 0) { - message(sprintf("Still missing after install attempt: %s", paste(still_missing, collapse = ", "))) - quit(status = 1) -} - -message("All R report dependencies installed successfully.") \ No newline at end of file diff --git a/tests/cli/cli_test.py b/tests/cli/cli_test.py index 9103b38..33be945 100644 --- a/tests/cli/cli_test.py +++ b/tests/cli/cli_test.py @@ -147,8 +147,7 @@ def test_valid_analyse_only_dispatches(self, tmp_path): 'aspect_ratio,eccentricity,membrane_volume,lumen_volume,surface_area,' 'is_enclosed,closure_fill_ratio\ntomo1.mrc,1,80.0,90.0,70.0,1.28,0.5,50000,30000,20000,True,0.4\n' ) - with patch('evaluator.commands.plot.plot.resolve_rscript', return_value='Rscript'), \ - patch('evaluator.commands.plot.plot.dispatch'): + with patch('evaluator.commands.plot.plot._run_section'): result = runner.invoke(evaluator, ['plot', '--analyse', str(analyse), '-o', str(tmp_path), '--section', 'qc']) assert result.exit_code == 0 diff --git a/tests/integration/analyse_integration_test.py b/tests/integration/analyse_integration_test.py index 88e7455..1eaa846 100644 --- a/tests/integration/analyse_integration_test.py +++ b/tests/integration/analyse_integration_test.py @@ -13,11 +13,12 @@ EXPECTED_COLUMNS = [ "tomogram", "label", "equiv_diameter_nm", "major_axis_diameter", "minor_axis_diameter", - "aspect_ratio", "eccentricity", + "xy_diameter_nm", "aspect_ratio", "eccentricity", "membrane_volume", "lumen_volume", "surface_area", "is_enclosed", "closure_fill_ratio", "sphere_rmse_rel", "solidity", "arc_coverage", "bbox_extent", - "is_vesicle_like", "qc_flags", + "is_vesicle_like", "touches_z_boundary", "qc_flags", + "membrane_thickness_nm", "membrane_thickness_source", "voxel_size_nm", "measurement_units", ] @@ -84,13 +85,16 @@ def test_closure_fill_ratio_above_threshold(self, analyse_csv): class TestAnalyseMeasurements: def test_equiv_diameter_within_generator_range(self, analyse_csv): """ - Synthetic EVs have outer diameters 42.9–128.6 nm (r_outer 40–120 vox). - equiv_diameter is derived from shell volume, not geometric diameter, - so it will be somewhat smaller — but must remain within the filter bounds. + Synthetic EVs have outer diameters 42.9-128.6 nm (r_outer 40-120 vox). + equiv_diameter is derived from shell volume inversion using the default membrane thickness of 7nm (which differs slightly from true shell thickness of 5.4nm). """ df = pd.read_csv(analyse_csv) - assert (df["equiv_diameter_nm"] > 20.0).all() - assert (df["equiv_diameter_nm"] < 500.0).all() + assert (df["equiv_diameter_nm"] > 35.0).all() + assert (df["equiv_diameter_nm"] < 150.0).all() + def test_membrane_thickness_columns(self, analyse_csv): + df = pd.read_csv(analyse_csv) + assert (df["membrane_thickness_source"] == "estimated").all() + assert (df["membrane_thickness_nm"] > 0).all() def test_major_axis_geq_minor_axis(self, analyse_csv): df = pd.read_csv(analyse_csv) assert (df["major_axis_diameter"] >= df["minor_axis_diameter"]).all() @@ -111,8 +115,8 @@ def test_measurement_units_are_nm(self, analyse_csv): def test_aspect_ratio_near_one_for_spheres(self, analyse_csv): """Synthetic perfect spheres must have aspect ratios close to 1.0.""" df = pd.read_csv(analyse_csv) - assert (df["aspect_ratio"] > 0.5).all() - assert (df["aspect_ratio"] < 2.0).all() + assert (df["aspect_ratio"] > 0.98).all() + assert (df["aspect_ratio"] < 1.02).all() def test_qc_metric_columns_finite(self, analyse_csv): df = pd.read_csv(analyse_csv) for col in ["sphere_rmse_rel", "solidity", "arc_coverage", "bbox_extent"]: @@ -157,6 +161,16 @@ def test_binary_segmentation_input_accepted(self, seg_path, tmp_path): df = pd.read_csv(csvs[0]) assert len(df) == N_EVS + def test_fitted_mrc_input_skipped_not_measured(self, tmp_path): + from evaluator.utils.mrc import FITTED_MRC_LABEL, writeMRCFile + fitted = np.zeros((10, 10, 10), dtype=np.uint16) + fitted[2:8, 2:8, 2:8] = 1 + fitted_path = tmp_path / "fitted.mrc" + writeMRCFile(fitted, 1.0, fitted_path, label=FITTED_MRC_LABEL) + analyse(fitted_path, tmp_path) + csv_dir = tmp_path / "evaluator" / "analyse" + assert not csv_dir.exists() or list(csv_dir.glob("*.csv")) == [] + def test_directory_input_processes_single_file(self, labelled_path, tmp_path_factory): """ Passing a directory containing one labelled MRC must produce a CSV diff --git a/tests/integration/label_integration_test.py b/tests/integration/label_integration_test.py index 4c630c9..3dc6eb9 100644 --- a/tests/integration/label_integration_test.py +++ b/tests/integration/label_integration_test.py @@ -1,6 +1,5 @@ # -- Import external dependencies ------ -import pytest -import numpy as np +import mrcfile, numpy as np, pytest from pathlib import Path # -- Import internal dependencies ------ @@ -87,6 +86,12 @@ def test_high_min_arc_coverage_excludes_all(self, seg_path, tmp_path): out_file = next(out_dir.glob("*.mrc")) data, _ = readMRCFile(out_file) assert int(data.max()) == 0 + def test_exclusion_tally_reports_size_floor_breakdown(self, seg_path, tmp_path): + label_components(seg_path, tmp_path, min_arc_coverage=1.01) + out_dir = tmp_path / "evaluator" / "label" + tally_file = next(out_dir.glob("*_exclusions.txt")) + tally = tally_file.read_text() + assert "too small to reach threshold" not in tally def test_tight_diameter_range_excludes_all(self, seg_path, tmp_path): label_components(seg_path, tmp_path, minimum_diameter_nm=1e6, maximum_diameter_nm=1e6 + 1) out_dir = tmp_path / "evaluator" / "label" @@ -98,4 +103,37 @@ def test_default_thresholds_retain_all_evs(self, seg_path, tmp_path): out_dir = tmp_path / "evaluator" / "label" out_file = next(out_dir.glob("*.mrc")) data, _ = readMRCFile(out_file) - assert int(data.max()) == N_EVS \ No newline at end of file + assert int(data.max()) == N_EVS + +# -- Define merging relabelling tests +class TestLabelComponentsMergeSurvivesToOutput: + @pytest.fixture + def split_shell_seg_path(self, tmp_path): + # Create a hollow sphere and simulate missing wedge + shape = (80, 80, 80) + centre = np.array([40.0, 40.0, 40.0]) + r_outer, r_inner = 20, 12 + zz, yy, xx = np.indices(shape) + dist = np.sqrt((zz - centre[0]) ** 2 + (yy - centre[1]) ** 2 + (xx - centre[2]) ** 2) + shell = (dist >= r_inner) & (dist <= r_outer) + # Remove the equatorial band to break 6-connectivity + band = np.abs(zz - centre[0]) <= 4 + split = shell & ~band + assert split.sum() > 0 + seg_path = tmp_path / "split_shell_segmentation.mrc" + with mrcfile.new(str(seg_path), overwrite=True) as mrc: + mrc.set_data(split.astype(np.float32)) + mrc.voxel_size = 5.36 + return seg_path + def test_split_vesicle_written_as_single_label(self, split_shell_seg_path, tmp_path): + label_components(split_shell_seg_path, tmp_path, min_arc_coverage=0.0, minimum_diameter_nm=0.0) + out_dir = tmp_path / "evaluator" / "label" + out_file = next(out_dir.glob("*.mrc")) + data, _ = readMRCFile(out_file) + from scipy import ndimage + struc = ndimage.generate_binary_structure(3, 1) + raw_data, _ = readMRCFile(split_shell_seg_path) + _, n_raw_components = ndimage.label(raw_data.astype(bool), structure=struc) + assert n_raw_components == 2, "test fixture must produce two disjoint pieces pre-merge" + assert int(data.max()) == 1 + assert sorted(np.unique(data).tolist()) == [0, 1] diff --git a/tests/integration/model_integration_test.py b/tests/integration/model_integration_test.py index df4a0ac..53d1aaf 100644 --- a/tests/integration/model_integration_test.py +++ b/tests/integration/model_integration_test.py @@ -1,14 +1,13 @@ # ==================== # Import external dependencies # ==================== -import json, tomllib -import pytest +import json, numpy as np, pytest, tomllib # ==================== # Import internal dependencies # ==================== from evaluator.commands.model.model import model_evs -from evaluator.utils.mrc import readMRCFile +from evaluator.utils.mrc import readMRCFile, writeMRCFile # ==================== # Define constants @@ -147,6 +146,17 @@ def test_overrides_recorded_in_params_toml(self, labelled_path, tmp_path): assert data['rmse_relative_max'] == 0.5 assert data['min_points'] == 5 + def test_binary_mask_input_warns(self, tmp_path, monkeypatch): + import evaluator.commands.model.model as model_cmd + binary = np.zeros((20, 20, 20), dtype=np.float32) + binary[5:15, 5:15, 5:15] = 1.0 + mask_path = tmp_path / "binary_mask.mrc" + writeMRCFile(binary, 1.0, mask_path) + warnings = [] + monkeypatch.setattr(model_cmd.lg, "warning", lambda msg: warnings.append(msg)) + model_evs(mask_path, tmp_path) + assert any("only one non-zero label found" in w for w in warnings) + def test_strict_min_points_override_fails_reliability(self, labelled_path, tmp_path): ''' A deliberately unreachable min_points threshold must flip every synthetic EV to unreliable, confirming the override reaches _assess_reliability rather than only being recorded in provenance. diff --git a/tests/integration/plot_integration_test.py b/tests/integration/plot_integration_test.py index beaa616..5285870 100644 --- a/tests/integration/plot_integration_test.py +++ b/tests/integration/plot_integration_test.py @@ -1,18 +1,23 @@ # ==================== # Import external dependencies # ==================== -import csv, pytest, shutil +import csv, pytest from openpyxl import load_workbook # ==================== # Import internal dependencies # ==================== from evaluator.commands.plot.plot import run_plot +from evaluator.commands.plot.utils.export import check_kaleido -pytestmark = pytest.mark.skipif( - shutil.which('Rscript') is None, - reason='Rscript not on PATH', -) +def _kaleido_ok() -> bool: + try: + check_kaleido() + return True + except RuntimeError: + return False + +pytestmark = pytest.mark.skipif(not _kaleido_ok(), reason='kaleido/Chrome unavailable') # ==================== # Define constants @@ -20,21 +25,21 @@ ANALYSE_COLUMNS = [ 'tomogram', 'label', 'equiv_diameter_nm', 'major_axis_diameter', 'minor_axis_diameter', 'aspect_ratio', 'eccentricity', 'membrane_volume', 'lumen_volume', 'surface_area', - 'is_enclosed', 'closure_fill_ratio', + 'is_enclosed', 'closure_fill_ratio', 'sphere_rmse_rel', 'is_vesicle_like', 'qc_flags', ] # ==================== # Define helper functions # ==================== -def _write_analyse_csv(path, n=12): +def _write_analyse_csv(path, n=12, tomogram='tomo1.mrc', shift=0.0): with path.open('w', newline='') as fh: writer = csv.writer(fh) writer.writerow(ANALYSE_COLUMNS) for i in range(n): writer.writerow([ - 'tomo1.mrc', i, 60.0 + i, 70.0 + i, 50.0 + i, 1.2, 0.4, - 40000 + i * 100, 25000 + i * 50, 18000 + i * 80, - bool(i % 3), 0.02 + 0.01 * i, + tomogram, i, 60.0 + i + shift, 70.0 + i + shift, 50.0 + i, 1.2 + 0.01 * i, 0.4 + 0.01 * i, + 40000 + i * 100, 25000 + i * 50 + shift * 100, 18000 + i * 80, + bool(i % 3), 0.02 + 0.01 * i, 0.05 + 0.005 * i, i != 0, 'sphere_rmse' if i == 0 else '', ]) # ==================== @@ -50,9 +55,12 @@ def analyse_csv(tmp_path_factory): def plot_output(analyse_csv, tmp_path_factory): out = tmp_path_factory.mktemp('plot_output') run_plot( - analyse_csv, None, out, - sections=['distributions', 'qc'], all_sections=False, - overwrite=False, rscript=None, + analyse_csv, + None, + out, + sections=['distributions', 'qc', 'scatter', 'correlations'], + all_sections=False, + overwrite=False, ) return out / 'evaluator' / 'plot' @@ -73,13 +81,46 @@ def test_qc_summary_written(self, plot_output): row = list(wb.active.iter_rows(min_row=2, max_row=2, values_only=True))[0] assert row[0] == 'tomo1.mrc' # tomogram column - def test_density_svg_written_per_feature(self, plot_output): - assert (plot_output / 'distributions' / 'equiv_diameter_nm_density.svg').exists() + def test_histogram_svg_written_per_feature(self, plot_output): + assert (plot_output / 'distributions' / 'equiv_diameter_nm_histogram.svg').exists() def test_index_md_lists_both_sections(self, plot_output): index = (plot_output / 'index.md').read_text() assert 'distributions/' in index and 'qc/' in index + def test_sphere_rmse_rel_reported(self, plot_output): + assert (plot_output / 'distributions' / 'sphere_rmse_rel_histogram.svg').exists() + + def test_debris_excluded_by_default(self, plot_output): + wb = load_workbook(plot_output / 'qc' / 'qc_summary.xlsx') + assert list(wb.active.iter_rows(min_row=2, max_row=2, values_only=True))[0][1] == 11 # 12 rows, 1 flagged + + def test_scatter_and_correlations_written(self, plot_output): + assert (plot_output / 'scatter' / 'equiv_diameter_nm_vs_lumen_volume.svg').exists() + assert (plot_output / 'correlations' / 'correlation_spearman.svg').exists() + +class TestPlotMultiRun: + def test_all_samples_reflected_in_per_sample_output(self, tmp_path): + csv_a, csv_b = tmp_path / 'a.csv', tmp_path / 'b.csv' + _write_analyse_csv(csv_a, tomogram='a.mrc') + _write_analyse_csv(csv_b, tomogram='b.mrc') + sheet = tmp_path / 'samples.tsv' + sheet.write_text(f'sample_id\tpath\nsample_a\t{csv_a}\nsample_b\t{csv_b}\n') + out = tmp_path / 'out' + run_plot(sheet, None, out, sections=['distributions'], all_sections=False, overwrite=False) + section_dir = out / 'evaluator' / 'plot' / 'distributions' + for sample_id in ('sample_a', 'sample_b'): + assert (section_dir / sample_id / 'summary_stats.xlsx').exists() + assert (section_dir / sample_id / 'equiv_diameter_nm_histogram.svg').exists() + + def test_include_flagged_keeps_debris(self, tmp_path): + csv_a = tmp_path / 'a.csv' + _write_analyse_csv(csv_a) + out = tmp_path / 'out' + run_plot(csv_a, None, out, sections=['qc'], all_sections=False, overwrite=False, include_flagged=True) + wb = load_workbook(out / 'evaluator' / 'plot' / 'qc' / 'qc_summary.xlsx') + assert list(wb.active.iter_rows(min_row=2, max_row=2, values_only=True))[0][1] == 12 + class TestPlotConcordance: def test_concordance_requires_matching_rows(self, tmp_path, analyse_csv): model_json = tmp_path / 'model.json' @@ -89,10 +130,59 @@ def test_concordance_requires_matching_rows(self, tmp_path, analyse_csv): "'reliability': {'is_reliable': true}}]}" ) out = tmp_path / 'out' - run_plot( - analyse_csv, model_json, out, - sections=['concordance'], all_sections=False, - overwrite=False, rscript=None, - ) + run_plot(analyse_csv, model_json, out, sections=['concordance'], all_sections=False, overwrite=False) section_dir = out / 'evaluator' / 'plot' / 'concordance' assert not (section_dir / 'concordance_stats.xlsx').exists() + + def test_concordance_and_reliability_with_matching_rows(self, tmp_path): + csv_a = tmp_path / 'a.csv' + _write_analyse_csv(csv_a, n=6) + model = tmp_path / 'model.json' + results = [ + {'label_id': i, 'source_file': 'tomo1.mrc', 'chosen_model': 'sphere', 'radius': (60 + i) / 2, 'rmse_nm': 1.0 + i, + 'reliability': {'is_reliable': bool(i % 2)}} + for i in range(1, 6) + ] + model.write_text(__import__('json').dumps({'results': results, 'parameters': {}})) + out = tmp_path / 'out' + run_plot(csv_a, model, out, sections=['concordance', 'reliability'], all_sections=False, overwrite=False) + base = out / 'evaluator' / 'plot' + assert (base / 'concordance' / 'model_vs_equiv_diameter_nm.svg').exists() + assert (base / 'concordance' / 'concordance_stats.xlsx').exists() + assert (base / 'reliability' / 'reliability.svg').exists() + + def test_agreement_plots_requested_pairs(self, tmp_path): + csv_a = tmp_path / 'a.csv' + _write_analyse_csv(csv_a, n=6) + model = tmp_path / 'model.json' + results = [ + {'label_id': i, 'source_file': 'tomo1.mrc', 'chosen_model': 'sphere', 'radius': (60 + i) / 2 + 0.3 * i, 'rmse_nm': 1.0, + 'reliability': {'is_reliable': True}} + for i in range(1, 6) + ] + model.write_text(__import__('json').dumps({'results': results, 'parameters': {}})) + out = tmp_path / 'out' + run_plot( + csv_a, model, out, sections=['qc'], all_sections=False, overwrite=False, + agreement=[('equiv_diameter_nm', 'model_diameter_nm'), ('major_axis_diameter', 'model_diameter_nm')], + ) + base = out / 'evaluator' / 'plot' / 'agreement' + assert (base / 'equiv_diameter_nm_vs_model_diameter_nm.svg').exists() + assert (base / 'major_axis_diameter_vs_model_diameter_nm.svg').exists() + assert (base / 'agreement_stats.xlsx').exists() + +class TestPlotCompare: + def test_compare_writes_group_outputs(self, tmp_path): + a, b = tmp_path / 'a.csv', tmp_path / 'b.csv' + _write_analyse_csv(a, tomogram='a.mrc') + _write_analyse_csv(b, tomogram='b.mrc', shift=5.0) + c = tmp_path / 'c.csv' + _write_analyse_csv(c, tomogram='c.mrc', shift=12.0) # pca needs three tomograms varying in two or more features + sheet = tmp_path / 'samples.tsv' + sheet.write_text(f'sample_id\tpath\tgroup\ns1\t{a}\tcontrol\ns2\t{b}\ttreated\ns3\t{c}\ttreated\n') + out = tmp_path / 'out' + run_plot(sheet, None, out, sections=['compare'], all_sections=False, overwrite=False) + base = out / 'evaluator' / 'plot' / 'compare' + assert (base / 'violin_equiv_diameter_nm.svg').exists() + assert (base / 'pca.svg').exists() + assert (base / 'group_summary.xlsx').exists() diff --git a/tests/unit/analyse_unit_test.py b/tests/unit/analyse_unit_test.py index 93b5b7b..e7102e1 100644 --- a/tests/unit/analyse_unit_test.py +++ b/tests/unit/analyse_unit_test.py @@ -7,7 +7,7 @@ from evaluator.commands.analyse.utils.filtering import checkEnclosed, morphologicalDilation from evaluator.commands.analyse.utils.geometry import arcCoverage, deriveAxes, measureAxes, shellVolume, sphereFitResidual from evaluator.commands.analyse.utils.io import saveResultsCSV -from evaluator.commands.analyse.utils.measurement import classifyVesicle, computeSurfaceArea, measureEccentricityAspectRatio, measureLumenVolume, measureMembraneVolumeDiameter +from evaluator.commands.analyse.utils.measurement import classifyVesicle, computeSurfaceArea, estimateMembraneThickness, measureEccentricityAspectRatio, measureLumenVolume, measureMembraneVolumeDiameter, measureXYDiameter, touchesZBoundary from tests.unit.conftest import R_INNER, R_OUTER, VOX_NM_SMALL # -- Define shell volume test ---------- @@ -49,17 +49,45 @@ def test_fills_single_missing_voxel(self, hollow_sphere_vol): closed = morphologicalDilation(punctured) assert closed[40, 40, 20], "Gap voxel was not restored by morphological dilation" def test_all_original_voxels_preserved(self, hollow_sphere_vol): - """Closing must not remove any pre-existing True voxels.""" + ''' + Closing must not remove any pre-existing True voxels. + ''' closed = morphologicalDilation(hollow_sphere_vol) - assert np.all(closed[hollow_sphere_vol]), ( - "morphologicalDilation removed voxels that were present in the original mask" - ) + assert np.all(closed[hollow_sphere_vol]), "morphologicalDilation removed voxels that were present in the original mask" def test_returns_bool_array(self, hollow_sphere_vol): closed = morphologicalDilation(hollow_sphere_vol) assert closed.dtype == bool or closed.dtype == np.bool_ def test_same_shape(self, hollow_sphere_vol): closed = morphologicalDilation(hollow_sphere_vol) assert closed.shape == hollow_sphere_vol.shape + def test_dilation_not_clipped_at_box_face(self): + ''' + Components touching mask boundary should be dilated outward instead of being clipped at the box face. + ''' + vol = np.zeros((10, 10, 10), dtype=bool) + vol[0, 5, 5] = True + dilated = morphologicalDilation(vol) + assert dilated[2, 5, 5], "Inward dilation did not reach the expected extent" + def test_symmetric_dilation_extent_matches_interior_point(self): + ''' + The mask boundary should not asymmetrically clip results. + ''' + boundary_vol = np.zeros((10, 10, 10), dtype=bool) + boundary_vol[0, 5, 5] = True + interior_vol = np.zeros((10, 10, 10), dtype=bool) + interior_vol[5, 5, 5] = True + dilated_boundary = morphologicalDilation(boundary_vol) + dilated_interior = morphologicalDilation(interior_vol) + # count of dilated voxels on the non-clipped (y, x) plane through the seed voxel must match + assert dilated_boundary[0, :, :].sum() == dilated_interior[5, :, :].sum() + +# -- Define touches-z-boundary test -------- +class TestTouchesZBoundary: + def test_component_away_from_boundary_not_flagged(self, sphere_regionprops): + assert touchesZBoundary(sphere_regionprops, volume_shape=(200, 200, 200)) is False + def test_component_touching_last_z_slice_flagged(self, sphere_regionprops): + # Shrink the volume's z-extent so sphere's max-z edge sits on the last slice + assert touchesZBoundary(sphere_regionprops, volume_shape=(sphere_regionprops.bbox[3], 200, 200)) is True # -- Define check enclosed test -------- class TestCheckEnclosed: @@ -111,14 +139,54 @@ def test_volume_scales_cubically(self, sphere_regionprops): def test_diameter_is_positive(self, sphere_regionprops): _, diam = measureMembraneVolumeDiameter(sphere_regionprops, scale=VOX_NM_SMALL) assert diam > 0 - def test_diameter_plausible_for_sphere(self, sphere_regionprops): - """ - equiv_diameter is derived from shell volume, not geometric diameter - For r_outer=20, r_inner=12 at scale=1: shell_vol ≈ 26,266 vox. - equiv_diam = (6 x 26266 / π)^(1/3) ≈ 36.9 vox - """ - _, diam = measureMembraneVolumeDiameter(sphere_regionprops, scale=1.0) - assert 30.0 < diam < 45.0, f"Unexpected equiv_diameter {diam:.2f} vox" + def test_diameter_matches_true_outer_diameter_given_correct_thickness(self, sphere_regionprops): + ''' + Synthetic shell has r_outer=20, r_inner=12, thickness=8 vox. + Shell-volume-inversion formula should recover the true outer diameter (2*r_outer=40 vox). + ''' + _, diam = measureMembraneVolumeDiameter(sphere_regionprops, scale=1.0, membrane_thickness_nm=8.0) + assert 38.0 < diam < 42.0, f"Unexpected outer diameter {diam:.2f} vox (expected ~40 vox)" + def test_solid_sphere_fallback_reproduces_old_biased_value_when_no_thickness(self, sphere_regionprops): + vol_nm3, diam = measureMembraneVolumeDiameter(sphere_regionprops, scale=1.0, membrane_thickness_nm=0.0) + assert np.isclose(diam, (6 * vol_nm3 / np.pi) ** (1 / 3)) + def test_thicker_assumed_membrane_gives_smaller_outer_diameter(self, sphere_regionprops): + _, diam_thin = measureMembraneVolumeDiameter(sphere_regionprops, scale=1.0, membrane_thickness_nm=4.0) + _, diam_thick = measureMembraneVolumeDiameter(sphere_regionprops, scale=1.0, membrane_thickness_nm=12.0) + assert diam_thick < diam_thin + +# -- Define xy diameter test -------- +class TestMeasureXYDiameter: + def _shell(self, radius=20, size=60): + zz, yy, xx = np.indices((size, size, size)) - size // 2 + dist = np.sqrt(zz ** 2 + yy ** 2 + xx ** 2) + return (dist >= radius - 2) & (dist <= radius) + def test_sphere_matches_outer_diameter(self): + d = measureXYDiameter(self._shell(), scale=1.0) + assert abs(d - 41) <= 2 + def test_scales_by_voxel_size(self): + shell = self._shell() + assert measureXYDiameter(shell, scale=0.5) == pytest.approx(measureXYDiameter(shell, scale=1.0) / 2) + def test_unchanged_when_poles_clipped(self): + shell = self._shell() + clipped = shell[15:45] # keeps the equator, removes both polar caps + assert measureXYDiameter(clipped, scale=1.0) == measureXYDiameter(shell, scale=1.0) + +# -- Define membrane thickness estimate test ----- +class TestEstimateMembraneThickness: + def test_recovers_known_shell_thickness(self, hollow_sphere_vol): + # fixture shell is 8 vox thick; the ridge estimate reads up to ~half a voxel high + est = estimateMembraneThickness(hollow_sphere_vol, voxel_size_nm=None) + assert abs(est - 8.0) <= 1.0 + def test_scales_with_voxel_size(self, hollow_sphere_vol): + est_vox = estimateMembraneThickness(hollow_sphere_vol, voxel_size_nm=None) + est_nm = estimateMembraneThickness(hollow_sphere_vol, voxel_size_nm=VOX_NM_SMALL) + assert np.isclose(est_nm, est_vox * VOX_NM_SMALL, rtol=1e-6) + def test_empty_mask_returns_none(self): + assert estimateMembraneThickness(np.zeros((20, 20, 20), dtype=bool), voxel_size_nm=1.0) is None + def test_too_few_ridge_voxels_returns_none(self): + mask = np.zeros((20, 20, 20), dtype=bool) + mask[10, 10, 10] = True + assert estimateMembraneThickness(mask, voxel_size_nm=1.0) is None # -- Define lumen measurement test ----- class TestMeasureLumenVolume: @@ -172,35 +240,53 @@ def test_empty_volume_returns_nan(self): # -- Define axis derivation test ------- class TestDeriveAxes: - def test_identity_tensor_returns_ones(self): - """Identity inertia tensor → all eigenvalues=1 → inv_sqrt=[1,1,1]""" + def test_identity_tensor_returns_equal_axes(self): + """Identity inertia tensor (i.e. isotropic) gives equal relative axis lengths""" result = deriveAxes(np.eye(3)) - assert np.allclose(result, [1.0, 1.0, 1.0]) - def test_diagonal_tensor_known_values(self): + assert np.allclose(result, result[0]) + def test_known_ellipsoid_semi_axes_recovered(self): """ - Diagonal tensor [1, 4, 9] → eigvalsh ascending [1,4,9] - → inv_sqrt = [1/1, 1/2, 1/3]. + Solid ellipsoid with known semi-axes a=3,b=2,c=1 (mass=1) has principal moments I_a=(b^2+c^2)/5, I_b=(a^2+c^2)/5, I_c=(a^2+b^2)/5. + deriveAxes must recover axis lengths proportional to the true (3, 2, 1) ratio. """ - tensor = np.diag([1.0, 4.0, 9.0]) + a, b, c = 3.0, 2.0, 1.0 + i_a = (b ** 2 + c ** 2) / 5.0 + i_b = (a ** 2 + c ** 2) / 5.0 + i_c = (a ** 2 + b ** 2) / 5.0 + tensor = np.diag(sorted([i_a, i_b, i_c])) # eigvalsh returns ascending order result = deriveAxes(tensor) - assert np.allclose(result, [1.0, 0.5, 1.0 / 3.0], rtol=1e-6) - def test_zero_eigenvalue_gives_zero_not_inf(self): - """A zero eigenvalue must produce 0.0 in inv_sqrt, not NaN or inf""" - tensor = np.diag([0.0, 1.0, 4.0]) + result_ratio = result / result[-1] + true_ratio = np.array([a, b, c]) / c + assert np.allclose(result_ratio, true_ratio, rtol=1e-6), f"Recovered axis ratio {result_ratio} does not match true ratio {true_ratio}" + def test_zero_eigenvalue_gives_zero_not_nan_or_inf(self): + """A zero eigenvalue must produce 0.0, not NaN or inf.""" + tensor = np.diag([0.0, 4.0, 4.0]) result = deriveAxes(tensor) - assert result[0] == 0.0 assert not np.any(np.isnan(result)) assert not np.any(np.isinf(result)) def test_output_is_descending(self): - """inv_sqrt must be in descending order (largest first)""" - tensor = np.diag([1.0, 4.0, 9.0]) + """Axis lengths must be in descending order (largest first)""" + a, b, c = 3.0, 2.0, 1.0 + i_a, i_b, i_c = (b ** 2 + c ** 2) / 5.0, (a ** 2 + c ** 2) / 5.0, (a ** 2 + b ** 2) / 5.0 + tensor = np.diag(sorted([i_a, i_b, i_c])) result = deriveAxes(tensor) assert result[0] >= result[1] >= result[2] def test_symmetric_tensor_gives_equal_axes(self): - """An isotropic (spherical) inertia tensor must give equal inv_sqrt values""" + """An isotropic (spherical) inertia tensor must give equal axis lengths""" tensor = 5.0 * np.eye(3) result = deriveAxes(tensor) assert np.allclose(result[0], result[1]) and np.allclose(result[1], result[2]) + def test_no_artificial_oblate_ceiling(self): + """ + A strongly oblate ellipsoid (a=b=5, c=1) must give an axis ratio beyond inverse square root's ~1.414 ceiling. + """ + a, b, c = 5.0, 5.0, 1.0 + i_a = (b ** 2 + c ** 2) / 5.0 + i_b = (a ** 2 + c ** 2) / 5.0 + i_c = (a ** 2 + b ** 2) / 5.0 + tensor = np.diag(sorted([i_a, i_b, i_c])) + result = deriveAxes(tensor) + assert (result[0] / result[-1]) > 2.0 # -- Define axes measurements test ----- class TestMeasureAxes: @@ -219,6 +305,36 @@ def test_both_axes_positive(self, sphere_regionprops): _, equiv_diam = measureMembraneVolumeDiameter(sphere_regionprops, scale=1.0) major, minor = measureAxes(sphere_regionprops, equiv_diam) assert major > 0 and minor > 0 + def test_known_ellipsoid_axis_ratio_recovered(self): + """ + Measurement of ellipsoid with semi-axes 12, 6, 3 vox recovers major/minor axis ratio close to 4.0. + """ + from skimage import measure as skmeasure + shape = (30, 30, 30) + zz, yy, xx = np.indices(shape) + cz, cy, cx = 15, 15, 15 + a, b, c = 12.0, 6.0, 3.0 + vol = (((zz - cz) / a) ** 2 + ((yy - cy) / b) ** 2 + ((xx - cx) / c) ** 2) <= 1.0 + props = skmeasure.regionprops(vol.astype(np.int32))[0] + _, equiv_diam = measureMembraneVolumeDiameter(props, scale=1.0) + major, minor = measureAxes(props, equiv_diam) + assert np.isclose(major / minor, a / c, rtol=0.15), f"Recovered axis ratio {major / minor:.2f} too far from true ratio {a / c:.2f}" + def test_single_voxel_component_returns_nan_not_crash(self): + """A degenerate single-voxel component returns NaN.""" + from skimage import measure as skmeasure + vol = np.zeros((5, 5, 5), dtype=np.int32) + vol[2, 2, 2] = 1 + props = skmeasure.regionprops(vol)[0] + major, minor = measureAxes(props, equiv_diameter_nm=10.0) + assert np.isnan(major) and np.isnan(minor) + def test_line_component_returns_nan_not_crash(self): + """A 1-voxel-wide line component returns NaN.""" + from skimage import measure as skmeasure + vol = np.zeros((10, 5, 5), dtype=np.int32) + vol[2:8, 2, 2] = 1 + props = skmeasure.regionprops(vol)[0] + major, minor = measureAxes(props, equiv_diameter_nm=10.0) + assert np.isnan(major) and np.isnan(minor) # -- Define eccentricity & aspect ratio test class TestMeasureEccentricityAspectRatio: @@ -307,6 +423,9 @@ def test_threshold_override_flips_result(self): args = (5000, 0.1, 2.0, 0.3, 0.9, True) assert classifyVesicle(*args, _QC_DEFAULTS)[0] is True assert classifyVesicle(*args, {**_QC_DEFAULTS, "qc_max_aspect_ratio": 1.5})[0] is False + def test_z_boundary_not_a_classification_input(self): + import inspect + assert 'touches_z_boundary' not in inspect.signature(classifyVesicle).parameters # -- Define results saving test -------- class TestSaveResultsCSV: diff --git a/tests/unit/config_unit_test.py b/tests/unit/config_unit_test.py index 72e2f72..095e84b 100644 --- a/tests/unit/config_unit_test.py +++ b/tests/unit/config_unit_test.py @@ -48,7 +48,11 @@ def _write_config(path: Path, content: str | None = None) -> Path: # Define tests # ==================== class TestSchema: - '''The Pydantic Config model and the packaged default are in lockstep''' + '''The Pydantic Config model and the packaged default match''' + def test_merge_centre_tol_factor_below_neighbour_ceiling(self): + from evaluator.utils.config import LabelConfig + assert LabelConfig().merge_centre_tol_factor <= 0.573 + def test_packaged_default_validates(self): '''The packaged config.toml must validate against Config''' from importlib.resources import files @@ -437,7 +441,7 @@ def bad_prompt(key, default, **kwargs): class TestPlotConfig: def test_default_construction(self): cfg = PlotConfig() - assert cfg.default_sections == ['distributions', 'qc', 'scatter'] + assert cfg.default_sections == ['distributions', 'qc', 'scatter', 'correlations'] assert cfg.fill_ratio_flag_threshold == 0.05 def test_extra_key_forbidden(self): @@ -454,7 +458,7 @@ def test_packaged_config_toml_has_plot_section(self): import tomllib data = tomllib.loads((files('evaluator') / 'config.toml').read_text(encoding='utf-8')) config = Config.model_validate(data) - assert config.plot.default_sections == ['distributions', 'qc', 'scatter'] + assert config.plot.default_sections == ['distributions', 'qc', 'scatter', 'correlations'] class TestQCAndReliabilityConfig: def test_analyse_qc_defaults(self): diff --git a/tests/unit/figures_unit_test.py b/tests/unit/figures_unit_test.py new file mode 100644 index 0000000..2c1b3e8 --- /dev/null +++ b/tests/unit/figures_unit_test.py @@ -0,0 +1,29 @@ +''' +Unit tests for the shared figure package +''' + +# ==================== +# Import external dependencies +# ==================== +import subprocess, sys + +# ==================== +# Import evaluator figure helpers +# ==================== +from evaluator.utils.figures import theme as figtheme + +# ==================== +# Define tests +# ==================== +class TestFigures: + def test_plots_import_without_streamlit(self): + code = "import sys; sys.modules['streamlit'] = None; import evaluator.utils.figures.plots" + subprocess.run([sys.executable, '-c', code], check=True) + + def test_resolve_default_theme_has_ui_and_roles(self): + theme = figtheme.resolve() + assert set(figtheme.UI_KEYS) <= set(theme) and set(figtheme.ROLE_KEYS) <= set(theme) + assert theme['paper_bg'] == '#FFFFFF' + + def test_resolve_named_theme(self): + assert figtheme.resolve('Neon')['base'] == figtheme.THEMES['Neon']['base'] diff --git a/tests/unit/label_unit_test.py b/tests/unit/label_unit_test.py index 523acc4..41676e3 100644 --- a/tests/unit/label_unit_test.py +++ b/tests/unit/label_unit_test.py @@ -1,14 +1,14 @@ # ==================== # Import external dependencies # ==================== -import numpy as np -import pytest +import numpy as np, pytest # ==================== # Import internal dependencies # ==================== from evaluator.commands.label.utils.geometric_proxies import estimateCentroidRadius, estimateArcCoverage from evaluator.commands.label.utils.merge import findMergeGroups, applyMerges +from evaluator.utils.config import LabelConfig # ==================== # Define point-cloud generator helpers @@ -86,8 +86,9 @@ def test_split_pair_same_source_merges(self): assert sorted(groups[0]) == [1, 2] def test_separate_vesicles_not_merged(self): - pts_a = _sphere_points(centre=(0, 0, 0), radius=5.0) - pts_b = _sphere_points(centre=(500, 500, 500), radius=20.0) + pts_a = _sphere_points(centre=(0, 0, 0), radius=10.0) + # centre_tol_factor=1.5 * (10+10) = 30 is the merge boundary + pts_b = _sphere_points(centre=(31, 0, 0), radius=10.0) groups = findMergeGroups({1: pts_a, 2: pts_b}, centre_tol_factor=1.5, radius_tol_pct=0.30) assert sorted(sorted(g) for g in groups) == [[1], [2]] @@ -100,6 +101,37 @@ def test_chain_merge_transitive(self): ) assert sorted(groups[0]) == [1, 2, 3] + def test_concentric_pair_not_merged(self): + # 23% radius difference passes merge_radius_tol_pct (0.30) but is concentric + outer = _sphere_points(centre=(0, 0, 0), radius=13.0) + inner = _sphere_points(centre=(0, 0, 0), radius=10.0) + groups = findMergeGroups({1: outer, 2: inner}, centre_tol_factor=1.5, radius_tol_pct=0.30) + assert sorted(sorted(g) for g in groups) == [[1], [2]] + + def test_concentric_same_radius_still_merges(self): + a = _sphere_points(centre=(0, 0, 0), radius=10.0) + b = _sphere_points(centre=(0, 0, 0), radius=10.0) + groups = findMergeGroups({1: a, 2: b}, centre_tol_factor=1.5, radius_tol_pct=0.30) + assert sorted(groups[0]) == [1, 2] + + def test_concentric_guard_can_be_relaxed(self): + outer = _sphere_points(centre=(0, 0, 0), radius=13.0) + inner = _sphere_points(centre=(0, 0, 0), radius=10.0) + groups = findMergeGroups({1: outer, 2: inner}, centre_tol_factor=1.5, radius_tol_pct=0.30, concentric_radius_pct=0.50) + assert sorted(groups[0]) == [1, 2] + + def test_default_centre_tol_does_not_fuse_neighbours(self): + pts_a = _sphere_points(centre=(0, 0, 0), radius=5.0) + pts_b = _sphere_points(centre=(7, 0, 0), radius=5.0) + groups = findMergeGroups({1: pts_a, 2: pts_b}, centre_tol_factor=LabelConfig().merge_centre_tol_factor, radius_tol_pct=0.30) + assert sorted(sorted(g) for g in groups) == [[1], [2]] + + def test_concentric_centre_frac_zero_disables_guard(self): + outer = _sphere_points(centre=(0, 0, 0), radius=13.0) + inner = _sphere_points(centre=(0, 0, 0), radius=10.0) + groups = findMergeGroups({1: outer, 2: inner}, centre_tol_factor=1.5, radius_tol_pct=0.30, concentric_centre_frac=0.0) + assert sorted(groups[0]) == [1, 2] + class TestApplyMerges: def test_voxel_count_preserved(self): vol = np.zeros((10, 10, 10), dtype=np.int32) diff --git a/tests/unit/model_unit_test.py b/tests/unit/model_unit_test.py index 37c4f47..74e5b24 100644 --- a/tests/unit/model_unit_test.py +++ b/tests/unit/model_unit_test.py @@ -11,6 +11,7 @@ fit_sphere_least_squares, fit_ellipsoid, fit_vesicle, + _beam_axis_flag, ) from evaluator.commands.model.utils.reconstruction import build_fitted_mrc @@ -23,24 +24,24 @@ def _sphere_points(centre=(50.0, 50.0, 50.0), radius=20.0, n=800, noise=0.0, see i = np.arange(n) phi = np.arccos(1 - 2 * (i + 0.5) / n) # polar angle: 0 -> pi (full sphere) theta = np.pi * (1 + 5 ** 0.5) * i # golden angle azimuth - x = radius * np.sin(phi) * np.cos(theta) - y = radius * np.sin(phi) * np.sin(theta) z = radius * np.cos(phi) - pts = np.stack([x, y, z], axis=1) + np.asarray(centre) + y = radius * np.sin(phi) * np.sin(theta) + x = radius * np.sin(phi) * np.cos(theta) + pts = np.stack([z, y, x], axis=1) + np.asarray(centre) if noise: pts = pts + rng.normal(0.0, noise, size=pts.shape) return pts -# -- _ellipsoid_points: n points over the full ellipsoid surface, optionally rotated +# -- _ellipsoid_points: n points over the full ellipsoid surface, optionally rotated. radii is given as (r_z, r_y, r_x). def _ellipsoid_points(centre=(50.0, 50.0, 50.0), radii=(10.0, 20.0, 30.0), rotation=None, n=800, noise=0.0, seed=0): rng = np.random.default_rng(seed) i = np.arange(n) phi = np.arccos(1 - 2 * (i + 0.5) / n) theta = np.pi * (1 + 5 ** 0.5) * i local = np.stack([ - radii[0] * np.sin(phi) * np.cos(theta), + radii[0] * np.cos(phi), radii[1] * np.sin(phi) * np.sin(theta), - radii[2] * np.cos(phi), + radii[2] * np.sin(phi) * np.cos(theta), ], axis=1) if rotation is not None: local = local @ rotation.T @@ -49,28 +50,17 @@ def _ellipsoid_points(centre=(50.0, 50.0, 50.0), radii=(10.0, 20.0, 30.0), rotat pts = pts + rng.normal(0.0, noise, size=pts.shape) return pts -# -- _rotation_about_x: right-handed rotation matrix about the x-axis (index 0) -def _rotation_about_x(angle_deg): - a = np.radians(angle_deg) - return np.array([ - [1, 0, 0], - [0, np.cos(a), -np.sin(a)], - [0, np.sin(a), np.cos(a)], - ]) - -# -- _narrow_band_points: points confined to a ~20 degree latitude band, -# simulating a vesicle with most of its surface missing +# -- _narrow_band_points: points confined to a ~20 degree latitude band about beam, simulating a vesicle with most of its surface missing def _narrow_band_points(centre=(50.0, 50.0, 50.0), radius=20.0, n=200, seed=0): rng = np.random.default_rng(seed) phi = np.radians(rng.uniform(80.0, 100.0, size=n)) # polar angle near equator theta = rng.uniform(0.0, 2 * np.pi, size=n) - x = radius * np.sin(phi) * np.cos(theta) - y = radius * np.sin(phi) * np.sin(theta) z = radius * np.cos(phi) - return np.stack([x, y, z], axis=1) + np.asarray(centre) + y = radius * np.sin(phi) * np.sin(theta) + x = radius * np.sin(phi) * np.cos(theta) + return np.stack([z, y, x], axis=1) + np.asarray(centre) -# -- _flat_ring_points: points on a flat annulus (constant z), a degenerate -# input for ellipsoid fitting +# -- _flat_ring_points: points on a flat annulus (constant z), a degenerate input for ellipsoid fitting def _flat_ring_points(centre=(50.0, 50.0, 50.0), radius=15.0, n=40): theta = np.linspace(0.0, 2 * np.pi, n, endpoint=False) x = radius * np.cos(theta) @@ -138,22 +128,16 @@ def test_perfect_sphere_selects_sphere(self): assert result['orientation'] is None def test_elongated_off_axis_selects_ellipsoid(self): ''' - Major axis (radius 30) rotated 90 degrees about x so it points - along global y, i.e. 90 degrees from the beam (z) axis - well - outside the beam-axis guard tolerance. + Major axis (radius 30) placed along global x, i.e. 90 degrees from the beam (z) axis - well outside the beam-axis guard tolerance. ''' - pts = _ellipsoid_points( - radii=(10.0, 10.0, 30.0), - rotation=_rotation_about_x(90.0), - noise=0.02, seed=5, - ) + pts = _ellipsoid_points(radii=(10.0, 10.0, 30.0), noise=0.02, seed=5) result = fit_vesicle(pts) assert result['chosen_model'] == 'ellipsoid' assert result['radii'] is not None assert result['orientation'] is not None def test_elongated_along_beam_axis_selects_sphere_beam_axis(self): - '''Major axis (radius 30) left aligned with global z (the beam axis).''' - pts = _ellipsoid_points(radii=(10.0, 10.0, 30.0), noise=0.02, seed=6) + '''Major axis (radius 30) aligned with global z (beam axis).''' + pts = _ellipsoid_points(radii=(30.0, 10.0, 10.0), noise=0.02, seed=6) result = fit_vesicle(pts) assert result['chosen_model'] == 'sphere (beam-axis)' assert result['radii'] is None @@ -180,6 +164,18 @@ def test_output_contains_expected_keys(self): 'beam_axis', 'sphere_fit', 'ellipsoid_fit', } assert expected_keys.issubset(result.keys()) + def test_beam_axis_flag_uses_index_0_not_2(self): + ''' + major_axis_angle_from_z_deg should be measured against column 0 (z), so an ellipsoid genuinely elongated along z should be flagged (angle near 0), while one elongated along x (90 degrees off-beam) should not be. + ''' + result_z = fit_ellipsoid(_ellipsoid_points(radii=(30.0, 10.0, 10.0), noise=0.02, seed=13)) + flag_z = _beam_axis_flag(result_z['evecs'], result_z['radii'], tol_deg=25.0) + assert flag_z['major_axis_angle_from_z_deg'] < 10.0 + assert flag_z['beam_axis_flagged'] is True + result_x = fit_ellipsoid(_ellipsoid_points(radii=(10.0, 10.0, 30.0), noise=0.02, seed=14)) + flag_x = _beam_axis_flag(result_x['evecs'], result_x['radii'], tol_deg=25.0) + assert flag_x['major_axis_angle_from_z_deg'] > 80.0 + assert flag_x['beam_axis_flagged'] is False # ==================== @@ -214,6 +210,22 @@ def test_span_threshold_is_configurable(self): assert fit_vesicle(pts)['reliability']['span_ok'] is False loosened = fit_vesicle(pts, min_latitude_span_deg=5.0)['reliability'] assert loosened['span_ok'] is True + def test_missing_wedge_caps_fail_span_gate_not_180(self): + '''A shell with missing wedge shouldn't report a full 180 degree latitude span.''' + rng = np.random.default_rng(15) + n = 400 + phi = np.radians(rng.uniform(70.0, 110.0, size=n)) + theta = rng.uniform(0.0, 2 * np.pi, size=n) + radius = 20.0 + z = radius * np.cos(phi) + y = radius * np.sin(phi) * np.sin(theta) + x = radius * np.sin(phi) * np.cos(theta) + pts = np.stack([z, y, x], axis=1) + np.array([50.0, 50.0, 50.0]) + result = fit_vesicle(pts) + rel = result['reliability'] + assert rel['lat_span_deg'] < 179.0 + assert rel['span_ok'] is False + assert rel['is_reliable'] is False def test_noisy_points_fail_rmse_gate(self): '''Noise comparable to the radius pushes relative RMSE over threshold.''' pts = _sphere_points(radius=20.0, n=800, noise=6.0, seed=12) diff --git a/tests/unit/mrc_unit_test.py b/tests/unit/mrc_unit_test.py index 6bac2f2..f59da1b 100644 --- a/tests/unit/mrc_unit_test.py +++ b/tests/unit/mrc_unit_test.py @@ -4,9 +4,12 @@ import pytest # -- Import internal dependencies ------ +from evaluator.utils import mrc as mrcutil from evaluator.utils.mrc import ( + FITTED_MRC_LABEL, labelComponents, readMRCFile, + readMRCLabels, validateMRCFile, writeMRCFile, ) @@ -61,6 +64,24 @@ def test_data_values_preserved(self, tmp_path): original = self._write_simple(p, shape=(2, 3, 4)) read, _ = readMRCFile(p) assert np.allclose(read, original) + def test_anisotropic_voxel_size_warns(self, tmp_path, monkeypatch): + p = tmp_path / "aniso.mrc" + data = np.ones((4, 4, 4), dtype=np.float32) + with mrcfile.new(str(p)) as mrc: + mrc.set_data(data) + mrc.voxel_size = (5.36, 8.0, 5.36) + warnings = [] + monkeypatch.setattr(mrcutil.lg, "warning", lambda msg: warnings.append(msg)) + _, vox_nm = readMRCFile(p) + assert np.isclose(vox_nm, 0.536, rtol=1e-3) + assert any("anisotropic" in w for w in warnings) + def test_isotropic_voxel_size_does_not_warn(self, tmp_path, monkeypatch): + p = tmp_path / "iso.mrc" + self._write_simple(p, vox_a=5.36) + warnings = [] + monkeypatch.setattr(mrcutil.lg, "warning", lambda msg: warnings.append(msg)) + readMRCFile(p) + assert not any("anisotropic" in w for w in warnings) # -- Define write data to MRC file test class TestWriteMRCFile: @@ -85,6 +106,15 @@ def test_creates_file(self, tmp_path): assert not p.exists() writeMRCFile(np.zeros((2, 2, 2), dtype=np.float32), 1.0, p) assert p.exists() + def test_no_fitted_label_by_default(self, tmp_path): + # mrcfile writes a "Created by mrcfile.py" label so only FITTED_MRC_LABEL marker should be absent + p = tmp_path / "plain.mrc" + writeMRCFile(np.zeros((2, 2, 2), dtype=np.float32), 1.0, p) + assert FITTED_MRC_LABEL not in readMRCLabels(p) + def test_label_roundtrips(self, tmp_path): + p = tmp_path / "labelled.mrc" + writeMRCFile(np.zeros((2, 2, 2), dtype=np.float32), 1.0, p, label=FITTED_MRC_LABEL) + assert readMRCLabels(p) == [FITTED_MRC_LABEL] # -- Define label components test ------ class TestLabelComponents: diff --git a/tests/unit/plot_unit_test.py b/tests/unit/plot_unit_test.py index 12f1a1b..0e66678 100644 --- a/tests/unit/plot_unit_test.py +++ b/tests/unit/plot_unit_test.py @@ -5,7 +5,7 @@ # ==================== # Import external dependencies # ==================== -import pytest, subprocess, tomllib +import pandas as pd, plotly.graph_objects as go, pytest, tomllib from pathlib import Path from unittest.mock import MagicMock, patch @@ -13,16 +13,14 @@ # Import evaluator plot functions # ==================== from evaluator.commands.plot.plot import run_plot -from evaluator.commands.plot.utils.dispatch import ( - resolve_rscript, - dispatch, - RscriptNotFoundError, - RscriptError, -) +from evaluator.commands.plot import sections as sectionutil +from evaluator.commands.plot.utils import data as datautil +from evaluator.commands.plot.utils import export as exportutil from evaluator.commands.plot.utils.input import ( PlotRun, - resolve_plot_inputs, available_sections, + parse_agreement, + resolve_plot_inputs, ) # ==================== @@ -50,6 +48,13 @@ def _make_analyse_csv(tmp_path): ) return path +@pytest.fixture +def fake_sections(): + names = ['distributions', 'qc', 'scatter', 'correlations', 'concordance', 'reliability', 'agreement'] + mocks = {n: MagicMock() for n in names} + with patch.dict('evaluator.commands.plot.sections.SECTIONS', mocks, clear=True), patch('evaluator.commands.plot.plot.exportutil.check_kaleido'), patch('evaluator.commands.plot.plot.sectionutil.load'): + yield mocks + # ==================== # Define tests for plot/utils/input.py # ==================== @@ -130,7 +135,7 @@ def test_missing_sample_id_column_raises(self, tmp_path): # == Define tests for available_sections == def test_analyse_only_gives_core_sections(self): runs = [PlotRun(sample_id='s', analyse_path=Path('a.csv'), model_path=None)] - assert available_sections(runs, False) == ['distributions', 'qc', 'scatter'] + assert available_sections(runs, False) == ['distributions', 'qc', 'scatter', 'correlations'] def test_model_only_gives_no_sections(self): '''No section currently runs on model output alone''' @@ -141,6 +146,8 @@ def test_analyse_and_model_adds_concordance(self): runs = [PlotRun(sample_id='s', analyse_path=Path('a.csv'), model_path=Path('m.json'))] sections = available_sections(runs, False) assert 'concordance' in sections + assert 'reliability' in sections + assert 'agreement' in sections def test_single_group_excludes_compare(self): runs = [ @@ -160,120 +167,228 @@ def test_single_run_mode_never_adds_compare(self): '''compare requires multi_run even if (somehow) groups differ''' runs = [PlotRun(sample_id='s', analyse_path=Path('a.csv'), model_path=None, group='g')] assert 'compare' not in available_sections(runs, False) - -# ==================== -# Define tests for plot/utils/dispatch.py -# ==================== -class TestPlotUtilsDispatch: - # == Define tests for resolve_rscript == - def test_explicit_path_that_exists(self, tmp_path): - fake_bin = tmp_path / 'Rscript' - fake_bin.touch() - assert resolve_rscript(fake_bin) == fake_bin - - def test_explicit_path_missing_raises(self, tmp_path): - missing = tmp_path / 'does_not_exist' - with pytest.raises(RscriptNotFoundError): - resolve_rscript(missing) - - def test_falls_back_to_path_lookup(self, monkeypatch): - monkeypatch.setattr('shutil.which', lambda name: '/usr/bin/Rscript') - assert resolve_rscript(None) == Path('/usr/bin/Rscript') - - def test_not_on_path_raises(self, monkeypatch): - monkeypatch.setattr('shutil.which', lambda name: None) - with pytest.raises(RscriptNotFoundError): - resolve_rscript(None) - - # == Define tests for dispatch == - def test_builds_expected_command(self): - with patch('subprocess.run') as mock_run: - mock_run.return_value = MagicMock(returncode=0, stdout='', stderr='') - dispatch(Path('/usr/bin/Rscript'), 'qc.R', ['/out', '/analyse.csv', 'single', 0.05]) - cmd = mock_run.call_args.args[0] - assert cmd[0] == '/usr/bin/Rscript' - assert cmd[1].endswith('r/qc.R') - assert cmd[2:] == ['/out', '/analyse.csv', 'single', '0.05'] - - def test_nonzero_exit_raises_rscript_error_with_stderr(self): - with patch('subprocess.run') as mock_run: - mock_run.return_value = MagicMock(returncode=1, stdout='', stderr='boom') - with pytest.raises(RscriptError, match='boom'): - dispatch(Path('/usr/bin/Rscript'), 'qc.R', []) - - def test_zero_exit_does_not_raise(self): - with patch('subprocess.run') as mock_run: - mock_run.return_value = MagicMock(returncode=0, stdout='', stderr='') - dispatch(Path('/usr/bin/Rscript'), 'qc.R', []) + + def test_parse_agreement_repeatable_and_space_separated(self): + assert parse_agreement(['a,b', 'c,d e,f']) == [('a', 'b'), ('c', 'd'), ('e', 'f')] + + @pytest.mark.parametrize('bad', ['a', 'a,b,c', ',b', 'a,']) + def test_parse_agreement_rejects_malformed_pairs(self, bad): + with pytest.raises(ValueError): + parse_agreement([bad]) # ==================== # Define tests for plot/plot.py # ==================== class TestPlot: # == Define tests for run_plot section gating == - def test_requested_but_unrunnable_section_is_skipped(self, tmp_path): - '''--section concordance with only --analyse given is not runnable and should be skipped''' + def test_requested_but_unrunnable_section_is_skipped(self, tmp_path, fake_sections): analyse = _make_analyse_csv(tmp_path) - with patch('evaluator.commands.plot.plot.resolve_rscript', return_value='Rscript'), \ - patch('evaluator.commands.plot.plot.dispatch') as mock_dispatch: - run_plot(analyse, None, tmp_path, sections=['concordance'], all_sections=False, overwrite=False, rscript=None) - mock_dispatch.assert_not_called() + run_plot(analyse, None, tmp_path, sections=['concordance'], all_sections=False, overwrite=False) + fake_sections['concordance'].assert_not_called() - def test_all_sections_runs_every_runnable_section(self, tmp_path): + def test_all_sections_runs_every_runnable_section(self, tmp_path, fake_sections): analyse = _make_analyse_csv(tmp_path) - with patch('evaluator.commands.plot.plot.resolve_rscript', return_value='Rscript'), \ - patch('evaluator.commands.plot.plot.dispatch') as mock_dispatch: - run_plot(analyse, None, tmp_path, sections=None, all_sections=True, overwrite=False, rscript=None) - called_scripts = {call.args[1] for call in mock_dispatch.call_args_list} - assert called_scripts == {'distributions.R', 'qc.R', 'scatter.R'} + run_plot(analyse, None, tmp_path, sections=None, all_sections=True, overwrite=False) + called = {n for n, m in fake_sections.items() if m.called} + assert called == {'distributions', 'qc', 'scatter', 'correlations'} - def test_qc_passes_fill_ratio_threshold_as_fourth_arg(self, tmp_path): + def test_include_flagged_passed_to_loader(self, tmp_path, fake_sections): analyse = _make_analyse_csv(tmp_path) - with patch('evaluator.commands.plot.plot.resolve_rscript', return_value='Rscript'), \ - patch('evaluator.commands.plot.plot.dispatch') as mock_dispatch: - run_plot(analyse, None, tmp_path, sections=['qc'], all_sections=False, overwrite=False, rscript=None) - args = mock_dispatch.call_args.args[2] - assert args[3] == 0.05 # PlotConfig default fill_ratio_flag_threshold - - def test_section_failure_is_not_fatal(self, tmp_path): - '''One section raising RscriptError must not prevent other sections from running''' + with patch('evaluator.commands.plot.plot.sectionutil.load') as mock_load: + run_plot(analyse, None, tmp_path, sections=['qc'], all_sections=False, overwrite=False, include_flagged=True) + assert mock_load.call_args.kwargs['include_flagged'] is True + + def test_section_failure_is_not_fatal(self, tmp_path, fake_sections): analyse = _make_analyse_csv(tmp_path) - with patch('evaluator.commands.plot.plot.resolve_rscript', return_value='Rscript'), \ - patch('evaluator.commands.plot.plot.dispatch', side_effect=RscriptError('boom')): - run_plot(analyse, None, tmp_path, sections=['distributions', 'qc'], all_sections=False, overwrite=False, rscript=None) # must not raise + fake_sections['distributions'].side_effect = ValueError('boom') + run_plot(analyse, None, tmp_path, sections=['distributions', 'qc'], all_sections=False, overwrite=False) # must not raise + fake_sections['qc'].assert_called_once() + + # == Define tests for multi-sample output == + def test_multi_sample_writes_subdirectory_per_sample(self, tmp_path, fake_sections): + a1, a2 = _make_analyse_csv(tmp_path), _write(tmp_path / 'a2.csv', 'x') + sheet = _sheet(tmp_path / 'sheet.tsv', [{'sample_id': 's1', 'path': a1}, {'sample_id': 's2', 'path': a2}]) + run_plot(sheet, None, tmp_path, sections=['qc'], all_sections=False, overwrite=False) + assert fake_sections['qc'].call_count == 2 + qc_dir = tmp_path / 'evaluator' / 'plot' / 'qc' + assert (qc_dir / 's1').is_dir() and (qc_dir / 's2').is_dir() # == Define tests for --overwrite option == - def test_existing_section_dir_skipped_without_overwrite(self, tmp_path): + def test_existing_section_dir_skipped_without_overwrite(self, tmp_path, fake_sections): analyse = _make_analyse_csv(tmp_path) - with patch('evaluator.commands.plot.plot.resolve_rscript', return_value='Rscript'), \ - patch('evaluator.commands.plot.plot.dispatch') as mock_dispatch: - run_plot(analyse, None, tmp_path, sections=['qc'], all_sections=False, overwrite=False, rscript=None) - mock_dispatch.reset_mock() - run_plot(analyse, None, tmp_path, sections=['qc'], all_sections=False, overwrite=False, rscript=None) - mock_dispatch.assert_not_called() - - def test_overwrite_reruns_existing_section(self, tmp_path): + run_plot(analyse, None, tmp_path, sections=['qc'], all_sections=False, overwrite=False) + fake_sections['qc'].reset_mock() + run_plot(analyse, None, tmp_path, sections=['qc'], all_sections=False, overwrite=False) + fake_sections['qc'].assert_not_called() + + def test_overwrite_reruns_existing_section(self, tmp_path, fake_sections): analyse = _make_analyse_csv(tmp_path) - with patch('evaluator.commands.plot.plot.resolve_rscript', return_value='Rscript'), \ - patch('evaluator.commands.plot.plot.dispatch') as mock_dispatch: - run_plot(analyse, None, tmp_path, sections=['qc'], all_sections=False, overwrite=False, rscript=None) - mock_dispatch.reset_mock() - run_plot(analyse, None, tmp_path, sections=['qc'], all_sections=False, overwrite=True, rscript=None) - mock_dispatch.assert_called_once() - - # == Define tests for _write_index == - def test_index_md_lists_completed_sections(self, tmp_path): + run_plot(analyse, None, tmp_path, sections=['qc'], all_sections=False, overwrite=False) + fake_sections['qc'].reset_mock() + run_plot(analyse, None, tmp_path, sections=['qc'], all_sections=False, overwrite=True) + fake_sections['qc'].assert_called_once() + + # == Define tests for _write_index and params == + def test_index_md_lists_completed_sections(self, tmp_path, fake_sections): analyse = _make_analyse_csv(tmp_path) - with patch('evaluator.commands.plot.plot.resolve_rscript', return_value='Rscript'), \ - patch('evaluator.commands.plot.plot.dispatch'): - run_plot(analyse, None, tmp_path, sections=['qc'], all_sections=False, overwrite=False, rscript=None) - index = (tmp_path / 'evaluator' / 'plot' / 'index.md').read_text() - assert 'qc/' in index + run_plot(analyse, None, tmp_path, sections=['qc'], all_sections=False, overwrite=False) + assert 'qc/' in (tmp_path / 'evaluator' / 'plot' / 'index.md').read_text() - def test_params_toml_written(self, tmp_path): + def test_params_toml_written(self, tmp_path, fake_sections): analyse = _make_analyse_csv(tmp_path) - with patch('evaluator.commands.plot.plot.resolve_rscript', return_value='Rscript'), \ - patch('evaluator.commands.plot.plot.dispatch'): - run_plot(analyse, None, tmp_path, sections=['qc'], all_sections=False, overwrite=False, rscript=None) + run_plot(analyse, None, tmp_path, sections=['qc'], all_sections=False, overwrite=False) params = tomllib.loads((tmp_path / 'evaluator' / 'plot' / 'params.toml').read_text()) - assert params['fill_ratio_flag_threshold'] == 0.05 + assert params['fill_ratio_flag_threshold'] == 0.05 and params['theme'] == 'Okabe-Ito' + + def test_agreement_pairs_add_the_section(self, tmp_path, fake_sections): + analyse = _make_analyse_csv(tmp_path) + model = _write(tmp_path / 'model.json', '{}') + run_plot(analyse, model, tmp_path, sections=['qc'], all_sections=False, overwrite=False, agreement=[('a', 'b')]) + fake_sections['agreement'].assert_called_once() + +# ==================== +# Define tests for plot/utils/data.py +# ==================== +class TestPlotUtilsData: + def test_aggregate_tomograms_takes_median_and_counts_evs(self): + df = pd.DataFrame({ + 'group': ['a'] * 4, 'sample_id': ['s'] * 4, 'tomogram': ['t1', 't1', 't1', 't2'], + 'x': [1.0, 2.0, 100.0, 5.0], + }) + out = datautil.aggregate_tomograms(df, ['x']) + assert out.set_index('tomogram')['x'].to_dict() == {'t1': 2.0, 't2': 5.0} + assert out.set_index('tomogram')['n_evs'].to_dict() == {'t1': 3, 't2': 1} + + def test_tomogram_note_states_tomogram_n(self): + tomo = pd.DataFrame({'sample_id': ['s1', 's2', 's2'], 'tomogram': ['a', 'b', 'c']}) + note = datautil.tomogram_note(tomo) + assert 'n = 3 tomogram(s)' in note and '2 sample(s)' in note + + def test_attach_fitted_diameter_masks_unreliable(self): + analyse = pd.DataFrame({'tomogram': ['t.mrc'] * 3, 'label': [1, 2, 3]}) + model = pd.DataFrame({ + 'source_file': ['/x/t.mrc', '/x/t.mrc'], 'label_id': [1, 2], 'model_diameter_nm': [100.0, 50.0], + 'reliability.is_reliable': [True, False], + }) + out = datautil.attach_fitted_diameter(analyse, model) + assert out['model_diameter_nm'].tolist()[0] == 100.0 + assert out['model_diameter_nm'].isna().tolist() == [False, True, True] # unreliable fit and unmatched row both NaN + + def test_attach_fitted_diameter_keeps_row_count_with_duplicate_model_rows(self): + analyse = pd.DataFrame({'tomogram': ['t.mrc'], 'label': [1]}) + model = pd.DataFrame({'source_file': ['/x/t.mrc', '/y/t.mrc'], 'label_id': [1, 1], 'model_diameter_nm': [10.0, 20.0]}) + assert len(datautil.attach_fitted_diameter(analyse, model)) == 1 + + def test_populations_split_by_qc_and_z_clip(self): + df = pd.DataFrame({ + 'is_vesicle_like': [True, True, False, True], + 'touches_z_boundary': [False, True, False, False], + }) + pops = datautil.populations(df) + assert {k: len(v) for k, v in pops.items()} == {'all': 4, 'vesicle_like': 3, 'vesicle_like_no_z_clip': 2} + + def test_populations_without_qc_column_is_all_only(self): + assert list(datautil.populations(pd.DataFrame({'x': [1, 2]}))) == ['all'] + + def _csv(self, tmp_path, rows): + header = 'tomogram,label,equiv_diameter_nm,major_axis_diameter,minor_axis_diameter,aspect_ratio,eccentricity,membrane_volume,lumen_volume,surface_area,is_enclosed,closure_fill_ratio,is_vesicle_like,qc_flags\n' + body = ''.join(f'tomo1.mrc,{i},80,90,70,1.2,0.5,5,3,2,True,0.4,{vl},{flags}\n' for i, (vl, flags) in enumerate(rows)) + return _write(tmp_path / 'a.csv', header + body) + + def test_flagged_rows_dropped_by_default(self, tmp_path): + path = self._csv(tmp_path, [(True, ''), (False, 'sphere_rmse'), (False, 'sphere_rmse|x')]) + assert len(datautil.load_analyse(path)) == 1 + + def test_include_flagged_keeps_everything(self, tmp_path): + path = self._csv(tmp_path, [(True, ''), (False, 'sphere_rmse')]) + assert len(datautil.load_analyse(path, include_flagged=True)) == 2 + + def test_missing_vesicle_like_column_is_not_filtered(self, tmp_path): + path = _make_analyse_csv(tmp_path) + assert len(datautil.load_analyse(path)) == 1 + + def test_missing_required_column_raises(self, tmp_path): + path = _write(tmp_path / 'a.csv', 'tomogram,label\ntomo1.mrc,1\n') + with pytest.raises(ValueError, match='missing required'): + datautil.load_analyse(path) + + def test_drop_reasons_counted_per_flag(self, tmp_path, capsys): + from evaluator.utils.settings import logger + messages = [] + sink = logger.add(lambda m: messages.append(str(m)), level='INFO') + try: + path = self._csv(tmp_path, [(True, ''), (False, '"sphere_rmse,open_shell"'), (False, 'sphere_rmse')]) + datautil.load_analyse(path) + finally: + logger.remove(sink) + assert any('sphere_rmse=2' in m and 'open_shell=1' in m for m in messages) + + def test_join_matches_on_basename_and_label_and_counts_unmatched(self): + analyse = pd.DataFrame({'tomogram': ['t.mrc', 't.mrc', 't.mrc'], 'label': [1, 2, 3]}) + model = pd.DataFrame({'source_file': ['/x/t.mrc', '/x/t.mrc'], 'label_id': [1, 9], 'radius': [10.0, 5.0]}) + joined, a_only, m_only = datautil.join_analyse_model(analyse, model) + assert list(joined['label']) == [1] and (a_only, m_only) == (2, 1) + + def test_pooling_note_states_true_n(self): + df = pd.DataFrame({'tomogram': ['a', 'a', 'b']}) + note = datautil.pooling_note(df, n_samples=2) + assert 'n = 3 EVs' in note and '2 tomogram(s)' in note and '2 sample(s)' in note + +# ==================== +# Define tests for plot/utils/export.py +# ==================== +def _kaleido_ok() -> bool: + try: + exportutil.check_kaleido() + return True + except RuntimeError: + return False + +class TestPlotUtilsExport: + def test_vectorise_replaces_scattergl(self): + fig = go.Figure(go.Scattergl(x=[1, 2], y=[1, 2], mode='markers')) + assert fig.data[0].type == 'scattergl' + assert exportutil._vectorise(fig).data[0].type == 'scatter' + + def test_write_xlsx_first_sheet_active(self, tmp_path): + from openpyxl import load_workbook + path = tmp_path / 't.xlsx' + exportutil.write_xlsx(path, {'main': pd.DataFrame({'a': [1, 2]}), 'notes': pd.DataFrame({'note': ['x']})}) + wb = load_workbook(path) + assert wb.sheetnames == ['main', 'notes'] and wb.active.title == 'main' + + @pytest.mark.skipif(not _kaleido_ok(), reason='kaleido/Chrome unavailable') + def test_write_figure_svg(self, tmp_path): + path = tmp_path / 'f.svg' + exportutil.write_figure(go.Figure(go.Scatter(x=[1, 2], y=[1, 2])), path) + assert path.read_text().lstrip().startswith(('