Metadata-Version: 2.4
Name: otica
Version: 0.1.3
Summary: Optimal Transport Independent Component Analysis
Author: Félix Laplante
License-Expression: GPL-3.0-only
Project-URL: Source, https://github.com/felixlaplante0/otica
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering :: Artificial Intelligence
Classifier: Topic :: Scientific/Engineering :: Mathematics
Classifier: Operating System :: POSIX :: Linux
Classifier: Operating System :: MacOS
Classifier: Operating System :: Microsoft :: Windows
Requires-Python: >=3.11
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy
Requires-Dist: scipy
Requires-Dist: scikit-learn
Provides-Extra: docs
Requires-Dist: furo; extra == "docs"
Requires-Dist: sphinx; extra == "docs"
Provides-Extra: test
Requires-Dist: pytest>=8; extra == "test"
Requires-Dist: pytest-cov>=5; extra == "test"
Dynamic: license-file

# 📊 Contrast-Free ICA

[![codecov](https://codecov.io/gh/felixlaplante0/otica/graph/badge.svg)](https://codecov.io/gh/felixlaplante0/otica)

**otica** is a Python package for linear independent component analysis (ICA) based on optimal transport. It recovers latent sources by maximizing their empirical squared 2-Wasserstein distances to the standard Gaussian, using a fixed non-Gaussianity criterion that requires no user-chosen contrast function or nonlinearity.

---

## ✨ Features

- **Contrast-free source separation**: Uses the squared 2-Wasserstein distance to the standard Gaussian as a fixed non-Gaussianity criterion.
- **Exact empirical objective**: Computes the one-dimensional Wasserstein criterion directly from ordered samples and Gaussian quantiles, without density estimation.
- **Riemannian optimization**: Optimizes the whitened ICA objective on the orthogonal group with a Picard-style limited-memory BFGS method and Armijo backtracking.
- **Dimension reduction**: Supports extraction of a specified number of components through principal-component whitening.
- **Flexible initialization**: Accepts FastICA, random, or user-provided initial unmixing matrices through `w_init`.
- **scikit-learn integration**: Native `BaseEstimator` integration with the standard transformer API, including `fit`, `transform`, `fit_transform`, and `inverse_transform`.

---

## ⚡ Method

For observations generated by the linear ICA model $X = S A^\top$, **otica** first centers and whitens the data. Let $Z \in \mathbb{R}^{d}$ be the whitened random vector. For an orthogonal unmixing matrix $W \in \mathbb{R}^{d \times d}$, the population objective is

$$
F(W) = \sum_{k = 1}^{d} \mathcal{W}_2\left( (Z W^\top)_k, \mathcal{N}(0, 1) \right)^2, \quad W W^\top = I_d.
$$

Given $n$ whitened observations collected as the rows of $Z \in \mathbb{R}^{n \times d}$, let $Y = Z W^\top$. OTICA maximizes

$$
\widehat{F}_n(W) = \sum_{k = 1}^{d} \mathcal{W}_2\left( \frac{1}{n} \sum_{i = 1}^{n} \delta_{Y_{ik}}, \mathcal{N}(0, 1) \right)^2.
$$

For each component $k$, the implementation sorts the entries of the $k$-th column as $Y_{(1)k} \leq \cdots \leq Y_{(n)k}$ and matches them with the corresponding standard-Gaussian rank statistics. Under the usual ICA assumptions, including mutually independent sources with at most one Gaussian component, the population objective identifies the sources up to permutation, sign, and scale.

---

## 🚀 Installation

```bash
pip install otica
```

## 🔧 Usage

### Example

The following example generates three independent non-Gaussian signals, mixes them linearly, and recovers them with `OTICA`. Because ICA is identifiable only up to permutation and sign, recovery is evaluated using the best absolute correlation for each true source.

```python
import matplotlib.pyplot as plt
import numpy as np
from otica import OTICA

rng = np.random.default_rng(42)
n_samples = 5_000
time = np.linspace(0.0, 8.0, n_samples)

# Generate independent, non-Gaussian latent sources.
sources = np.column_stack(
    [
        rng.laplace(size=n_samples),
        rng.uniform(-np.sqrt(3.0), np.sqrt(3.0), size=n_samples),
        rng.standard_t(df=5, size=n_samples) * np.sqrt(3.0 / 5.0),
    ]
)

# Mix the sources into the observed signals.
mixing = np.array(
    [
        [1.0, 0.5, -0.2],
        [0.2, 1.0, 0.4],
        [-0.4, 0.1, 1.0],
    ]
)
X = sources @ mixing.T

# Fit OTICA and recover the latent components.
model = OTICA(random_state=42)
estimated_sources = model.fit_transform(X)

correlations = np.corrcoef(sources.T, estimated_sources.T)[:3, 3:]
best_indices = np.abs(correlations).argmax(axis=1)
best_correlations = correlations[np.arange(3), best_indices]
recovered_sources = estimated_sources[:, best_indices] * np.sign(best_correlations)
print("Best absolute correlation per source:", np.abs(best_correlations))

fig, axes = plt.subplots(3, 2, sharex=True, figsize=(12, 6))
for component in range(3):
    axes[component, 0].plot(time[:500], sources[:500, component])
    axes[component, 1].plot(
        time[:500], recovered_sources[:500, component], color="tab:orange"
    )
    axes[component, 0].set_ylabel(f"Source {component + 1}")

axes[0, 0].set_title("True sources")
axes[0, 1].set_title("Recovered sources")
axes[-1, 0].set_xlabel("Time")
axes[-1, 1].set_xlabel("Time")
fig.tight_layout()
plt.show()
```

---

## 📖 Learn More

For the mathematical formulation, configuration details, and API reference, visit [otica's documentation](https://felixlaplante0.github.io/otica).
