multi_file_pipeline.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
The hard case: TWO meshes per timepoint.

Big samples are often exported as several overlapping meshes (e.g. an abaxial and an adaxial
view of the same tissue). That means one physical cell can appear in two meshes under two
different labels, neighbours can sit across the seam, and a cell's parent can live in the
*other* mesh. This example runs on examples/data_multi and shows how Teal handles all three.

The toy tissue: meshes A (file_n=0) and B (file_n=1), sharing a strip of cells.

    T0:  A[1 2 3 4]      B[101 102 103]        4 == 101   (the same physical cell)
    T1:  A[11..15]       B[111..114]          15 == 111
    T2:  A[21..26]       B[121..125]          26 == 121

Four ideas, in order:
    1. overlaps      -- one cell, two labels, ONE node
    2. neighbors     -- spatial adjacency across the seam between meshes
    3. parents       -- one parent file per mesh
    4. crossparents  -- a cell whose parent lives in the other mesh

Run it from anywhere:  python examples/multi_file_pipeline.py
"""

import os
import sys

sys.path.append(os.path.join(os.path.dirname(__file__), ".."))

import teal as tl
from teal import files_by_t

DATA = os.path.join(os.path.dirname(__file__), 'data_multi')


# ============================================================================
# 1. Gather the files
# ============================================================================

# Two meshes per timepoint means TWO files behind each timepoint key. Their order is what
# defines file_n: sorted by filename, so T0_A.ply -> file_n=0 and T0_B.ply -> file_n=1.
# Every other per-mesh file (parents, attributes) follows that same order.
graph_files = files_by_t(f'{DATA}/graphs', 0)      # {0: [T0_A.ply, T0_B.ply], ...}
parent_files = files_by_t(f'{DATA}/parents', 0)    # {1: [T1_A_parents.csv, T1_B_parents.csv], ...}
size_files = files_by_t(f'{DATA}/size', 0)

# overlaps: which cell of mesh A is the same physical cell as which cell of mesh B.
# one file per PAIR of meshes -- with 2 meshes that is 1 file; with 3 it would be 3.
overlaps = files_by_t(f'{DATA}/overlaps', 0)

# neighbors: adjacency the meshes cannot see on their own, because the two cells sit in
# different meshes. Format is `label_a, "label_b label_c ..."`.
neighbors = files_by_t(f'{DATA}/neighbors', 0)

# crossparents: for a cell whose parent is in the OTHER mesh. The filename carries the
# routing:  T<t>_<parent_mesh><child_mesh>
#   T1_01.csv -> at T1, parent is in mesh 0, child is in mesh 1
#   T2_10.csv -> at T2, parent is in mesh 1, child is in mesh 0
cross_parent_files = files_by_t(f'{DATA}/crossparents', 0)


# ============================================================================
# 2. Build
# ============================================================================

"""
Building prints a few warnings. They are worth reading rather than silencing:

  t1 [T1_B_parents.csv]: 1 cells have no value: [114]
  t2 [T2_A_parents.csv]: 1 cells have no value: [25]

That is precisely the cross-parent signal -- those two cells have no parent inside their own
mesh's file, because their parent is in the other mesh. The crossparents files below fill them
in (report() confirms it: 0 unparented cells).

  t2 [<dataframe>]: 1 cells have no value: [26]

This one is benign. 26 is a merged cell, so it has a row in each mesh; the mesh-A row cannot
resolve its parent but the mesh-B row can, and they are the same node. Use teal.verbose(False)
once you know which gaps are expected.
"""

tmlps = tl.make_timelapse(
    graph_files,
    overlaps=overlaps,
    neighbors=neighbors,
    parent_files=parent_files,
    cross_parent_files=cross_parent_files,
    sample_name='multi',
)

print(tmlps)
for t in range(tmlps.time):
    print(f'  t{t}: meshes={sorted(tmlps[t].file_ns)}  cells={len(tmlps[t].ids)}')


# ============================================================================
# 3. Overlaps -- one cell, two labels, ONE node
# ============================================================================

# A cell seen in both meshes is NOT two cells. Teal merges them into a single node that
# remembers both (label, file_n) pairs -- which is also how heat() writes a value back out
# to both meshes later.
print('\ncells that were merged across the two meshes:')
for node, attrs in tmlps.data.nodes(data=True):
    if len(attrs.get('labels', [])) > 1:
        print(f'  node {node}  <-  {attrs["labels"]}')

# The consequence: a merged cell gets measured TWICE, once per mesh. Loading with
# overwrite='add' keeps both values as a list, and a Filter decides what they mean.
tmlps.load_attr('SIZE', size_files, overwrite='add')

merged = next(n for n, d in tmlps.data.nodes(data=True) if len(d.get('labels', [])) > 1)
print(f'\nmerged cell {merged} measured in both meshes -> SIZE = {tmlps.data.nodes[merged]["SIZE"]}')

tmlps.filter_attr('SIZE', 'avg_all')   # reconcile the two measurements
print(f'after filter_attr(avg_all)                    -> SIZE = {tmlps.data.nodes[merged]["SIZE"]}')


# ============================================================================
# 4. Neighbors -- adjacency across the seam
# ============================================================================

# label -> node, so we can talk about cells by their mesh label
node_of = {lab: n for n, d in tmlps.data.nodes(data=True) for (lab, _) in d.get('labels', [])}

# cell 3 lives in mesh A, cell 102 in mesh B. Neither mesh knows they touch: only the
# neighbors file does. Without it the two meshes would be spatially disconnected.
a, b = node_of[3], node_of[102]
touching = tmlps.data.has_edge(a, b) or tmlps.data.has_edge(b, a)
print(f'\nspace edge between mesh-A cell 3 and mesh-B cell 102: {touching}')


# ============================================================================
# 5. Cross-parents -- a parent in the other mesh
# ============================================================================

# Cell 114 only exists in mesh B, but its parent (cell 3) only exists in mesh A. Mesh B's
# own parent file cannot express that, so it goes in crossparents/T1_01.csv instead.
# Same story in reverse for cell 25 (mesh A) whose parent 112 is in mesh B.
print('\ncross-parented cells:')
print(f'  mesh-A cell 3   -> mesh-B cell 114 : {tmlps.data.has_edge(node_of[3], node_of[114])}')
print(f'  mesh-B cell 112 -> mesh-A cell 25  : {tmlps.data.has_edge(node_of[112], node_of[25])}')

# If a cross-parent file were misnumbered, the parent lookup would fall back to a raw label
# and invent a cell that never existed. That is exactly what phantom_cells() catches, and
# why report() should always show 0 of them.
print('\nreport:')
tmlps.report(attrs=['SIZE'])