Refactor ReaDDy processing and all visualization

.gitignore

@@ -123,3 +123,4 @@ ENV/
 *.cym
 *.simularium
 **/analysis_outputs/**
+*.h5

docs/conf.py

@@ -58,7 +58,7 @@
 
 # List of modules to be mocked up. Useful when some external dependencies are
 # not met at build time and break the building process.
-autodoc_mock_imports = []
+autodoc_mock_imports = ["readdy"]
 
 # Controls how to represent typehints.
 autodoc_typehints = "signature"

docs/index.rst

@@ -16,6 +16,7 @@
 
    Simulations <simulation>
    Analysis <analysis>
+   Visualization <visualization>
 
 .. autosummary::
    :toctree: _summary
@@ -25,6 +26,7 @@
 
    subcell_pipeline.simulation
    subcell_pipeline.analysis
+   subcell_pipeline.visualization
 
 .. toctree::
    :hidden:

docs/visualization.rst

@@ -0,0 +1,5 @@
+Visualization workflow notebooks
+================================
+
+.. include:: ../subcell_pipeline/visualization/README.md
+   :parser: myst_parser.sphinx_

pyproject.toml

@@ -27,6 +27,7 @@ dependencies = [
     "setuptools>=70.0.0",
     "io-collection>=0.10.2",
     "python-dotenv>=1.0.1",
+    "imageio",
 ]
 
 [project.urls]
@@ -95,6 +96,7 @@ extend-ignore = [
     "D100", # Missing docstring in public module
     "D101", # Missing docstring in public class
     "D103", # Missing docstring in public function
+    "D105", # Missing docstring in magic method
     "D107", # Missing docstring in __init__
     "D202", # Blank lines between the function body and the function docstring
     "D203", # 1 blank line required before class docstring

subcell_pipeline/analysis/compression_metrics/README.md

@@ -2,6 +2,6 @@
 
 ## Metrics for comparing traces of compressed fibers
 
-Analysis combines compression simulations from Cytosim and Readdy and calculates various metrics to compare the compressed fibers.
+Analysis combines compression simulations from Cytosim and Readdy and calculates various compression metrics metrics to compare fibers.
 
-- **Calculate compression metrics** ([source](https://github.com/simularium/subcell-pipeline/blob/main/subcell_pipeline/analysis/compression_metrics/_compare_compression_metrics.py) | [notebook](https://simularium.github.io/subcell-pipeline/_notebooks/analysis/compression_metrics/_compare_compression_metrics.html))
+- **Compare compression metrics between simulators** ([source](https://github.com/simularium/subcell-pipeline/blob/main/subcell_pipeline/analysis/compression_metrics/_compare_compression_metrics.py) | [notebook](https://simularium.github.io/subcell-pipeline/_notebooks/analysis/compression_metrics/_compare_compression_metrics.html))

subcell_pipeline/analysis/compression_metrics/_compare_compression_metrics.py

@@ -1,5 +1,5 @@
 # %% [markdown]
-# # Compare metrics across simulators
+# # Compare compression metrics between simulators
 
 # %% [markdown]
 """
@@ -39,12 +39,14 @@
 
 Defines the `COMPRESSION_VELOCITY` simulation series, which compresses a single
 500 nm actin fiber at four different velocities (4.7, 15, 47, and 150 μm/s) with
-five replicates each (random seeds 1, 2, 3, 4, and 5).
+five replicates each and the baseline `NO_COMPRESSION` simulation series, which
+simulates a single actin fiber with a free barbed end across five replicates.
 """
 
 # %%
 # Name of the simulation series
-series_name: str = "COMPRESSION_VELOCITY"
+compression_series_name: str = "COMPRESSION_VELOCITY"
+no_compression_series_name: str = "NO_COMPRESSION"
 
 # S3 bucket Cytosim for input and output files
 cytosim_bucket: str = "s3://cytosim-working-bucket"
@@ -89,39 +91,61 @@
 """
 
 # %%
-cytosim_metrics = get_compression_metric_data(
+cytosim_metrics_compression = get_compression_metric_data(
     bucket=cytosim_bucket,
-    series_name=series_name,
+    series_name=compression_series_name,
     condition_keys=condition_keys,
     random_seeds=random_seeds,
     metrics=metrics,
     recalculate=recalculate,
 )
-cytosim_metrics["simulator"] = "cytosim"
+cytosim_metrics_compression["simulator"] = "cytosim"
+
+# %%
+cytosim_metrics_no_compression = get_compression_metric_data(
+    bucket=cytosim_bucket,
+    series_name=no_compression_series_name,
+    condition_keys=[""],
+    random_seeds=random_seeds,
+    metrics=metrics,
+    recalculate=recalculate,
+)
+cytosim_metrics_no_compression["simulator"] = "cytosim"
 
 # %% [markdown]
 """
 ## Calculate metrics for ReaDDy data
 """
 
 # %%
-readdy_metrics = get_compression_metric_data(
+readdy_metrics_compression = get_compression_metric_data(
     bucket=readdy_bucket,
-    series_name=series_name,
+    series_name=f"ACTIN_{compression_series_name}",
     condition_keys=condition_keys,
     random_seeds=random_seeds,
     metrics=metrics,
     recalculate=recalculate,
 )
-readdy_metrics["simulator"] = "readdy"
+readdy_metrics_compression["simulator"] = "readdy"
+
+# %%
+readdy_metrics_no_compression = get_compression_metric_data(
+    bucket=readdy_bucket,
+    series_name=f"ACTIN_{no_compression_series_name}",
+    condition_keys=[""],
+    random_seeds=random_seeds,
+    metrics=metrics,
+    recalculate=recalculate,
+)
+readdy_metrics_no_compression["simulator"] = "readdy"
 
 # %% [markdown]
 """
 ## Combine metrics from both simulators
 """
 
 # %%
-combined_metrics = pd.concat([cytosim_metrics, readdy_metrics])
+combined_metrics = pd.concat([cytosim_metrics_compression, readdy_metrics_compression])
 combined_metrics["repeat"] = combined_metrics["seed"] - 1
 combined_metrics["velocity"] = combined_metrics["key"].astype("int") / 10
 
@@ -135,7 +159,6 @@
     combined_metrics, str(save_location), "actin_compression_combined_metrics.csv"
 )
 
-
 # %% [markdown]
 """
 ## Plot metrics vs time

subcell_pipeline/analysis/compression_metrics/compression_analysis.py

@@ -131,7 +131,10 @@ def calculate_compression_metrics(
                 polymer_trace=polymer_trace, **options
             )
 
-    return df_metrics.reset_index().rename(columns={"index": "time"})
+    metrics = df_metrics.reset_index().rename(columns={"index": "time"})
+    metrics["normalized_time"] = metrics["time"] / metrics["time"].max()
+
+    return metrics
 
 
 def save_compression_metrics(

subcell_pipeline/analysis/compression_metrics/compression_metric.py

@@ -61,6 +61,58 @@ def label(self: Enum) -> str:
         }
         return labels.get(self.value, "")
 
+    def description(self: Enum) -> str:
+        """
+        Return the description for the compression metric.
+
+        Parameters
+        ----------
+        self
+            the CompressionMetric object
+
+        Returns
+        -------
+        :
+            The description (and units) for the compression metric.
+        """
+        units = {
+            CompressionMetric.NON_COPLANARITY.value: "3rd component variance from PCA",
+            CompressionMetric.PEAK_ASYMMETRY.value: "normalized peak distance",
+            CompressionMetric.SUM_BENDING_ENERGY.value: "sum of bending energy",
+            CompressionMetric.AVERAGE_PERP_DISTANCE.value: "distance (nm)",
+            CompressionMetric.TOTAL_FIBER_TWIST.value: "total fiber twist",
+            CompressionMetric.CALC_BENDING_ENERGY.value: "energy",
+            CompressionMetric.CONTOUR_LENGTH.value: "filament contour length (nm)",
+            CompressionMetric.COMPRESSION_RATIO.value: "compression ratio",
+        }
+        return units.get(self.value, "")
+
+    def bounds(self: Enum) -> tuple[float, float]:
+        """
+        Return the default bounds for the compression metric.
+
+        Parameters
+        ----------
+        self
+            the CompressionMetric object
+
+        Returns
+        -------
+        :
+            The default bounds for the compression metric.
+        """
+        bounds = {
+            CompressionMetric.NON_COPLANARITY.value: (0, 0.03),
+            CompressionMetric.PEAK_ASYMMETRY.value: (0, 0.5),
+            CompressionMetric.SUM_BENDING_ENERGY.value: (0, 0),  # TODO
+            CompressionMetric.AVERAGE_PERP_DISTANCE.value: (0, 85.0),
+            CompressionMetric.TOTAL_FIBER_TWIST.value: (0, 0),  # TODO
+            CompressionMetric.CALC_BENDING_ENERGY.value: (0, 10),
+            CompressionMetric.CONTOUR_LENGTH.value: (480, 505),
+            CompressionMetric.COMPRESSION_RATIO.value: (0, 1),  # TODO
+        }
+        return bounds.get(self.value, (0, 0))
+
     def calculate_metric(
         self, polymer_trace: np.ndarray, **options: dict[str, Any]
     ) -> Union[float, np.floating[Any]]:

subcell_pipeline/analysis/compression_metrics/polymer_trace.py

@@ -1,6 +1,6 @@
 """Methods to calculate metrics from polymer trace data."""
 
-from typing import Any, Dict, Tuple, Union
+from typing import Any, Dict, Tuple
 
 import numpy as np
 from sklearn.decomposition import PCA
@@ -61,7 +61,7 @@ def get_end_to_end_axis_distances_and_projections(
 def get_average_distance_from_end_to_end_axis(
     polymer_trace: np.ndarray,
     **options: Dict[str, Any],
-) -> Union[float, np.floating[Any]]:
+) -> float:
     """
     Calculate the average perpendicular distance of polymer trace points from
     the end-to-end axis.
@@ -152,7 +152,7 @@ def get_pca_polymer_trace_projection(
 def get_contour_length_from_trace(
     polymer_trace: np.ndarray,
     **options: Dict[str, Any],
-) -> Union[float, np.floating[Any]]:
+) -> float:
     """
     Calculate the sum of inter-monomer distances in the trace.
 
@@ -172,13 +172,13 @@ def get_contour_length_from_trace(
     total_distance = np.float_(0)
     for i in range(len(polymer_trace) - 1):
         total_distance += np.linalg.norm(polymer_trace[i] - polymer_trace[i + 1])
-    return total_distance
+    return total_distance.item()
 
 
 def get_bending_energy_from_trace(
     polymer_trace: np.ndarray,
     **options: Dict[str, Any],
-) -> Union[float, np.floating[Any]]:
+) -> float:
     """
     Calculate the bending energy per monomer of a polymer trace.
 
@@ -207,6 +207,14 @@ def get_bending_energy_from_trace(
         vec1 = polymer_trace[ind + 1] - polymer_trace[ind]
         vec2 = polymer_trace[ind + 2] - polymer_trace[ind + 1]
 
+        if np.isclose(np.linalg.norm(vec1), 0.0) or np.isclose(
+            np.linalg.norm(vec2), 0.0
+        ):
+            # TODO handle this differently?
+            cos_angle[ind] = 0.0
+            print("Warning: zero vector in bending energy calculation.")
+            continue
+
         cos_angle[ind] = (
             np.dot(vec1, vec2) / np.linalg.norm(vec1) / np.linalg.norm(vec2)
         )
@@ -215,7 +223,7 @@ def get_bending_energy_from_trace(
     # the type checker is unable to infer its type
     energy = bending_constant * (1 - np.nanmean(cos_angle))
 
-    return energy
+    return energy.item()
 
 
 def get_total_fiber_twist(
@@ -454,7 +462,7 @@ def get_sum_bending_energy(
 def get_compression_ratio(
     polymer_trace: np.ndarray,
     **options: Dict[str, Any],
-) -> Union[float, np.floating[Any]]:
+) -> float:
     """
     Calculate the compression ratio of a polymer trace.
 
@@ -474,5 +482,5 @@ def get_compression_ratio(
     :
         The compression ratio of the polymer trace.
     """
-    end_to_end_axis_length = np.linalg.norm(polymer_trace[-1] - polymer_trace[0])
+    end_to_end_axis_length = np.linalg.norm(polymer_trace[-1] - polymer_trace[0]).item()
     return 1 - end_to_end_axis_length / get_contour_length_from_trace(polymer_trace)

subcell_pipeline/analysis/dimensionality_reduction/_run_pacmap_on_compression_simulations.py

@@ -67,7 +67,9 @@
 """
 
 # %%
-readdy_data = get_merged_data(readdy_bucket, series_name, condition_keys, random_seeds)
+readdy_data = get_merged_data(
+    readdy_bucket, f"ACTIN_{series_name}", condition_keys, random_seeds
+)
 readdy_data["simulator"] = "readdy"
 
 # %%