Metadata-Version: 2.4
Name: MLBpy
Version: 0.2.0
Summary: Multi-layer Bundling of features or samples using spectral clustering consensus
Author: Mehran Fazli
License-Expression: MIT
Project-URL: Homepage, https://github.com/mehranfazli/MLBpy
Project-URL: Repository, https://github.com/mehranfazli/MLBpy
Project-URL: Paper, https://doi.org/10.1007/s11538-024-01335-8
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.24
Requires-Dist: scipy>=1.10
Requires-Dist: scikit-learn>=1.3
Provides-Extra: dataframe
Requires-Dist: pandas>=2.0; extra == "dataframe"
Provides-Extra: plot
Requires-Dist: matplotlib>=3.7; extra == "plot"
Requires-Dist: networkx>=3.4; extra == "plot"
Requires-Dist: plotly>=6.1; extra == "plot"
Requires-Dist: kaleido>=1.0; extra == "plot"
Provides-Extra: simulation
Requires-Dist: matplotlib>=3.7; extra == "simulation"
Requires-Dist: networkx>=3.0; extra == "simulation"
Provides-Extra: reconstruction
Requires-Dist: networkx>=3.0; extra == "reconstruction"
Provides-Extra: test
Requires-Dist: pytest>=7.4; extra == "test"
Requires-Dist: pandas>=2.0; extra == "test"
Requires-Dist: matplotlib>=3.7; extra == "test"
Requires-Dist: networkx>=3.4; extra == "test"
Requires-Dist: plotly>=6.1; extra == "test"
Provides-Extra: all
Requires-Dist: pandas>=2.0; extra == "all"
Requires-Dist: matplotlib>=3.7; extra == "all"
Requires-Dist: networkx>=3.4; extra == "all"
Requires-Dist: plotly>=6.1; extra == "all"
Requires-Dist: kaleido>=1.0; extra == "all"
Dynamic: license-file

# MLBpy

**MLBpy** is a Python implementation of Multi-layer Bundling (MLB), a method
for finding correlation structure at multiple scales in high-dimensional data.
Instead of requiring the user to choose one spectral-clustering solution, MLB
combines several prominent solutions and identifies groups of items whose
relationships persist across them. The analyzed items may be features, such as
genes or proteins, or samples, such as patients or experimental subjects.

## How Multi-layer Bundling works

MLB begins with an item-by-item affinity matrix. In the original method this is
the absolute Pearson correlation matrix; MLBpy also supports other similarity
measures and signed associations. The normalized graph-Laplacian spectrum is
then calculated. A large gap between consecutive eigenvalues indicates a
prominent choice for the number of spectral clusters. MLB ranks these
eigengaps and computes a separate spectral-clustering partition, called a
**clustering regime**, for each selected cluster count.

The clustering regimes are alternative views of the same data. They are
ordered by eigengap prominence, not by their number of clusters, so they do not
have to be nested and their cluster counts do not have to increase from one
regime to the next. MLB creates a hierarchy by intersecting these regimes.

For an item `x`, let `c_i(x)` be its cluster label in clustering regime `i`.
At bundle layer `l`, two items `x` and `y` belong to the same bundle exactly
when

```text
c_i(x) = c_i(y) for every regime i from 1 through l.
```

Equivalently, each layer groups items having the same accumulated membership
signature `(c_1, c_2, ..., c_l)`. Layer 1 therefore uses the most prominent
clustering regime. Layer 2 retains items together only if they co-cluster in
both of the first two regimes, and each later layer adds one more co-clustering
requirement.

This cumulative intersection gives the layers their global-to-local structure.
A bundle may remain unchanged or split when a regime is added, but separate
bundles can never merge. Consequently, the number of bundles cannot decrease
with layer number, while individual bundles cannot grow. Items that remain
together through several layers form a persistent bundle: their relationship
is supported across several prominent spectral projections rather than only
one selected cluster count.

