.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples\06_matching\plot_06_downstream_tables.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_06_matching_plot_06_downstream_tables.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_06_matching_plot_06_downstream_tables.py:


What matching changes in a cohort feature table
===============================================

**Background.** Habitat features are usually analysed as a table with one
row per patient and one column per habitat. That only works if a column
refers to the same kind of habitat in every patient.

**Purpose.** You build one cohort table from five demo subjects with raw
and with prototype-matched ids, compare the spread of one column, and see
how ``max_distance`` can leave poorly matching habitats out of the shared
columns.

**Key terms.**

* **prototype** -- see
  :doc:`/auto_examples/06_matching/plot_03_prototype_steps`.
* **volume fraction** -- the share of ROI voxels that belong to one
  habitat.
* **max_distance** -- an optional cutoff: a habitat farther than this from
  every free prototype keeps a subject-local id instead of a shared one.

A cohort habitat table has one column per habitat id: ``H1_volume``,
``H2_mean_enhancement``, ... A column is only one variable if ``H2``
means the same habitat in every row. With per-subject (``one_step``)
fits it does not, until the ids are matched.

This page builds the same table twice -- raw per-subject ids, then
prototype-matched ids -- and shows the spread of one column across
subjects. It ends with ``max_distance``: how a habitat that resembles no
prototype can be left out of the shared columns instead of being forced
into one.

.. GENERATED FROM PYTHON SOURCE LINES 36-39

Per-subject habitats for five subjects
--------------------------------------
sphinx_gallery_thumbnail_number = 1

.. GENERATED FROM PYTHON SOURCE LINES 39-105

.. code-block:: Python

    from pathlib import Path

    import matplotlib.pyplot as plt
    import numpy as np
    import pandas as pd

    from habit.contracts import Cohort, cohort_from_directory
    from habit.datasets import fetch_demo
    from habit.habitat_model import KMeansHabitatModelFitter
    from habit.kernels import habitat_volume_fractions
    from habit.pipeline import voxel_units
    from habit.precision import align_habitat_maps_to_prototypes
    from habit.viz import plot_prototype_matching
    from habit.voxel_features import ExpressionVoxelFeatures

    # Change DATA / MODALITIES / ROI and the expressions to your own layout.
    DATA = fetch_demo()
    MODALITIES = ("pre_contrast", "LAP", "PVP")
    ROI = "LAP"
    cohort = cohort_from_directory(DATA, modalities=MODALITIES, roi=ROI)
    extractor = ExpressionVoxelFeatures(
        features={
            "rel_enh_lap": "(LAP - pre_contrast) / (pre_contrast + eps)",
            "rel_enh_pvp": "(PVP - pre_contrast) / (pre_contrast + eps)",
        },
        roi=ROI,
    )

    maps, models = [], []
    for subject in cohort:
        units = voxel_units(extractor(subject))
        fitter = KMeansHabitatModelFitter(
            min_habitats=2, max_habitats=5, validation="silhouette", n_init=3
        )
        fitter.set_random_state(0)
        model = fitter.fit([units], cohort=Cohort([subject], name=subject.subject_id))
        maps.append(model.assigner()(units))
        models.append(model)
    # Shared names for all five subjects, from their fitted centroids.
    matched = align_habitat_maps_to_prototypes(maps, models=models)


    def cohort_table(habitat_maps, centroid_rows) -> pd.DataFrame:
        """One row per subject: volume fraction and mean arterial enhancement per id.

        Args:
            habitat_maps: One HabitatMap per subject (raw or matched ids).
            centroid_rows: Per subject, ``{habitat_id: centroid row}`` in the
                same ids as the map.

        Returns:
            pd.DataFrame: Columns ``H{k}_volume`` and ``H{k}_enh_lap``; a
            habitat a subject does not have is volume 0 and enhancement NaN.
        """
        ids = sorted({k for rows in centroid_rows for k in rows})
        records = []
        for habitat_map, rows in zip(habitat_maps, centroid_rows):
            fractions = habitat_volume_fractions(habitat_map.label_array, ids)
            record = {"subject": habitat_map.subject_id}
            for k in ids:
                record[f"H{k}_volume"] = round(fractions[k], 3)
                record[f"H{k}_enh_lap"] = round(float(rows[k][0]), 3) if k in rows else np.nan
            records.append(record)
        return pd.DataFrame(records).set_index("subject")






