Reuse results and compare analysis settings

# Open count stores and run CyteArc analyses.
import cytearc

# Keep routine logs and progress bars out of the results.
cytearc.configure_output(level="WARNING", progress=False)

# Download the prepared example store.
dataset = cytearc.cytebase.connect("cytearc_docs").download_dataset(
    "tabula_sapiens_100k",
    destination="cytearc_datasets",
    zarr=True,
)

# Open the count store for this analysis.
ds = cytearc.DataStore(f"{dataset}/data.zarr", nthreads=4, min_features_per_cell=-1)

# Open the saved analysis used to check reuse.
baseline_run = ds.pipeline.open(label="docs_default")

# Keep the cells used by the saved analysis.
cell_selection = baseline_run["analysis_cell_selection"]

# Keep the variable genes used by the saved analysis.
hvg_ref = baseline_run["highly_variable_features"]

# Inspect the opened store's cells and features.
ds

# Reopen the saved normalized expression.
normalized = baseline_run["normalized"]

# Reopen the saved PCA result.
pca = baseline_run["pca"]

# Keep the saved nearest-neighbor index.
ann = baseline_run["ann_index"]

# Keep the saved 21-neighbor result.
neighbors_k21 = baseline_run["neighbors"]

# Keep the baseline connectivity graph.
graph_k21 = baseline_run["connectivity_map"]

# Query the existing index with 31 neighbors per cell.
neighbors_k31 = ds.graph.neighbors(ann, k=31)

# Build connectivity for the new neighbor count.
graph_k31 = ds.graph.connectivity(neighbors_k31)

# Inspect which saved index the new neighbor query used.
query_status = ds.artifacts.inspect(neighbors_k31)

# Check which results were reused and which changed.
reuse_checks = {
    "normalization reused": ds.features.normalize(cell_selection, hvg_ref) == normalized,
    "PCA reused": ds.reduction.pca(normalized, dims=30) == pca,
    "ANN index reused": query_status.input_ref("ann_index") == ann,
    "neighbors changed": neighbors_k31 != neighbors_k21,
    "graph changed": graph_k31 != graph_k21,
}

# Stop if any result differs from the reuse behavior described here.
assert all(reuse_checks.values()), reuse_checks

# Display the verified comparison.
reuse_checks

# Fit a separate reduction with 40 principal components.
pca_dims40 = ds.reduction.pca(normalized, dims=40)

# Build a search index from the 40-component reduction.
ann_dims40 = ds.graph.ann_index(pca_dims40)

# Find 21 neighbors using the new reduction.
neighbors_dims40 = ds.graph.neighbors(ann_dims40, k=21)

# Build the connectivity graph for this reduction.
graph_dims40 = ds.graph.connectivity(neighbors_dims40)

# Check which results were reused and which changed.
branch_checks = {
    "PCA changed": pca_dims40 != pca,
    "ANN index changed": ann_dims40 != ann,
    "neighbors changed": neighbors_dims40 != neighbors_k21,
    "graph changed": graph_dims40 != graph_k21,
}

# Stop if the new branch unexpectedly shares one of these results.
assert all(branch_checks.values()), branch_checks

# Display the verified comparison.
branch_checks

# Revisit the original references after building the 40-component branch.
{
    "baseline PCA": pca,
    "baseline neighbors": neighbors_k21,
    "baseline graph": graph_k21,
}

# Trace the shared inputs and differences between all three graphs.
lineage = ds.artifacts.lineage(
    {
        "k21 graph": graph_k21,
        "k31 graph": graph_k31,
        "dims40 graph": graph_dims40,
    }
)

# Display the shared inputs and diverging analysis branches.
lineage

# Export the displayed lineage report as Markdown.
lineage_markdown = lineage.to_markdown()

# Preview the opening lines of the exported report.
lineage_markdown.splitlines()[:12]

# Inspect the saved chain from normalization to connectivity.
{
    "normalization": normalized,
    "PCA": pca,
    "neighbor index": ann,
    "neighbors": neighbors_k21,
    "graph": graph_k21,
}

# Find the exact saved graph from its neighbor input.
[reopened_graph] = ds.artifacts.list(
    from_assay="RNA",
    kind="connectivity_map",
    operation="build_connectivity_map",
    inputs={"neighbors": neighbors_k21},
    complete_only=True,
)

# Check that the catalog returned the baseline graph.
assert reopened_graph == graph_k21

# Inspect the saved graph's operation, parameters, and inputs.
status = ds.artifacts.inspect(reopened_graph)

# Inspect the graph operation and the exact inputs it records.
{
    "operation": status.operation,
    "parameters": status.parameters,
    "inputs": status.inputs,
    "complete": status.complete,
}