There is no universally correct final layer. Very deep layers can eventually
fragment the data into many small or singleton bundles. The user first chooses
a preliminary maximum layer, then examines the Laplacian spectrum, bundle-size
curves, and Sankey diagram. A useful layer is often near a stable region in the
bundle-size distribution, before further layers add little structure or cause
excessive fragmentation.

After selecting a layer, the bundle labels can be projected onto the original
item network. MLBpy can also calculate the mean affinity between every pair of
bundles and use those scores to reconstruct a bundle-level network. Bundle
formation and network reconstruction are separate steps: the former defines
which items persist together, while the latter summarizes relationships among
the resulting bundles.

MLBpy is based on:

> Fazli, M., Bertram, R. & Striegel, D.A. Multi-layer Bundling as a New
> Approach for Determining Multi-scale Correlations Within a High-Dimensional
> Dataset. *Bulletin of Mathematical Biology* **86**, 105 (2024).
> [https://doi.org/10.1007/s11538-024-01335-8](https://doi.org/10.1007/s11538-024-01335-8)

## Features

- Feature-level and sample-level bundling
- Pearson, Spearman, Euclidean, and cosine similarity
- Signed and unsigned associations
- Normalized graph-Laplacian analysis
- Eigengap-based cluster selection
- Spectral clustering and multi-layer consensus bundling
- Persistent bundle-size visualization
- Interactive bundle-persistence Sankey visualization
- ForceAtlas2 network visualization
- Strongest-neighbor bundle-network reconstruction
- Threshold-based bundle-network reconstruction
- Deterministic synthetic signed-network simulation
- Raw-data and precomputed-similarity input

## Installation

Install the package from the repository directory:

```bash
python -m pip install -e .
```

Install plotting, simulation, and reconstruction dependencies:

```bash
python -m pip install -e ".[all]"
```

For development and testing:

```bash
python -m pip install -e ".[all,test]"
python -m pytest
```

## Input format

Raw input data must always have samples in rows and features in columns:

| Sample | Feature A | Feature B | Feature C |
|---|---:|---:|---:|
| Sample 1 | 0.24 | 1.08 | 0.51 |
| Sample 2 | 0.31 | 0.92 | 0.66 |
| Sample 3 | 0.18 | 1.15 | 0.47 |

MLBpy can bundle either the columns (features) or rows (samples).

## Tutorial: bundle features

Feature analysis is the default. DataFrame column names are preserved in their
exact input order:

```python
from MLBpy import MLB

model = MLB(
    method="pearson",
    signed=False,
    n_layers=7,
    random_state=18,
)

model.fit(dataframe, axis="features")
```

For a NumPy array, provide feature names in the same order as its columns:

```python
model.fit(X, axis="features", feature_names=feature_names)
```

## Tutorial: bundle samples

To study relationships among samples, select the sample axis:

```python
model.fit(dataframe, axis="samples")
```

MLBpy transposes the data internally. DataFrame index values become the sample
names and remain in their original order.

For an array, provide sample names or subject codes in row order:

```python
model.fit(X, axis="samples", sample_names=sample_codes)
```

MLBpy rejects duplicate names and name lists whose length does not match the
corresponding data axis. The analyzed target and ordered names are stored in:

```python
model.analysis_axis_
model.item_names_in_
```

Depending on the target, the model also stores `feature_names_in_` or
`sample_names_in_`.

## Choose the preliminary maximum layer

The user first chooses the maximum number of MLB layers to examine:

```python
maximum_layer = 7

model = MLB(
    method="pearson",
    signed=False,
    n_layers=maximum_layer,
    random_state=18,
)
```

After fitting, the selected spectral cluster counts are available from
`model.cluster_counts_`.

## Review the preliminary figures

Inspect the Laplacian spectrum and bundle-size curves before selecting the
final bundle layer:

```python
import matplotlib.pyplot as plt
from MLBpy import plot_bundle_counts, plot_spectrum

figure, axes = plt.subplots(1, 2, figsize=(16, 6))
plot_spectrum(model, axes[0])
plot_bundle_counts(model, axes[1])
figure.tight_layout()
plt.show()
```

### Laplacian spectrum

The spectrum displays the ordered normalized graph-Laplacian eigenvalues and
highlights the eigengap-selected cluster counts. Its default maximum displayed
index is:

```python
min(item_number - 1, max(30, ceil(item_number / 10)))
```

The complete spectrum remains available in `model.eigenvalues_`.

### Persistent bundles by size

For each layer, the bundle plot counts bundles whose sizes are strictly greater
than:

```text
0, 1, 2, 5, 10, 20, 30, and 50
```

These figures provide the preliminary information for choosing a layer.

## Select a bundle layer

After reviewing the figures, select one layer between 1 and the preliminary
maximum:

```python
selected_layer = 5
bundles = model.get_bundles(selected_layer)
labels = model.get_labels(selected_layer)
```

Print the selected bundle list:

```python
for bundle_id, names in bundles.items():
    print(f"Bundle {bundle_id}, size {len(names)}")
    print(names)
```

For feature analysis, the lists contain feature names. For sample analysis,
they contain sample names or subject codes.

## Visualize persistence through the selected layer

Display how bundles split or persist across every consecutive layer from layer
1 through the user's selected layer:

```python
from MLBpy import plot_bundle_sankey

sankey = plot_bundle_sankey(model, selected_layer)
sankey.show()
```

Each node is a bundle at one layer. Each link connects bundles in consecutive
layers, and its width equals the number of features or samples shared by the
source and target bundles. Hovering over a node shows its bundle size;
hovering over a link shows its shared-item count.

Save the interactive figure as HTML:

```python
sankey.write_html("bundle_persistence_sankey.html")
```

The `plot` installation option includes Kaleido, so the same figure can be
exported as a high-resolution static image. Kaleido 1 also requires Chrome or
Chromium to be installed on the computer:

```python
sankey.write_image(
    "bundle_persistence_sankey.png",
    width=1000,
    height=500,
    scale=10,
)
```

Set `show_labels=False` in `plot_bundle_sankey()` to hide labels while retaining
the bundle and layer information in the interactive hover text.

## Draw the selected-layer network

Draw an item network and color it using the selected bundles. ForceAtlas2 is
the default layout:

```python
import matplotlib.pyplot as plt
from MLBpy import plot_bundle_network

figure, axis = plt.subplots(figsize=(14, 11))
plot_bundle_network(
    model,
    adjacency,
    layer=selected_layer,
    ax=axis,
    layout="forceatlas2",
    layout_seed=18,
)
figure.tight_layout()
plt.show()
```

The adjacency-matrix order must exactly match `model.item_names_in_`: column
order for feature analysis or row order for sample analysis.

The `layout` option accepts `"forceatlas2"` (default), `"spring"`,
`"kamada_kawai"`, `"circular"`, or `"spectral"`.

## Use a precomputed similarity matrix

Analyze a feature-by-feature similarity matrix:

```python
model.fit_similarity(
    feature_similarity,
    axis="features",
    feature_names=feature_names,
)
```

Analyze a sample-by-sample similarity matrix:

```python
model.fit_similarity(
    sample_similarity,
    axis="samples",
    sample_names=sample_codes,
)
```

The matrix must be square, symmetric, finite, and ordered consistently with
the supplied names.

## Bundle-network reconstruction

After choosing a layer, aggregate item similarities into bundle-pair scores:

```python
from MLBpy import bundle_similarity_scores, reconstruct_bundle_network

labels = model.get_labels(selected_layer)
bundle_scores, bundle_ids = bundle_similarity_scores(
    model.similarity_matrix_,
    labels,
)
```

### Option 1: strongest-neighbor reconstruction

This corresponds to `reversed_synthetic_bundle_net()` in the original code:

```python
strongest_network = reconstruct_bundle_network(
    bundle_scores,
    bundle_ids=bundle_ids,
    method="strongest",
)
```

It selects each bundle's strongest remaining positive edge, preserves tied
strongest scores, and connects separate components using the strongest
remaining cross-component scores. It generally produces a sparse, connected
bundle network.

The explicit function `reconstruct_strongest_network()` is also available.

### Option 2: threshold reconstruction

This corresponds to `reversed_synthetic_bundle_net_thr()`:

```python
threshold_network = reconstruct_bundle_network(
    bundle_scores,
    bundle_ids=bundle_ids,
    method="threshold",
    threshold=0.5,
)
```

It adds every pair whose score is strictly greater than the threshold truncated
to four decimal places. The resulting network is not forced to be connected.
The explicit function `reconstruct_threshold_network()` is also available.

Persistent-homology threshold selection is not yet part of the package; the
threshold is currently supplied by the user.

### Draw the reconstructed bundle networks

Both reconstructed graphs use ForceAtlas2 by default. Node colors correspond
to bundle IDs, node areas represent bundle sizes, and edge widths represent
bundle-pair similarity. The same alternative `layout` values are available:

```python
import matplotlib.pyplot as plt
from MLBpy import plot_reconstructed_bundle_network

bundle_sizes = {
    bundle_id: len(names)
    for bundle_id, names in model.get_bundles(selected_layer).items()
}

figure, axes = plt.subplots(1, 2, figsize=(18, 8))
plot_reconstructed_bundle_network(
    strongest_network,
    axes[0],
    bundle_sizes=bundle_sizes,
    layer=selected_layer,
    layout="forceatlas2",
    layout_seed=18,
)
plot_reconstructed_bundle_network(
    threshold_network,
    axes[1],
    bundle_sizes=bundle_sizes,
    layer=selected_layer,
    layout="forceatlas2",
    layout_seed=18,
)
figure.tight_layout()
plt.show()
```

Overlay both reconstructions in one network. Threshold edges are drawn first
in gray with `alpha=0.4`; strongest-neighbor edges are drawn above them in red.
Node positions come from the strongest-neighbor network using the selected
layout. Edge thickness increases with bundle-pair similarity in both layers:

```python
from MLBpy import plot_combined_bundle_network

figure, axis = plt.subplots(figsize=(11, 9))
plot_combined_bundle_network(
    strongest_network,
    threshold_network,
    axis,
    bundle_sizes=bundle_sizes,
    layer=selected_layer,
    layout="forceatlas2",
    layout_seed=18,
)
figure.tight_layout()
plt.show()
```

## Synthetic-network example

Generate the extended-network simulation used during development:

```python
from MLBpy import simulate_network

simulation = simulate_network(
    n_samples=400,
    n_features=500,
    beta=0.6,
    noise_sd=0.2,
    seed=18,
    extra_edges=10,
)
```

The result provides:

```python
simulation.data
simulation.affinity
simulation.adjacency
simulation.signed_adjacency
simulation.feature_names
simulation.sample_names
```

Run feature-level MLB on the simulated affinity:

```python
model.fit_similarity(
    simulation.affinity,
    axis="features",
    feature_names=simulation.feature_names,
)
```

Run sample-level MLB on the simulated data:

```python
sample_model = MLB(
    method="pearson",
    signed=False,
    n_layers=7,
    random_state=18,
)
sample_model.fit(
    simulation.data,
    axis="samples",
    sample_names=simulation.sample_names,
)
```

## Complete Jupyter workflow

Open `examples/MLBpy_complete_workflow.ipynb` in Jupyter Notebook or JupyterLab
and run its cells in order. It contains two complete examples.

The first example asks for simulation parameters and the preliminary maximum
layer, generates a 400 × 500 dataset, displays the preliminary figures, asks
for a specific layer, displays a Sankey plot, draws the ForceAtlas2-colored
item network, reconstructs both bundle-network options, and saves all results.

The second example loads scikit-learn's Wisconsin Diagnostic Breast Cancer
dataset and demonstrates sample-level MLB. The input is a DataFrame with 569
samples in rows and 30 measurements in columns. Ordered sample codes are stored
in the DataFrame index, the measurements are standardized, and the model is fit
with `axis="samples"` using Euclidean similarity. Diagnosis is excluded from
the fit and used only afterward to summarize malignant and benign membership
within the discovered bundles. Two separate PCA plots display the same samples
colored by the selected MLB bundles and by malignant/benign diagnosis for
post-hoc visual comparison.

Both examples retain the bundle assignments from layer 1 through the
preliminary maximum, the selected bundle lists, diagnostic figures, and
interactive Sankey diagrams.

## Fitted results

Important fitted attributes include:

```python
model.analysis_axis_
model.item_names_in_
model.similarity_matrix_
model.eigenvalues_
model.cluster_counts_
model.partition_labels_
model.consensus_matrices_
model.layer_labels_
model.bundle_matrices_
model.bundles_by_layer_
model.labels_
model.bundles_
```

## Public API

| Name | Purpose |
|---|---|
| `MLB` | Fit the multi-layer bundling model to features, samples, or a precomputed similarity matrix. |
| `compute_similarity()` | Calculate Pearson, Spearman, Euclidean, or cosine item similarity. |
| `normalized_laplacian_spectrum()` | Calculate the ordered normalized graph-Laplacian spectrum. |
| `eigengap_candidates()` | Rank cluster counts using the largest eigengaps. |
| `simulate_network()` | Generate deterministic signed-network test data with samples in rows and features in columns. |
| `bundle_similarity_scores()` | Aggregate item similarities into mean between-bundle scores. |
| `reconstruct_strongest_network()` | Build the sparse strongest-neighbor bundle network. |
| `reconstruct_threshold_network()` | Connect all bundle pairs strictly above a user threshold. |
| `reconstruct_bundle_network()` | Select either reconstruction strategy through one interface. |
| `plot_spectrum()` | Plot the preliminary Laplacian spectrum with a limited eigenvalue index. |
| `plot_bundle_counts()` | Plot persistent bundle counts for the MLB_core size limits. |
| `plot_bundle_sankey()` | Show bundle persistence from layer 1 through the selected layer. |
| `plot_bundle_network()` | Draw the item network colored by selected-layer bundles. |
| `plot_reconstructed_bundle_network()` | Draw one weighted reconstructed bundle network. |
| `plot_combined_bundle_network()` | Overlay threshold edges behind strongest-neighbor edges using strongest-network positions. |

## Reproducibility

Set `random_state` for spectral clustering and `seed` for simulation. Identical
inputs, parameters, and seeds produce reproducible results.

## Citation

If you use MLB or MLBpy in your research, please cite:

Fazli, M., Bertram, R. & Striegel, D.A. Multi-layer Bundling as a New Approach
for Determining Multi-scale Correlations Within a High-Dimensional Dataset.
*Bulletin of Mathematical Biology* **86**, 105 (2024).
[https://doi.org/10.1007/s11538-024-01335-8](https://doi.org/10.1007/s11538-024-01335-8)

### BibTeX

```bibtex
@article{Fazli2024MultiLayerBundling,
  author  = {Fazli, M. and Bertram, R. and Striegel, D. A.},
  title   = {Multi-layer Bundling as a New Approach for Determining Multi-scale Correlations Within a High-Dimensional Dataset},
  journal = {Bulletin of Mathematical Biology},
  volume  = {86},
  pages   = {105},
  year    = {2024},
  doi     = {10.1007/s11538-024-01335-8},
  url     = {https://doi.org/10.1007/s11538-024-01335-8}
}
```

## Project status

MLBpy provides a tested implementation of the core MLB workflow and associated
simulation, visualization, and bundle-network reconstruction tools. Users
should report the input orientation, similarity method, sign treatment,
preliminary maximum layer, selected layer, and reconstruction method.

## License

MLBpy is distributed under the MIT License.