.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    HABIT demo data (cached)
    DATA (preprocessed root): C:\Users\dongm\.habit_data\demo-data-v1\preprocessed

    On-disk inventory of this folder:
      subjects (5): subj001, subj002, subj003, subj004, subj005
      image series: LAP, PVP, delay_3min, pre_contrast
      mask keys:    LAP, PVP, delay_3min, pre_contrast
      example image: images/subj001/delay_3min/WATER__BH_Ax_LAVA_Flex_3min_Series0012.nrrd
      example mask:  masks/subj001/delay_3min/WATER__BH_Ax_LAVA_Flex_10min_Series0017_mask.nrrd

    Your own data must use the same folder tree (change IDs / series names):

      DATA/
        images/<subject_id>/<modality>/<one image file>
        masks/<subject_id>/<roi>/<one mask file>

    Then load it with the same call the demos use:

      cohort = cohort_from_directory(DATA, modalities=("LAP",), roi="LAP")

    Swap DATA / modalities / roi to match your tree. Mask key is often the
    same as one image series (here LAP).




.. GENERATED FROM PYTHON SOURCE LINES 106-110

The same table, raw ids vs matched ids
--------------------------------------
Raw: habitat ``k`` is row ``k`` of each subject's own centroids.
Matched: the assignment table says which prototype id each habitat got.

.. GENERATED FROM PYTHON SOURCE LINES 110-124

.. code-block:: Python

    raw_rows = [{k + 1: row for k, row in enumerate(np.asarray(m.centroids))} for m in models]
    matched_rows = []
    for subject_id, model in zip(matched.assignments.subject_id.unique(), models):
        part = matched.assignments[matched.assignments.subject_id == subject_id]
        centroids = np.asarray(model.centroids)
        matched_rows.append(
            {int(p): centroids[int(h) - 1] for h, p in zip(part.habitat_id, part.prototype_id)}
        )

    raw_table = cohort_table(maps, raw_rows)
    matched_table = cohort_table(matched.habitat_maps, matched_rows)
    print("raw ids\n", raw_table.to_string(), "\n")
    print("matched ids\n", matched_table.to_string())





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    raw ids
              H1_volume  H1_enh_lap  H2_volume  H2_enh_lap  H3_volume  H3_enh_lap
    subject                                                                     
    subj001      0.206       1.484      0.406       0.980      0.387       1.616
    subj002      0.609       2.195      0.391       0.858      0.000         NaN
    subj003      0.383       0.618      0.617       0.297      0.000         NaN
    subj004      0.508       0.836      0.492       1.425      0.000         NaN
    subj005      0.342       3.211      0.658       2.040      0.000         NaN 

    matched ids
              H1_volume  H1_enh_lap  H2_volume  H2_enh_lap  H3_volume  H3_enh_lap
    subject                                                                     
    subj001      0.406       0.980      0.387       1.616      0.206       1.484
    subj002      0.391       0.858      0.609       2.195      0.000         NaN
    subj003      0.617       0.297      0.383       0.618      0.000         NaN
    subj004      0.508       0.836      0.492       1.425      0.000         NaN
    subj005      0.000         NaN      0.658       2.040      0.342       3.211




.. GENERATED FROM PYTHON SOURCE LINES 125-131

Spread of the enhancement columns across subjects. With raw ids a column
mixes weakly and strongly enhancing habitats. After matching a column
gathers the habitats closest to one prototype; what spread remains is
between-subject variation plus habitats that resemble no prototype well
but still had to take a name (one-to-one matching never merges two
habitats of a subject). The last section shows how to find those.

.. GENERATED FROM PYTHON SOURCE LINES 131-157

.. code-block:: Python

    Path("out").mkdir(exist_ok=True)
    fig, axes = plt.subplots(1, 2, figsize=(8.4, 3.6), sharey=True)
    for ax, (label, table) in zip(axes, (("raw ids", raw_table), ("matched ids", matched_table))):
        columns = [c for c in table.columns if c.endswith("_enh_lap")]
        for x, column in enumerate(columns):
            values = table[column].dropna().to_numpy()
            ax.scatter(np.full(values.size, x), values, s=36, color="#0072B2", alpha=0.8)
            ax.plot([x - 0.25, x + 0.25], [values.mean()] * 2, color="black", linewidth=1.5)
        ax.set_xticks(range(len(columns)))
        ax.set_xticklabels([c.split("_")[0] for c in columns])
        ax.set_title(f"{label}: column spread")
        ax.grid(True, axis="y", alpha=0.25)
    axes[0].set_ylabel("mean arterial enhancement (per subject)")
    fig.tight_layout()
    fig.savefig("out/downstream_column_spread.png", dpi=150, bbox_inches="tight")
    plt.show()

    spread = pd.DataFrame(
        {
            "raw_ids_sd": raw_table.filter(like="_enh_lap").std(),
            "matched_ids_sd": matched_table.filter(like="_enh_lap").std(),
        }
    ).round(3)
    print(spread)
    spread




.. image-sg:: /auto_examples/06_matching/images/sphx_glr_plot_06_downstream_tables_001.png
   :alt: raw ids: column spread, matched ids: column spread
   :srcset: /auto_examples/06_matching/images/sphx_glr_plot_06_downstream_tables_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

                raw_ids_sd  matched_ids_sd
    H1_enh_lap       1.059           0.304
    H2_enh_lap       0.653           0.621
    H3_enh_lap         NaN           1.221


.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>raw_ids_sd</th>
          <th>matched_ids_sd</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>H1_enh_lap</th>
          <td>1.059</td>
          <td>0.304</td>
        </tr>
        <tr>
          <th>H2_enh_lap</th>
          <td>0.653</td>
          <td>0.621</td>
        </tr>
        <tr>
          <th>H3_enh_lap</th>
          <td>NaN</td>
          <td>1.221</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 158-165

Leaving outliers unnamed with ``max_distance``
----------------------------------------------
By default every habitat is named, however far it is from its
prototype. With ``max_distance`` (Euclidean, in feature units) a habitat
farther than that from every free prototype keeps a subject-local id
above ``K`` and ``prototype_id`` NA. It then belongs to no shared
column. Off by default: when set, report how many habitats it removed.

.. GENERATED FROM PYTHON SOURCE LINES 165-177

.. code-block:: Python

    strict = align_habitat_maps_to_prototypes(maps, models=models, max_distance=1.0)
    unnamed = strict.assignments[strict.assignments.prototype_id.isna()]
    print(unnamed.to_string(index=False))
    K = strict.prototypes.shape[0]
    for habitat_map in strict.habitat_maps:
        local = [k for k in habitat_map.habitat_ids if k > K]
        if local:
            print(f"{habitat_map.subject_id}: ids {local} are subject-local (not a cohort column)")

    fig = plot_prototype_matching([np.asarray(m.centroids) for m in models], strict)
    fig.savefig("out/downstream_max_distance.png", dpi=150, bbox_inches="tight")
    plt.show()



.. image-sg:: /auto_examples/06_matching/images/sphx_glr_plot_06_downstream_tables_002.png
   :alt: Habitats named by shared prototypes
   :srcset: /auto_examples/06_matching/images/sphx_glr_plot_06_downstream_tables_002.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    subject_id  habitat_id  prototype_id  distance
       subj003           2          <NA>       NaN
       subj005           1          <NA>       NaN
    subj003: ids [4] are subject-local (not a cohort column)
    subj005: ids [4] are subject-local (not a cohort column)





.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 33.703 seconds)


.. _sphx_glr_download_auto_examples_06_matching_plot_06_downstream_tables.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_06_downstream_tables.ipynb <plot_06_downstream_tables.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_06_downstream_tables.py <plot_06_downstream_tables.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_06_downstream_tables.zip <plot_06_downstream_tables.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_