diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md new file mode 100644 index 0000000..6a348f9 --- /dev/null +++ b/.github/CHANGELOG.md @@ -0,0 +1,5 @@ +# Version 0.1.0a1 + +#### New Features + +* Initial release of `kalelinear`: knowledge-aware linear and kernel methods for multi-source/multi-view learning, including transfer learning, domain adaptation, manifold regularization, and group-aware estimators and transformers. diff --git a/.github/workflows/changelog.yml b/.github/workflows/changelog.yml new file mode 100644 index 0000000..e006387 --- /dev/null +++ b/.github/workflows/changelog.yml @@ -0,0 +1,31 @@ +# This workflow will generate a log of changes automatically upon a new release. +# See https://github.com/marketplace/actions/changelog-ci + +name: changelog + +on: + pull_request: + types: [opened] + +jobs: + log-changes: + name: Log changes + runs-on: ubuntu-latest + steps: + - name: Checkout code + uses: actions/checkout@v4 + with: + repository: ${{ github.event.pull_request.head.repo.full_name }} + ref: ${{ github.event.pull_request.head.ref }} + # Keep checkout shallow; changelog-ci handles unshallow internally. + # Using fetch-depth: 0 causes changelog-ci to fail with: + # "fatal: --unshallow on a complete repository does not make sense" + fetch-depth: 1 + token: ${{ secrets.GITHUB_TOKEN }} + - name: Run changelog + uses: saadmk11/changelog-ci@v1.2.0 + with: + changelog_filename: .github/CHANGELOG.md + config_file: .github/changelog-ci-config.json + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} diff --git a/.github/workflows/codeql-analysis.yml b/.github/workflows/codeql-analysis.yml new file mode 100644 index 0000000..95e4c21 --- /dev/null +++ b/.github/workflows/codeql-analysis.yml @@ -0,0 +1,70 @@ +# For most projects, this workflow file will not need changing; you simply need +# to commit it to your repository. +# +# You may wish to alter this file to override the set of languages analyzed, +# or to provide custom queries or build logic. +# +# ******** NOTE ******** +# We have attempted to detect the languages in your repository. Please check +# the `language` matrix defined below to confirm you have the correct set of +# supported CodeQL languages. +# +name: codeql-analysis + +on: + push: + branches: [ main ] + pull_request: + # The branches below must be a subset of the branches above + branches: [ main ] + schedule: + - cron: '15 18 * * 5' + +jobs: + analyze: + name: Analyze + runs-on: ubuntu-latest + permissions: + actions: read + contents: read + security-events: write + + strategy: + fail-fast: false + matrix: + language: [ 'python' ] + # CodeQL supports [ 'cpp', 'csharp', 'go', 'java', 'javascript', 'python', 'ruby' ] + # Learn more about CodeQL language support at https://git.io/codeql-language-support + + steps: + - name: Checkout repository + uses: actions/checkout@v4 + + # Initializes the CodeQL tools for scanning. + - name: Initialize CodeQL + uses: github/codeql-action/init@v3 + with: + languages: ${{ matrix.language }} + # If you wish to specify custom queries, you can do so here or in a config file. + # By default, queries listed here will override any specified in a config file. + # Prefix the list here with "+" to use these queries and those in the config file. + # queries: ./path/to/local/query, your-org/your-repo/queries@main + + # Autobuild attempts to build any compiled languages (C/C++, C#, or Java). + # If this step fails, then you should remove it and run the build manually (see below) + - name: Autobuild + uses: github/codeql-action/autobuild@v3 + + # ℹī¸ Command-line programs to run using the OS shell. + # 📚 https://git.io/JvXDl + + # ✏ī¸ If the Autobuild fails above, remove it and uncomment the following three lines + # and modify them (or add more) to build your code if your project + # uses a compiled language + + #- run: | + # make bootstrap + # make release + + - name: Perform CodeQL Analysis + uses: github/codeql-action/analyze@v3 diff --git a/.github/workflows/pre-commit.yml b/.github/workflows/pre-commit.yml new file mode 100644 index 0000000..5ba5112 --- /dev/null +++ b/.github/workflows/pre-commit.yml @@ -0,0 +1,28 @@ +# This workflow will run lint and many other pre-commit hooks. +# https://pre-commit.com/ + +name: pre-commit-check + +on: + push: + branches: [main] + pull_request: + branches: [main] + +jobs: + pre-commit-check: + name: Pre-commit checks including linting + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: "3.11" + - name: Install dependencies + run: | + pip install pre-commit isort + pre-commit install + - name: Run pre-commit checks including linting + run: | + pre-commit run --all-files diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml new file mode 100644 index 0000000..ebcc404 --- /dev/null +++ b/.github/workflows/release.yml @@ -0,0 +1,44 @@ +# This workflow will release a package on PyPI automatically when it is tagged/released. +# https://packaging.python.org/guides/publishing-package-distribution-releases-using-github-actions-ci-cd-workflows/ + +name: release + +on: + release: + types: [created] + +jobs: + build-n-publish: + name: Release on PyPI + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - name: Set up Python 3.11 + uses: actions/setup-python@v5 + with: + python-version: 3.11 + - name: Install pypa/build + run: >- + python -m + pip install + build + --user + - name: Build a binary wheel and a source tarball + run: >- + python -m + build + --sdist + --wheel + --outdir dist/ + . + - name: Publish distribution to Test PyPI + if: github.event.release.prerelease + uses: pypa/gh-action-pypi-publish@release/v1 + with: + password: ${{ secrets.TEST_PYPI_API_TOKEN }} + repository-url: https://test.pypi.org/legacy/ + - name: Publish distribution to PyPI + if: ${{ !github.event.release.prerelease }} + uses: pypa/gh-action-pypi-publish@release/v1 + with: + password: ${{ secrets.PYPI_API_TOKEN }} diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml new file mode 100644 index 0000000..1967143 --- /dev/null +++ b/.github/workflows/test.yml @@ -0,0 +1,95 @@ +# This workflow will install Python dependencies, run tests, and report the coverage with a variety of Python versions and OSs. +# For more information see: https://help.github.com/actions/language-and-framework-guides/using-python-with-github-actions + +name: test + +on: + push: + branches: [main] + pull_request: + branches: [main] + schedule: + # * is a special character in YAML, so you have to quote this string + - cron: "0 0 * * *" # every midnight + +jobs: + test: + name: Test (${{ matrix.os }}, python version ${{ matrix.python-version }}) + runs-on: ${{ matrix.os }} + strategy: + matrix: +# os: [ubuntu-latest, windows-latest] + os: [ubuntu-latest] + python-version: ["3.10", "3.11", "3.12"] # list of Python versions to test +# exclude: +# - os: windows-latest +# python-version: "3.10" + include: + - os: ubuntu-latest + path: ~/.cache/pip +# - os: windows-latest +# path: ~\AppData\Local\pip\Cache + + steps: + - uses: actions/checkout@v4 + - name: Set up Python using Miniconda + uses: conda-incubator/setup-miniconda@v3 + with: + auto-update-conda: true + python-version: ${{ matrix.python-version }} + miniconda-version: latest + + - name: Cache pip dependencies + id: cache_pip + uses: actions/cache@v4 + with: + path: ${{ matrix.path }} + key: ${{ runner.os }}-python${{ matrix.python-version }}-pip-20250302-${{ hashFiles('**/setup.py') }} + restore-keys: | + ${{ runner.os }}-python${{ matrix.python-version }}-pip-20250302- + # We have used a softer matching strategy for the full hash of setup.py, as recommended by GitHub. + # See: https://github.com/davronaliyev/Cache-dependencies-in-GitHub-Actions/blob/main/examples.md#python---pip + # This restores the cache first and then downloads any changed packages to avoid updating the cache with + # every change to the setup.py file, thus reducing the storage requirements of GitHub Action. + # We set a date tag to the cache key to show the updated date of the cache. We can update this date tag to + # generate new cache after every major changes in setup.py. + + - name: Install project and dev dependencies + run: | + pip install --no-build-isolation -e .[dev] + shell: bash -l {0} + + - name: Cache downloaded test data + id: cache_data + uses: actions/cache@v4 + with: + path: tests/test_data + key: ${{ runner.os }}-python${{ matrix.python-version }}-data-${{ hashFiles('tests/download_test_data.py') }} + restore-keys: | + ${{ runner.os }}-python${{ matrix.python-version }}-data-${{ hashFiles('tests/download_test_data.py') }} + # Use strict matching for the hash of download_test_data.py, as we want to update the cache whenever the file changes. + + - name: Download test data + if: steps.cache_data.outputs.cache-hit != 'true' + run: | + python tests/download_test_data.py + shell: bash -l {0} + + - name: Run tests with thread limits + id: run_tests + run: | + export OMP_NUM_THREADS=1 + export MKL_NUM_THREADS=1 + export NUMEXPR_NUM_THREADS=1 + pytest --nbmake --nbmake-timeout=3000 --cov=kalelinear + shell: bash -l {0} + + - name: Determine coverage + run: | + coverage xml + shell: bash -l {0} + + - name: Report coverage + uses: codecov/codecov-action@v4 + with: + token: ${{ secrets.CODECOV_TOKEN }} diff --git a/README.md b/README.md index f32c024..d972956 100644 --- a/README.md +++ b/README.md @@ -2,19 +2,18 @@ kalelinear logo

- - +[![tests](https://github.com/pykale/linear/workflows/test/badge.svg)](https://github.com/pykale/linear/actions/workflows/test.yml) +[![codecov](https://codecov.io/gh/pykale/linear/branch/main/graph/badge.svg)](https://codecov.io/gh/pykale/linear) [![GitHub license](https://img.shields.io/badge/license-MIT-blue.svg)](https://github.com/pykale/linear/blob/main/LICENSE) [![Python](https://img.shields.io/badge/python-3.10%20%7C%203.11%20%7C%203.12-blue)](https://www.python.org) [![PyPI version](https://img.shields.io/pypi/v/kalelinear?color=blue)](https://pypi.org/project/kalelinear/) [![PyPI downloads](https://pepy.tech/badge/kalelinear)](https://pepy.tech/project/kalelinear) - `kalelinear` is a Python library for learning harmonized or individualized models from multi-source/multi-view data in linear or reproducing kernel Hilbert spaces (RKHS). It provides NumPy-based methods for leveraging related data distributions and structural assumptions, including transfer learning, domain adaptation, manifold regularization, and group-aware learning, through a [`scikit-learn`](https://github.com/scikit-learn/scikit-learn) style API. The package is part of the [PyKale](https://github.com/pykale/pykale) ecosystem and focuses on linear and kernel methods for data characterized by covariates (e.g., domain labels, group labels, side information), unlabeled target samples, or tensor structures. -## Methods and features +## What's included - Transformer models for learning feature embeddings: - Multilinear Principal Component Analysis (MPCA) [[1](#references)] @@ -28,24 +27,9 @@ The package is part of the [PyKale](https://github.com/pykale/pykale) ecosystem - Covariate Independence Regularized Learning Framework (CoIRSVM, CoIRLS) [[8](#references)][[9](#references)] - Group-specific Discriminant Analysis (GSDA) [[9](#references)][[10](#references)] - NumPy-compatible inputs and outputs. -- scikit-learn style `fit`, `transform`, `predict`, `fit_transform`, and - `fit_predict` workflows where applicable. +- scikit-learn style `fit`, `transform`, `predict`, `fit_transform`, and `fit_predict` workflows where applicable. - Optional covariate encoding for categorical domain or group labels. -## Installation - -Install the released package from PyPI: - -```bash -pip install kalelinear -``` - -Install from a local checkout for development: - -```bash -pip install -e ".[dev]" -``` - `kalelinear` requires Python 3.10 or later. Core dependencies include: - [NumPy](http://www.numpy.org/) @@ -56,112 +40,23 @@ pip install -e ".[dev]" - [cvxopt](http://cvxopt.org/) - [osqp](https://osqp.org/) -## Quick Start - -### Learn a Domain-Invariant Embedding - -```python -import numpy as np -from kalelinear.transformer import TCA - -X = np.array( - [ - [-2.0, -1.8], - [-1.8, -2.1], - [1.9, 1.7], - [2.1, 2.0], - [-1.4, -1.2], - [-1.2, -1.1], - [1.2, 1.1], - [1.4, 1.3], - ] -) -domain_labels = np.array([0, 0, 0, 0, 1, 1, 1, 1]) - -transformer = TCA(n_components=2) -z = transformer.fit_transform(X, covariates=domain_labels, target_covariate=1) - -z_source = z[domain_labels == 0] -z_target = z[domain_labels == 1] -``` - -TCA, JDA, and BDA take domain labels through `covariates`. They do not accept -separate source and target arrays; stack samples into one array and use -`target_covariate` to identify the target domain. - -### Use MIDA with Categorical Covariates - -```python -import numpy as np -from kalelinear.transformer import MIDA - -x = np.random.default_rng(0).normal(size=(8, 4)) -y = np.array([0, 0, 1, 1, 0, 0, 1, 1]) -domains = np.array(["source", "source", "source", "source", "target", "target", "target", "target"]) - -transformer = MIDA(n_components=2, covariate_encoder="onehot") -z = transformer.fit_transform(x, y=y, covariates=domains) -``` - -### Train a Domain Adaptation Classifier +## Getting started -For ARSVM and ARRLS, pass all source and target samples in `x`, labels for the -source samples in `y`, and a covariate vector identifying the target domain. +### Installation -```python -import numpy as np -from kalelinear.estimator import ARSVM - -x = np.array( - [ - [-2.2, -1.9], - [-1.9, -2.1], - [1.8, 2.1], - [2.0, 1.9], - [-1.4, -1.2], - [-1.1, -1.3], - [1.3, 1.1], - [1.5, 1.2], - ] -) - -source_labels = np.array([0, 0, 1, 1]) -domains = np.array([0, 0, 0, 0, 1, 1, 1, 1]) -x_target = x[domains == 1] - -clf = ARSVM() -clf.fit(x, source_labels, covariates=domains, target_covariate=1) -y_pred = clf.predict(x_target) -``` - -### Train a Manifold-Regularized Classifier - -LapSVM and LapRLS can use labeled source samples together with unlabeled target -samples. The labels array may contain only the labeled source examples. - -```python -import numpy as np -from kalelinear.estimator import LapSVM - -x_source = np.array([[-2.0, -1.8], [-1.8, -2.1], [1.9, 1.7], [2.1, 2.0]]) -ys = np.array([0, 0, 1, 1]) -x_target = np.array([[-1.4, -1.2], [-1.2, -1.1], [1.2, 1.1], [1.4, 1.3]]) - -x_train = np.vstack((x_source, x_target)) +Install the released package from PyPI: -clf = LapSVM(kernel="linear") -clf.fit(x_train, ys) -y_pred = clf.predict(x_target) +```bash +pip install kalelinear ``` -## Public API +Install from a local checkout for development: -```python -from kalelinear.transformer import BDA, JDA, MIDA, MPCA, TCA -from kalelinear.estimator import ARRLS, ARSVM, CoIRLS, CoIRSVM, GSDA, LapRLS, LapSVM +```bash +pip install -e ".[dev]" ``` -## Development +### Development From the root of the repository, run the following commands in your terminal: @@ -190,6 +85,21 @@ From the root of the repository, run the following commands in your terminal: sphinx-build -b html docs/source docs/build/html ``` +### Public API + +```python +from kalelinear.transformer import BDA, JDA, MIDA, MPCA, TCA +from kalelinear.estimator import ARRLS, ARSVM, CoIRLS, CoIRSVM, GSDA, LapRLS, LapSVM +``` + +Worked examples for the main transformers and estimators are collected in +[Tutorials](TUTORIALS.md): + +- Learn a domain-invariant embedding with TCA +- Use MIDA with categorical covariates +- Train a domain adaptation classifier (ARSVM, ARRLS) +- Train a manifold-regularized classifier (LapSVM, LapRLS) + # References [1] Lu, H., Plataniotis, K.N. and Venetsanopoulos, A.N., 2008. [MPCA: Multilinear principal component analysis of tensor objects](https://ieeexplore.ieee.org/abstract/document/4359192/). _IEEE Transactions on Neural Networks_, 19(1), pp.18-39. diff --git a/TUTORIALS.md b/TUTORIALS.md new file mode 100644 index 0000000..320e125 --- /dev/null +++ b/TUTORIALS.md @@ -0,0 +1,103 @@ +# Tutorials + +Worked examples for the `kalelinear` transformers and estimators. For +installation instructions and an API overview, see the +[README](README.md). + +## Learn a Domain-Invariant Embedding + +### Use TCA for Two-Domain Adaptation + +```python +import numpy as np +from kalelinear.transformer import TCA + +X = np.array( + [ + [-2.0, -1.8], + [-1.8, -2.1], + [1.9, 1.7], + [2.1, 2.0], + [-1.4, -1.2], + [-1.2, -1.1], + [1.2, 1.1], + [1.4, 1.3], + ] +) +domain_labels = np.array([0, 0, 0, 0, 1, 1, 1, 1]) + +transformer = TCA(n_components=2) +z = transformer.fit_transform(X, covariates=domain_labels, target_covariate=1) + +z_source = z[domain_labels == 0] +z_target = z[domain_labels == 1] +``` + +TCA, JDA, and BDA take domain labels through `covariates`. They do not accept +separate source and target arrays; stack samples into one array and use +`target_covariate` to identify the target domain. + +### Use MIDA with Categorical Domain Covariates + +```python +import numpy as np +from kalelinear.transformer import MIDA + +X = np.random.default_rng(0).normal(size=(8, 4)) +y = np.array([0, 0, 1, 1, 0, 0, 1, 1]) +domains = np.array(["source1", "source1", "source2", "source2", "target", "target", "target", "target"]) + +transformer = MIDA(n_components=2, covariate_encoder="onehot") +z = transformer.fit_transform(X, y=y, covariates=domains) +``` + +## Train a Domain Adaptation Classifier + +For ARSVM and ARRLS, pass all source and target samples in `X`, labels for the +source samples in `y`, and a covariate vector identifying the target domain. + +```python +import numpy as np +from kalelinear.estimator import ARSVM + +X = np.array( + [ + [-2.2, -1.9], + [-1.9, -2.1], + [1.8, 2.1], + [2.0, 1.9], + [-1.4, -1.2], + [-1.1, -1.3], + [1.3, 1.1], + [1.5, 1.2], + ] +) + +source_labels = np.array([0, 0, 1, 1]) +domains = np.array([0, 0, 0, 0, 1, 1, 1, 1]) +X_target = X[domains == 1] + +clf = ARSVM() +clf.fit(X, source_labels, covariates=domains, target_covariate=1) +y_pred = clf.predict(X_target) +``` + +## Train a Manifold-Regularized Classifier + +LapSVM and LapRLS can use labeled source samples together with unlabeled target +samples. The labels array may contain only the labeled source examples. + +```python +import numpy as np +from kalelinear.estimator import LapSVM + +X_source = np.array([[-2.0, -1.8], [-1.8, -2.1], [1.9, 1.7], [2.1, 2.0]]) +ys = np.array([0, 0, 1, 1]) +X_target = np.array([[-1.4, -1.2], [-1.2, -1.1], [1.2, 1.1], [1.4, 1.3]]) + +X_train = np.vstack((X_source, X_target)) + +clf = LapSVM(kernel="linear") +clf.fit(X_train, ys) +y_pred = clf.predict(X_target) +``` diff --git a/docs/source/conf.py b/docs/source/conf.py index 4abcaa7..0932ac4 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -1,4 +1,4 @@ -"""Sphinx configuration for Kale-Linear documentation.""" +"""Sphinx configuration for kalelinear documentation.""" from __future__ import annotations @@ -8,7 +8,7 @@ sys.path.insert(0, os.path.abspath("../..")) -project = "Kale-Linear" +project = "kalelinear" author = "The PyKale team" copyright = f"{datetime.now().year}, {author}" diff --git a/docs/source/index.rst b/docs/source/index.rst index 90ee5fd..343bc78 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -1,4 +1,4 @@ -Kale-Linear Documentation +kalelinear Documentation ========================= Getting Started @@ -11,7 +11,7 @@ Getting Started installation tutorial -Kale-Linear API +kalelinear API --------------- .. toctree:: @@ -23,7 +23,7 @@ Kale-Linear API api_estimators api_utilities -Kale-Linear APIs above are ordered following the machine learning pipeline, +kalelinear APIs above are ordered following the machine learning pipeline, i.e., feature embedding transformers, predictive estimators, and reusable utilities, rather than alphabetically. diff --git a/docs/source/tutorial.rst b/docs/source/tutorial.rst index 2ac9c36..a83442d 100644 --- a/docs/source/tutorial.rst +++ b/docs/source/tutorial.rst @@ -13,7 +13,7 @@ label is the target domain. import numpy as np from kalelinear.transformer import TCA - x = np.array( + X = np.array( [ [-2.0, -1.8], [-1.8, -2.1], @@ -28,7 +28,7 @@ label is the target domain. domain_labels = np.array([0, 0, 0, 0, 1, 1, 1, 1]) transformer = TCA(n_components=2) - z = transformer.fit_transform(x, covariates=domain_labels, target_covariate=1) + z = transformer.fit_transform(X, covariates=domain_labels, target_covariate=1) z_source = z[domain_labels == 0] z_target = z[domain_labels == 1] @@ -41,19 +41,19 @@ Use MIDA with Categorical Covariates import numpy as np from kalelinear.transformer import MIDA - x = np.random.default_rng(0).normal(size=(8, 4)) + X = np.random.default_rng(0).normal(size=(8, 4)) y = np.array([0, 0, 1, 1, 0, 0, 1, 1]) domains = np.array( ["source", "source", "source", "source", "target", "target", "target", "target"] ) transformer = MIDA(n_components=2, covariate_encoder="onehot") - z = transformer.fit_transform(x, y=y, covariates=domains) + z = transformer.fit_transform(X, y=y, covariates=domains) Train a Domain Adaptation Classifier ------------------------------------ -For ARSVM and ARRLS, pass all source and target samples in ``x``, labels for the +For ARSVM and ARRLS, pass all source and target samples in ``X``, labels for the source samples in ``y``, and a covariate vector identifying the target domain. .. code-block:: python @@ -61,7 +61,7 @@ source samples in ``y``, and a covariate vector identifying the target domain. import numpy as np from kalelinear.estimator import ARSVM - x = np.array( + X = np.array( [ [-2.2, -1.9], [-1.9, -2.1], @@ -76,11 +76,11 @@ source samples in ``y``, and a covariate vector identifying the target domain. source_labels = np.array([0, 0, 1, 1]) domains = np.array([0, 0, 0, 0, 1, 1, 1, 1]) - x_target = x[domains == 1] + X_target = X[domains == 1] clf = ARSVM() - clf.fit(x, source_labels, covariates=domains, target_covariate=1) - y_pred = clf.predict(x_target) + clf.fit(X, source_labels, covariates=domains, target_covariate=1) + y_pred = clf.predict(X_target) Train a Manifold-Regularized Classifier --------------------------------------- @@ -93,12 +93,12 @@ samples. The labels array may contain only the labeled source examples. import numpy as np from kalelinear.estimator import LapSVM - x_source = np.array([[-2.0, -1.8], [-1.8, -2.1], [1.9, 1.7], [2.1, 2.0]]) + X_source = np.array([[-2.0, -1.8], [-1.8, -2.1], [1.9, 1.7], [2.1, 2.0]]) ys = np.array([0, 0, 1, 1]) - x_target = np.array([[-1.4, -1.2], [-1.2, -1.1], [1.2, 1.1], [1.4, 1.3]]) + X_target = np.array([[-1.4, -1.2], [-1.2, -1.1], [1.2, 1.1], [1.4, 1.3]]) - x_train = np.vstack((x_source, x_target)) + X_train = np.vstack((X_source, X_target)) clf = LapSVM(kernel="linear") - clf.fit(x_train, ys) - y_pred = clf.predict(x_target) + clf.fit(X_train, ys) + y_pred = clf.predict(X_target) diff --git a/docs/source/usage.rst b/docs/source/usage.rst index 10e0566..29f1116 100644 --- a/docs/source/usage.rst +++ b/docs/source/usage.rst @@ -15,7 +15,7 @@ label is the target domain. import numpy as np from kalelinear.transformer import TCA - x = np.array( + X = np.array( [ [-2.0, -1.8], [-1.8, -2.1], @@ -30,7 +30,7 @@ label is the target domain. domain_labels = np.array([0, 0, 0, 0, 1, 1, 1, 1]) transformer = TCA(n_components=2) - z = transformer.fit_transform(x, covariates=domain_labels, target_covariate=1) + z = transformer.fit_transform(X, covariates=domain_labels, target_covariate=1) z_source = z[domain_labels == 0] z_target = z[domain_labels == 1] @@ -38,7 +38,7 @@ label is the target domain. Domain Adaptation Estimators ---------------------------- -ARSVM and ARRLS use all source and target samples in ``x``, labels for the +ARSVM and ARRLS use all source and target samples in ``X``, labels for the source samples, and covariates that mark each sample's domain. .. code-block:: python @@ -46,7 +46,7 @@ source samples, and covariates that mark each sample's domain. import numpy as np from kalelinear.estimator import ARSVM - x = np.array( + X = np.array( [ [-2.2, -1.9], [-1.9, -2.1], @@ -60,8 +60,8 @@ source samples, and covariates that mark each sample's domain. ) source_labels = np.array([0, 0, 1, 1]) domains = np.array([0, 0, 0, 0, 1, 1, 1, 1]) - x_target = x[domains == 1] + X_target = X[domains == 1] clf = ARSVM() - clf.fit(x, source_labels, covariates=domains, target_covariate=1) - y_pred = clf.predict(x_target) + clf.fit(X, source_labels, covariates=domains, target_covariate=1) + y_pred = clf.predict(X_target) diff --git a/kalelinear/estimator/_artl.py b/kalelinear/estimator/_artl.py index f066d45..23c532b 100644 --- a/kalelinear/estimator/_artl.py +++ b/kalelinear/estimator/_artl.py @@ -107,7 +107,7 @@ def __init__( kernel="linear", lambda_=1.0, gamma_=0.0, - k_neighbour=5, + k_neighbors=5, solver="osqp", manifold_metric="cosine", knn_mode="distance", @@ -125,7 +125,7 @@ def __init__( MMD regulisation param, by default 1.0 gamma_ : float, optional manifold regulisation param, by default 0.0 - k_neighbour : int, optional + k_neighbors : int, optional number of nearest numbers for each sample in manifold regularisation, by default 5 solver : str, optional @@ -148,7 +148,7 @@ def __init__( self.C = C self.gamma_ = gamma_ self.solver = solver - self.k_neighbour = k_neighbour + self.k_neighbors = k_neighbors # self.alpha = None self.knn_mode = knn_mode self.manifold_metric = manifold_metric @@ -191,7 +191,7 @@ def fit(self, X, y, covariates=None, target_covariate=None, unlabeled_value=None y_ = self._lb.fit_transform(y) if self.gamma_ != 0: - lap_mat = lap_norm(X, n_neighbour=self.k_neighbour, mode=self.knn_mode) + lap_mat = lap_norm(X, n_neighbors=self.k_neighbors, mode=self.knn_mode) Q_ = unit_matrix + multi_dot([(self.lambda_ * M + self.gamma_ * lap_mat), x_kernel_matrix]) else: Q_ = unit_matrix + multi_dot([(self.lambda_ * M), x_kernel_matrix]) @@ -281,7 +281,7 @@ def __init__( lambda_=1.0, gamma_=0.0, sigma_=1.0, - k_neighbour=5, + k_neighbors=5, manifold_metric="cosine", knn_mode="distance", **kwargs, @@ -298,7 +298,7 @@ def __init__( manifold regularisation param, by default 0.0 sigma_ : float, optional l2 regularisation param, by default 1.0 - k_neighbour : int, optional + k_neighbors : int, optional number of nearest numbers for each sample in manifold regularisation, by default 5 manifold_metric : str, optional @@ -318,7 +318,7 @@ def __init__( self.lambda_ = lambda_ self.gamma_ = gamma_ self.sigma_ = sigma_ - self.k_neighbour = k_neighbour + self.k_neighbors = k_neighbors # self.coef_ = None self.knn_mode = knn_mode self.manifold_metric = manifold_metric @@ -360,7 +360,7 @@ def fit(self, X, y, covariates=None, target_covariate=None, unlabeled_value=None J[:nl, :nl] = np.eye(nl) if self.gamma_ != 0: - lap_mat = lap_norm(X, n_neighbour=self.k_neighbour, metric=self.manifold_metric, mode=self.knn_mode) + lap_mat = lap_norm(X, n_neighbors=self.k_neighbors, metric=self.manifold_metric, mode=self.knn_mode) Q_ = np.dot((J + self.lambda_ * M + self.gamma_ * lap_mat), x_kernel_matrix) + self.sigma_ * unit_matrix else: Q_ = np.dot((J + self.lambda_ * M), x_kernel_matrix) + self.sigma_ * unit_matrix diff --git a/kalelinear/estimator/_coir.py b/kalelinear/estimator/_coir.py index eb399fa..b6af2c4 100644 --- a/kalelinear/estimator/_coir.py +++ b/kalelinear/estimator/_coir.py @@ -46,7 +46,7 @@ def __init__( kernel="linear", lambda_=1.0, mu=0.0, - k_neighbour=3, + k_neighbors=3, manifold_metric="cosine", knn_mode="distance", solver="osqp", @@ -65,7 +65,7 @@ def __init__( param for covariate (side information) independence regularisation, by default 1 mu : float, optional param for manifold regularisation, by default 0 - k_neighbour : int, optional + k_neighbors : int, optional number of nearest numbers for each sample in manifold regularisation, by default 3 manifold_metric : str, optional @@ -94,7 +94,7 @@ def __init__( # self.support_vectors_ = None # self.n_support_ = None self.manifold_metric = manifold_metric - self.k_neighbour = k_neighbour + self.k_neighbors = k_neighbors self.knn_mode = knn_mode self.covariate_encoder = covariate_encoder self._lb = LabelBinarizer(pos_label=1, neg_label=-1) @@ -130,7 +130,7 @@ def fit(self, X, y, covariates=None): Q_ = unit_matrix.copy() if self.mu != 0: - lap_mat = lap_norm(X, n_neighbour=self.k_neighbour, metric=self.manifold_metric, mode=self.knn_mode) + lap_mat = lap_norm(X, n_neighbors=self.k_neighbors, metric=self.manifold_metric, mode=self.knn_mode) Q_ += np.dot( self.lambda_ / np.square(n - 1) * covariate_hsic_matrix + self.mu / np.square(n) * lap_mat, x_kernel_matrix, @@ -316,7 +316,7 @@ def fit(self, X, y, covariates=None): J[:nl, :nl] = np.eye(nl) if self.mu != 0: - lap_mat = lap_norm(X, n_neighbour=self.k, mode=self.knn_mode, metric=self.manifold_metric) + lap_mat = lap_norm(X, n_neighbors=self.k, mode=self.knn_mode, metric=self.manifold_metric) Q_ = self.sigma_ * unit_matrix + np.dot( J + self.lambda_ / np.square(n - 1) * covariate_hsic_matrix + self.mu / np.square(n) * lap_mat, x_kernel_matrix, diff --git a/kalelinear/estimator/_gsda.py b/kalelinear/estimator/_gsda.py index 293588f..5337e7a 100644 --- a/kalelinear/estimator/_gsda.py +++ b/kalelinear/estimator/_gsda.py @@ -307,7 +307,7 @@ def _lbfgs_solver(self, X, y, groups, target_idx=None): z += delta_theta[i] * (alphas[len(alphas) - 1 - i] - beta_i) else: z = np.eye(self.theta_.shape[0]) @ q - # Line search and update x + # Line search and update X # Implement a line search algorithm to find an appropriate step size # step_size = 1 # Placeholder # Update memory @@ -379,20 +379,20 @@ def _gd_solver(self, X, y, groups, target_idx=None): return self def compute_gsda_gradient(self, X, y, groups, target_idx=None): - n_sample = X.shape[0] + n_samples = X.shape[0] n_tgt = y.shape[0] if target_idx is None: - x_tgt = X[:n_tgt] + X_tgt = X[:n_tgt] else: - x_tgt = X[target_idx] + X_tgt = X[target_idx] - y_hat = expit(x_tgt @ self.theta_) + y_hat = expit(X_tgt @ self.theta_) # n_feature = X.shape[1] _simple_hsic = simple_hsic_grad_term(self.theta_, X, groups) - hsic_proba = expit(np.dot(self.theta_, _simple_hsic) / np.square(n_sample - 1)) - grad_hsic = (hsic_proba - 1) * _simple_hsic / np.square(n_sample - 1) + hsic_proba = expit(np.dot(self.theta_, _simple_hsic) / np.square(n_samples - 1)) + grad_hsic = (hsic_proba - 1) * _simple_hsic / np.square(n_samples - 1) - delta_grad = (x_tgt.T @ (y_hat - y)) / n_tgt + delta_grad = (X_tgt.T @ (y_hat - y)) / n_tgt if self.regularization is not None: delta_grad += self.theta_ * self.alpha delta_grad += self.lambda_ * grad_hsic diff --git a/kalelinear/estimator/_manifold_learn.py b/kalelinear/estimator/_manifold_learn.py index 62bf880..1410148 100644 --- a/kalelinear/estimator/_manifold_learn.py +++ b/kalelinear/estimator/_manifold_learn.py @@ -24,7 +24,7 @@ def __init__( kernel="linear", gamma_=1.0, solver="osqp", - k_neighbour=3, + k_neighbors=3, manifold_metric="cosine", knn_mode="distance", **kwargs, @@ -41,7 +41,7 @@ def __init__( param for manifold regularisation, by default 1.0 solver : str, optional quadratic programming solver, [cvxopt, osqp], by default 'osqp' - k_neighbour : int, optional + k_neighbors : int, optional number of nearest numbers for each sample in manifold regularisation, by default 3 manifold_metric : str, optional @@ -62,7 +62,7 @@ def __init__( self.solver = solver self.kwargs = kwargs self.manifold_metric = manifold_metric - self.k_neighbour = k_neighbour + self.k_neighbors = k_neighbors self.knn_mode = knn_mode self._lb = LabelBinarizer(pos_label=1, neg_label=-1) @@ -87,7 +87,7 @@ def fit(self, X, y): if self.gamma_ == 0: Q_ = centering_matrix else: - lap_mat = lap_norm(X, n_neighbour=self.k_neighbour, mode=self.knn_mode) + lap_mat = lap_norm(X, n_neighbors=self.k_neighbors, mode=self.knn_mode) Q_ = centering_matrix + self.gamma_ * np.dot(lap_mat, x_kernel_matrix) y_ = self._lb.fit_transform(y) @@ -171,7 +171,7 @@ def __init__( kernel="linear", gamma_=1.0, sigma_=1.0, - k_neighbour=5, + k_neighbors=5, manifold_metric="cosine", knn_mode="distance", **kwargs, @@ -186,7 +186,7 @@ def __init__( manifold regularisation param, by default 1.0 sigma_ : float, optional l2 regularisation param, by default 1.0 - k_neighbour : int, optional + k_neighbors : int, optional number of nearest numbers for each sample in manifold regularisation, by default 5 manifold_metric : str, optional @@ -205,7 +205,7 @@ def __init__( self.kernel = kernel self.gamma_ = gamma_ self.sigma_ = sigma_ - self.k_neighbour = k_neighbour + self.k_neighbors = k_neighbors # self.coef_ = None self.knn_mode = knn_mode self.manifold_metric = manifold_metric @@ -235,7 +235,7 @@ def fit(self, X, y): J[:nl, :nl] = np.eye(nl) if self.gamma_ != 0: - lap_mat = lap_norm(X, n_neighbour=self.k_neighbour, metric=self.manifold_metric, mode=self.knn_mode) + lap_mat = lap_norm(X, n_neighbors=self.k_neighbors, metric=self.manifold_metric, mode=self.knn_mode) Q_ = np.dot((J + self.gamma_ * lap_mat), x_kernel_matrix) + self.sigma_ * unit_matrix else: Q_ = np.dot(J, x_kernel_matrix) + self.sigma_ * centering_matrix diff --git a/kalelinear/estimator/base.py b/kalelinear/estimator/base.py index c0efaa5..f27c321 100644 --- a/kalelinear/estimator/base.py +++ b/kalelinear/estimator/base.py @@ -20,7 +20,7 @@ class BaseKaleEstimator(BaseEstimator, ClassifierMixin): def __init__( self, kernel="linear", - k_neighbour=5, + k_neighbors=5, manifold_metric="cosine", knn_mode="distance", pos_label=1, @@ -29,7 +29,7 @@ def __init__( ): super().__init__() self.kernel = kernel - self.k_neighbour = k_neighbour + self.k_neighbors = k_neighbors self.manifold_metric = manifold_metric self.knn_mode = knn_mode self.coef_ = None @@ -82,6 +82,15 @@ def _quadprog(cls, P, y, C, solver="osqp"): q = -1 * np.ones(n_labeled) upper_bound = C / n_labeled + # The semi-dual Hessian is only positive-semidefinite up to rounding + # (manifold/MMD terms such as ``L @ K`` are not PSD in general), which + # makes convex QP solvers such as osqp fail on some datasets. Project + # ``P`` onto the PSD cone so the QP is always convex and solvable. + P = P.astype(np.float64) + P = 0.5 * (P + P.T) + eigenvalues, eigenvectors = np.linalg.eigh(P) + P = (eigenvectors * np.clip(eigenvalues, 0, None)) @ eigenvectors.T + if solver == "cvxopt": G = np.zeros((2 * n_labeled, n_labeled)) G[:n_labeled, :] = -1 * np.eye(n_labeled) diff --git a/kalelinear/pipeline/__init__.py b/kalelinear/pipeline/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/kalelinear/pipeline/mpca_trainer.py b/kalelinear/pipeline/mpca_trainer.py new file mode 100644 index 0000000..cd1ce2d --- /dev/null +++ b/kalelinear/pipeline/mpca_trainer.py @@ -0,0 +1,222 @@ +# ============================================================================= +# Author: Shuo Zhou, shuo.zhou@sheffield.ac.uk +# Haiping Lu, h.lu@sheffield.ac.uk or hplu@ieee.org +# ============================================================================= + +"""Implementation of MPCA->Feature Selection->Linear SVM/LogisticRegression Pipeline + +References: + [1] Swift, A. J., Lu, H., Uthoff, J., Garg, P., Cogliano, M., Taylor, J., ... & Kiely, D. G. (2020). A machine + learning cardiac magnetic resonance approach to extract disease features and automate pulmonary arterial + hypertension diagnosis. European Heart Journal-Cardiovascular Imaging. + [2] Song, X., Meng, L., Shi, Q., & Lu, H. (2015, October). Learning tensor-based features for whole-brain fMRI + classification. In International Conference on Medical Image Computing and Computer-Assisted Intervention + (pp. 613-620). Springer, Cham. + [3] Lu, H., Plataniotis, K. N., & Venetsanopoulos, A. N. (2008). MPCA: Multilinear principal component analysis of + tensor objects. IEEE Transactions on Neural Networks, 19(1), 18-39. +""" + +import logging + +import numpy as np +from sklearn.base import BaseEstimator, ClassifierMixin +from sklearn.feature_selection import f_classif +from sklearn.linear_model import LogisticRegression +from sklearn.model_selection import GridSearchCV +from sklearn.svm import LinearSVC, SVC +from sklearn.utils.validation import check_is_fitted + +from kalelinear.transformer import MPCA + +param_c_grids = list(np.logspace(-4, 2, 7)) +classifiers = { + "svc": [SVC, {"kernel": ["linear"], "C": param_c_grids, "max_iter": [50000]}], + "linear_svc": [LinearSVC, {"C": param_c_grids}], + "lr": [LogisticRegression, {"C": param_c_grids}], +} + +# k-fold cross-validation used for grid search, i.e. searching for optimal value of C +default_search_params = {"cv": 5} +default_mpca_params = {"explained_variance_ratio": 0.97, "vectorize": True} + + +class MPCATrainer(BaseEstimator, ClassifierMixin): + """Trainer of pipeline: MPCA->Feature selection->Classifier + + Args: + classifier (str, optional): Available classifier options: {"svc", "linear_svc", "lr"}, where "svc" trains a + support vector classifier, supports both linear and non-linear kernels, optimizes with library "libsvm"; + "linear_svc" trains a support vector classifier with linear kernel only, and optimizes with library + "liblinear", which suppose to be faster and better in handling large number of samples; and "lr" trains + a classifier with logistic regression. Defaults to "svc". + classifier_params (dict, optional): Parameters of classifier. Defaults to 'auto'. + classifier_param_grid (dict, optional): Grids for searching the optimal hyper-parameters. Works only when + classifier_params == "auto". Defaults to None by searching from the following hyper-parameter values: + 1. svc, {"kernel": ["linear"], "C": [0.0001, 0.001, 0.01, 0.1, 1, 10, 100], "max_iter": [50000]}, + 2. linear_svc, {"C": [0.0001, 0.001, 0.01, 0.1, 1, 10, 100]}, + 3. lr, {"C": [0.0001, 0.001, 0.01, 0.1, 1, 10, 100]} + mpca_params (dict, optional): Parameters of MPCA, e.g., {"explained_variance_ratio": 0.8}. Defaults to None, + i.e., using the default parameters + (https://kalelinear.readthedocs.io/en/latest/api_transformers.html#kalelinear.transformer.MPCA). + n_features (int, optional): Number of features for feature selection. Defaults to None, i.e., all features + after dimension reduction will be used. + search_params (dict, optional): Parameters of grid search, for more detail please see + https://scikit-learn.org/stable/modules/grid_search.html#grid-search. Defaults to None, i.e., using the + default params: {"cv": 5}. + """ + + def __init__( + self, + classifier="svc", + classifier_params="auto", + classifier_param_grid=None, + mpca_params=None, + n_features=None, + search_params=None, + ): + if classifier not in ["svc", "linear_svc", "lr"]: + error_msg = "Valid classifier should be 'svc', 'linear_svc', or 'lr', but given %s" % classifier + logging.error(error_msg) + raise ValueError(error_msg) + + self.classifier = classifier + # init mpca object + if mpca_params is None: + self.mpca_params = dict(default_mpca_params) + else: + self.mpca_params = dict(mpca_params) + self.mpca = MPCA(**self.mpca_params) + # init feature selection parameters + self.n_features = n_features + # init classifier object + if search_params is None: + self.search_params = dict(default_search_params) + else: + self.search_params = dict(search_params) + self.classifier_param_grid = classifier_param_grid + + self.auto_classifier_param = False + if classifier_params == "auto": + self.auto_classifier_param = True + if self.classifier_param_grid is None: + self.classifier_param_grid = { + param_name: list(values) for param_name, values in classifiers[classifier][1].items() + } + self.grid_search = GridSearchCV( + classifiers[classifier][0](), param_grid=self.classifier_param_grid, **self.search_params + ) + self.clf = None + elif isinstance(classifier_params, dict): + self.clf = classifiers[classifier][0](**classifier_params) + else: + error_msg = "Invalid classifier parameter type" + logging.error(error_msg) + raise ValueError(error_msg) + + if isinstance(classifier_params, dict): + self.classifier_params = dict(classifier_params) + self.clf = classifiers[classifier][0](**classifier_params) + elif classifier_params == "auto": + self.auto_classifier_param = True + self.classifier_params = "auto" + else: + error_msg = "Invalid classifier parameter type" + logging.error(error_msg) + raise ValueError(error_msg) + + def fit(self, X, y): + """Fit a pipeline with the given data X and labels y + + Args: + X (array-like tensor): input data, shape (n_samples, I_1, I_2, ..., I_N) + y (array-like): data labels, shape (n_samples, ) + + Returns: + self + """ + # fit mpca + self.mpca.fit(X) + self.mpca.set_params(**{"vectorize": True}) + X_transformed = self.mpca.transform(X) + + # feature selection + if self.n_features is None: + # MPCA.transform(vectorize=True) already returns features ordered by variance. + self.feature_order_ = np.arange(X_transformed.shape[1]) + self.n_features_ = X_transformed.shape[1] + else: + f_score, _ = f_classif(X_transformed, y) + self.feature_order_ = np.argsort(-f_score) + self.n_features_ = min(self.n_features, X_transformed.shape[1]) + X_transformed = X_transformed[:, self.feature_order_][:, : self.n_features_] + + # fit classifier + if self.auto_classifier_param: + param_grid = {name: list(values) for name, values in self.classifier_param_grid.items()} + extra_c = 1 / X.shape[0] + if extra_c not in param_grid["C"]: + param_grid["C"].append(extra_c) + self.grid_search = GridSearchCV( + classifiers[self.classifier][0](), param_grid=param_grid, **self.search_params + ) + self.grid_search.fit(X_transformed, y) + self.clf = self.grid_search.best_estimator_ + if self.classifier == "svc": + self.clf.set_params(**{"probability": True}) + + self.clf.fit(X_transformed, y) + return self + + def predict(self, X): + """Predict the labels for the given data X + + Args: + X (array-like tensor): input data, shape (n_samples, I_1, I_2, ..., I_N) + + Returns: + array-like: Predicted labels, shape (n_samples, ) + """ + features = self._extract_feature(X) + return self.clf.predict(features) + + def decision_function(self, X): + """Decision scores of each class for the given data X + + Args: + X (array-like tensor): input data, shape (n_samples, I_1, I_2, ..., I_N) + + Returns: + array-like: decision scores, shape (n_samples,) for binary case, else (n_samples, n_classes) + """ + features = self._extract_feature(X) + return self.clf.decision_function(features) + + def predict_proba(self, X): + """Probability of each class for the given data X. Not supported by "linear_svc". + + Args: + X (array-like tensor): input data, shape (n_samples, I_1, I_2, ..., I_N) + + Returns: + array-like: probabilities, shape (n_samples, n_classes) + """ + if self.classifier == "linear_svc": + error_msg = "Linear SVC does not support computing probability." + logging.error(error_msg) + raise ValueError(error_msg) + features = self._extract_feature(X) + return self.clf.predict_proba(features) + + def _extract_feature(self, X): + """Extracting features for the given data X with MPCA->Feature selection + + Args: + X (array-like tensor): input data, shape (n_samples, I_1, I_2, ..., I_N) + + Returns: + array-like: n_new, shape (n_samples, n_features) + """ + check_is_fitted(self) + X_transformed = self.mpca.transform(X) + + return X_transformed[:, self.feature_order_][:, : self.n_features_] diff --git a/kalelinear/transformer/_base.py b/kalelinear/transformer/_base.py index 2867804..f5de7d7 100644 --- a/kalelinear/transformer/_base.py +++ b/kalelinear/transformer/_base.py @@ -106,6 +106,17 @@ def _eigendecompose( """Compute eigenpairs for a kernel matrix or a generalized eigenproblem.""" a, b = _get_eigenproblem_matrices(eigenproblem) + # Generalized eigenproblems require a positive-definite ``b``. Constraint + # matrices such as centered kernel matrices are only positive-semidefinite + # up to floating-point rounding, which makes scipy's ``eigh(a, b)`` fail + # with a non-positive-definite error. Regularize ``b`` with a small ridge + # relative to its scale; the shift is negligible for the retained + # components. + if b is not None: + b = 0.5 * (b + b.T) + ridge = 10 * np.finfo(b.dtype).eps * max(1.0, np.max(np.abs(b))) * b.shape[0] + b = b + ridge * np.eye(b.shape[0], dtype=b.dtype) + n_components = _check_n_components((a, b), n_components) solver = _check_solver((a, b), n_components, solver, eigenvalue_order) @@ -307,15 +318,15 @@ def __init__( self.fit_inverse_transform = fit_inverse_transform self.augment = augment - def _fit_inverse_transform(self, x_transformed, X): + def _fit_inverse_transform(self, X_transformed, X): if hasattr(X, "tocsr"): raise NotImplementedError("Inverse transform not implemented for sparse matrices!") - n_samples = x_transformed.shape[0] - x_transformed_kernel_matrix = self._get_kernel(x_transformed) + n_samples = X_transformed.shape[0] + x_transformed_kernel_matrix = self._get_kernel(X_transformed) x_transformed_kernel_matrix.flat[:: n_samples + 1] += self.alpha self.dual_coef_ = la.solve(x_transformed_kernel_matrix, X, assume_a="pos", overwrite_a=True) - self.x_transformed_fit_ = x_transformed + self.X_transformed_fit_ = X_transformed @property def _n_features_out(self): @@ -350,7 +361,7 @@ def orig_coef_(self): w = self.eigenvectors_ if self.scale_components: w = _scale_eigenvectors(self.eigenvalues_, w) - return safe_sparse_dot(w.T, self.x_fit_) + return safe_sparse_dot(w.T, self.X_fit_) def _requires_covariates(self): return False @@ -546,8 +557,8 @@ def fit(self, X, y=None, covariates=None, **fit_params): elif hasattr(self, "_factor_validator"): delattr(self, "_factor_validator") - self.x_fit_ = self._augment_data(context.X_fit, context.covariates_fit) - self.x_fit_raw_ = context.X_fit + self.X_fit_ = self._augment_data(context.X_fit, context.covariates_fit) + self.X_fit_raw_ = context.X_fit self.covariates_input_ = raw_covariates self.covariates_fit_ = context.covariates_fit self.fit_context_ = context @@ -556,15 +567,15 @@ def fit(self, X, y=None, covariates=None, **fit_params): self.gamma_ = 1 / _num_features(X) if self.gamma is None else self.gamma self._centerer = KernelCenterer() - x_fit_kernel_matrix = self._get_kernel(self.x_fit_) - x_fit_kernel_matrix = self._centerer.fit_transform(x_fit_kernel_matrix) + X_fit_kernel_matrix = self._get_kernel(self.X_fit_) + X_fit_kernel_matrix = self._centerer.fit_transform(X_fit_kernel_matrix) - eigenproblem = self._make_eigenproblem(x_fit_kernel_matrix, context) + eigenproblem = self._make_eigenproblem(X_fit_kernel_matrix, context) self._fit_transform_in_place(eigenproblem) if self.fit_inverse_transform: - x_transformed = self.transform(context.X_input, covariates=context.covariates_input) - self._fit_inverse_transform(x_transformed, context.X_input) + X_transformed = self.transform(context.X_input, covariates=context.covariates_input) + self._fit_inverse_transform(X_transformed, context.X_input) return self @@ -593,21 +604,21 @@ def transform(self, X, covariates=None): covariates = self._factor_validator.transform(covariates) X_query = self._augment_data(X, covariates) - x_fit_kernel_matrix = self._get_kernel(X_query, self.x_fit_) - x_fit_kernel_matrix = self._centerer.transform(x_fit_kernel_matrix) + X_fit_kernel_matrix = self._get_kernel(X_query, self.X_fit_) + X_fit_kernel_matrix = self._centerer.transform(X_fit_kernel_matrix) w = self.eigenvectors_ if self.scale_components: w = _scale_eigenvectors(self.eigenvalues_, w) - z = safe_sparse_dot(x_fit_kernel_matrix, w) + z = safe_sparse_dot(X_fit_kernel_matrix, w) if self.augment == "post": z = np.hstack((z, covariates)) return z - def inverse_transform(self, z): + def inverse_transform(self, X): check_is_fitted(self) if not self.fit_inverse_transform: raise NotFittedError( @@ -616,7 +627,7 @@ def inverse_transform(self, z): "the inverse transform is not available." ) - k_z = self._get_kernel(z, self.x_transformed_fit_) + k_z = self._get_kernel(X, self.X_transformed_fit_) return safe_sparse_dot(k_z, self.dual_coef_) def fit_transform(self, X, y=None, covariates=None, **fit_params): diff --git a/kalelinear/transformer/_mpca.py b/kalelinear/transformer/_mpca.py index e75befe..07373eb 100644 --- a/kalelinear/transformer/_mpca.py +++ b/kalelinear/transformer/_mpca.py @@ -12,11 +12,37 @@ import warnings import numpy as np -import scipy.linalg as la +from numpy.linalg import eigvalsh +from scipy.linalg import eigh from sklearn.base import BaseEstimator, TransformerMixin from tensorly.base import fold, unfold from tensorly.tenalg import multi_mode_dot +_CHUNK_ELEMS = 4000000 # target number of float64 elements per processed chunk (~32 MB) + + +def _chunk_size(n_samples, elems_per_sample, chunk_elems=_CHUNK_ELEMS): + """Return the largest chunk keeping per-chunk element count under ``chunk_elems``. + + Parameters + ---------- + n_samples : int + Number of samples in the data. + elems_per_sample : int + Number of array elements contributed by a single sample (e.g. the + number of features of the input or unfolded tensor). + chunk_elems : int, default=4000000 + Memory budget in number of elements for one processed chunk. + + Returns + ------- + chunk_size : int + Number of samples to process at a time. + """ + if elems_per_sample <= 0: + return n_samples + return min(n_samples, max(1, chunk_elems // elems_per_sample)) + def _check_n_dim(X, n_dims): """Validate the number of dimensions. @@ -81,7 +107,7 @@ class MPCA(BaseEstimator, TransformerMixin): Parameters ---------- - var_ratio : float, default=0.97 + explained_variance_ratio : float, default=0.97 Target cumulative explained variance ratio per mode. max_iter : int, default=1 Maximum number of alternating optimization iterations. @@ -89,19 +115,28 @@ class MPCA(BaseEstimator, TransformerMixin): If ``True``, output projected tensors as vectors. n_components : int, optional Number of output features when ``vectorize=True``. + output_shape : tuple of int, optional + Number of components to keep per mode. If given, it overrides + ``explained_variance_ratio`` and the output dimensions are set exactly. + Default is None, i.e. the output shape is derived from + ``explained_variance_ratio``. Attributes ---------- - proj_mats : list of ndarray + proj_mats_ : list of ndarray Transposed projection matrices with shapes ``(P_i, I_i)``. - idx_order : ndarray + idx_order_ : ndarray Feature ranking indices by descending projected variance. mean_ : ndarray Per-feature empirical mean of the training data. - shape_in : tuple + input_shape_ : tuple Input per-sample tensor shape. - shape_out : tuple - Output per-sample tensor shape after projection. + output_shape_ : tuple + Output per-sample tensor shape after projection. Equals + ``output_shape`` when given, otherwise determined by + ``explained_variance_ratio``. + explained_variance_ratio_ : tuple of float + Achieved cumulative explained variance ratio per mode after fitting. References ---------- @@ -113,35 +148,46 @@ class MPCA(BaseEstimator, TransformerMixin): -------- >>> import numpy as np >>> from kalelinear.transformer import MPCA - >>> x = np.random.random((40, 20, 25, 20)) - >>> x.shape + >>> X = np.random.random((40, 20, 25, 20)) + >>> X.shape (40, 20, 25, 20) >>> mpca = MPCA() - >>> x_projected = mpca.fit_transform(x) - >>> x_projected.shape + >>> X_projected = mpca.fit_transform(X) + >>> X_projected.shape (40, 18, 23, 18) - >>> x_projected = mpca.transform(x) - >>> x_projected.shape + >>> X_projected = mpca.transform(X) + >>> X_projected.shape (40, 7452) - >>> x_projected = mpca.transform(x) - >>> x_projected.shape + >>> X_projected = mpca.transform(X, vectorize=True) + >>> X_projected.shape (40, 50) - >>> x_rec = mpca.inverse_transform(x_projected) - >>> x_rec.shape + >>> X_reconstructed = mpca.inverse_transform(X_projected) + >>> X_reconstructed.shape (40, 20, 25, 20) """ - def __init__(self, var_ratio=0.97, max_iter=1, vectorize=False, n_components=None): - self.var_ratio = var_ratio + def __init__( + self, explained_variance_ratio=0.97, max_iter=1, vectorize=False, n_components=None, output_shape=None + ): + self.explained_variance_ratio = explained_variance_ratio if max_iter > 0 and isinstance(max_iter, int): self.max_iter = max_iter else: msg = "Number of max iterations must be a positive integer but given %s" % max_iter logging.error(msg) raise ValueError(msg) - self.proj_mats = [] self.vectorize = vectorize self.n_components = n_components + if output_shape is None: + self.output_shape = None + elif isinstance(output_shape, (tuple, list)) and all( + isinstance(v, (int, np.integer)) and v > 0 for v in output_shape + ): + self.output_shape = tuple(int(v) for v in output_shape) + else: + msg = "output_shape must be None or a sequence of positive integers but given %s" % (output_shape,) + logging.error(msg) + raise ValueError(msg) def fit(self, X, y=None): """Fit MPCA to tensor data. @@ -176,88 +222,162 @@ def _fit(self, X): """ shape_ = X.shape # shape of input data + n_samples = shape_[0] n_dims = X.ndim - self.shape_in = shape_[1:] - self.mean_ = np.mean(X, axis=0) - X = X - self.mean_ + if n_samples < 2: + error_msg = "MPCA requires at least 2 samples to fit." + logging.error(error_msg) + raise ValueError(error_msg) - # init - shape_out = () - proj_mats = [] + self.input_shape_ = shape_[1:] - # get the output tensor shape based on the cumulative distribution of eigen values for each mode - for i in range(1, n_dims): - mode_data_mat = unfold(X, mode=i) - singular_vec_left, singular_val, singular_vec_right = la.svd(mode_data_mat, full_matrices=False) - eig_values = np.square(singular_val) - idx_sorted = (-1 * eig_values).argsort() - cum = eig_values[idx_sorted] - tot_var = np.sum(cum) - - for j in range(1, cum.shape[0] + 1): - if np.sum(cum[:j]) / tot_var > self.var_ratio: - shape_out += (j,) - break - proj_mats.append(singular_vec_left[:, idx_sorted][:, : shape_out[i - 1]].T) - - for i_iter in range(self.max_iter): - for i in range(1, n_dims): # ith mode - x_projected = multi_mode_dot( - X, - [proj_mats[m] for m in range(n_dims - 1) if m != i - 1], - modes=[m for m in range(1, n_dims) if m != i], - ) - mode_data_mat = unfold(x_projected, i) - - singular_vec_left, singular_val, singular_vec_right = la.svd(mode_data_mat, full_matrices=False) - eig_values = np.square(singular_val) - idx_sorted = (-1 * eig_values).argsort() - proj_mats[i - 1] = (singular_vec_left[:, idx_sorted][:, : shape_out[i - 1]]).T + # Samples are processed in chunks so that a centered copy of the full + # data never has to be materialized and disk-backed inputs (e.g. a + # memory-mapped array or a loader over file paths) are read one chunk + # at a time. + n_features_in = int(np.prod(shape_[1:])) + chunk_size = _chunk_size(n_samples, n_features_in) - x_projected = multi_mode_dot(X, proj_mats, modes=[m for m in range(1, n_dims)]) - x_proj_unfold = unfold(x_projected, mode=0) # unfold the tensor projection to shape (n_samples, n_features) - # x_proj_cov = np.diag(np.dot(x_proj_unfold.T, x_proj_unfold)) # covariance of unfolded features - x_proj_cov = np.sum(np.multiply(x_proj_unfold.T, x_proj_unfold.T), axis=1) # memory saving computing covariance - idx_order = (-1 * x_proj_cov).argsort() + mean_acc = np.zeros(shape_[1:], dtype=np.float64) + for start in range(0, n_samples, chunk_size): + mean_acc += X[start : start + chunk_size].sum(axis=0, dtype=np.float64) + self.mean_ = mean_acc / n_samples - self.proj_mats = proj_mats - self.idx_order = idx_order - self.shape_out = shape_out - self.n_dims = n_dims + # init: accumulate per-mode covariance over chunked unfoldings + covariance_matrices = {} + for i in range(1, n_dims): + mode_cov = np.zeros((shape_[i], shape_[i])) + for start in range(0, n_samples, chunk_size): + batch = X[start : start + chunk_size] - self.mean_ + batch_unfold = unfold(batch, mode=i) + mode_cov += batch_unfold @ batch_unfold.T + covariance_matrices[i] = mode_cov + + # get the output tensor shape: either user-specified or derived from the + # cumulative distribution of eigen values for each mode + if self.output_shape is not None: + output_shape = self.output_shape + if len(output_shape) != n_dims - 1: + error_msg = "output_shape must have length %s (one entry per mode) but has %s" % ( + n_dims - 1, + len(output_shape), + ) + logging.error(error_msg) + raise ValueError(error_msg) + for mode_i, n_comp in enumerate(output_shape, start=1): + if n_comp > shape_[mode_i]: + error_msg = "output_shape entry %s must not exceed the input size %s of mode %s but is %s" % ( + mode_i - 1, + shape_[mode_i], + mode_i, + n_comp, + ) + logging.error(error_msg) + raise ValueError(error_msg) + else: + output_shape = () + + proj_matrices = [] + explained_variance_ratios = [] + for mode_i in range(1, n_dims): + if self.output_shape is None: + eigenvalues = eigvalsh(covariance_matrices[mode_i]) + idx_sorted = np.argsort(eigenvalues)[::-1] + cum = eigenvalues[idx_sorted] + tot_var = np.sum(cum) + + cum_var = np.cumsum(cum) + mode_n_components = min( + int(np.searchsorted(cum_var, self.explained_variance_ratio * tot_var, side="right")) + 1, + shape_[mode_i], + ) + output_shape += (mode_n_components,) + else: + mode_n_components = output_shape[mode_i - 1] + + # Only the top mode_n_components eigenvectors are needed; scipy eigh + # returns them in ascending eigenvalue order for the requested range. + subset_eigenvalues, eigenvectors = eigh( + covariance_matrices[mode_i], subset_by_index=[shape_[mode_i] - mode_n_components, shape_[mode_i] - 1] + ) + tot_var = np.trace(covariance_matrices[mode_i]) + explained_variance_ratio = subset_eigenvalues.sum() / tot_var if tot_var > 0 else 0.0 + explained_variance_ratios.append(explained_variance_ratio) + proj_matrices.append(eigenvectors[:, ::-1].T) + + for _iter in range(self.max_iter): + for mode_i in range(1, n_dims): # ith mode + mode_cov_mat = np.zeros((shape_[mode_i], shape_[mode_i])) + proj_other = [proj_matrices[m] for m in range(n_dims - 1) if m != mode_i - 1] + modes_other = [m for m in range(1, n_dims) if m != mode_i] + for start in range(0, n_samples, chunk_size): + batch = X[start : start + chunk_size] - self.mean_ + batch_proj = multi_mode_dot(batch, proj_other, modes=modes_other) + batch_unfold = unfold(batch_proj, mode=mode_i) + mode_cov_mat += batch_unfold @ batch_unfold.T + + _, eigenvectors = eigh( + mode_cov_mat, + subset_by_index=[shape_[mode_i] - output_shape[mode_i - 1], shape_[mode_i] - 1], + ) + proj_matrices[mode_i - 1] = eigenvectors[:, ::-1].T + + # variance of the projected features, accumulated per chunk so the full + # unfolded projection never has to be materialized + x_proj_var = np.zeros(int(np.prod(output_shape))) + modes_all = [m for m in range(1, n_dims)] + for start in range(0, n_samples, chunk_size): + batch = X[start : start + chunk_size] - self.mean_ + batch_proj = multi_mode_dot(batch, proj_matrices, modes=modes_all) + batch_unfold = unfold(batch_proj, mode=0) # unfold the chunked projection to shape (n_chunk, n_features) + x_proj_var += np.einsum("ij,ij->j", batch_unfold, batch_unfold) + idx_order_ = np.argsort(-x_proj_var) + + self.proj_mats_ = proj_matrices + self.idx_order_ = idx_order_ + self.output_shape_ = output_shape + self.explained_variance_ratio_ = tuple(explained_variance_ratios) + self.n_dims_ = n_dims return self - def transform(self, X): + def transform(self, X, vectorize=None): """Project data to the MPCA subspace. Parameters ---------- X : ndarray of shape (n_samples, I_1, ..., I_N) or (I_1, ..., I_N) Input tensor data. + vectorize : bool, default=None + Whether to return the projected data as vectors. If ``None`` + (default), the value set in ``__init__`` (``self.vectorize``) is + used. Returns ------- - x_projected : ndarray + X_projected : ndarray Projected data. Shape is ``(n_samples, P_1, ..., P_N)`` when ``vectorize=False``. Otherwise returns vectorized features with optional truncation to ``n_components``. """ + if vectorize is None: + vectorize = self.vectorize # reshape X to shape (1, I_1, I_2, ..., I_N) if X in shape (I_1, I_2, ..., I_N), i.e. n_samples = 1 - if X.ndim == self.n_dims - 1: + if X.ndim == self.n_dims_ - 1: X = X.reshape((1,) + X.shape) - _check_tensor_dim_shape(X, self.n_dims, self.shape_in) + _check_tensor_dim_shape(X, self.n_dims_, self.input_shape_) X = X - self.mean_ # projected tensor in lower dimensions - x_projected = multi_mode_dot(X, self.proj_mats, modes=[m for m in range(1, self.n_dims)]) + X_projected = multi_mode_dot(X, self.proj_mats_, modes=[m for m in range(1, self.n_dims_)]) n_components = self.n_components - if self.vectorize: - x_projected = unfold(x_projected, mode=0) - x_projected = x_projected[:, self.idx_order] + if vectorize: + X_projected = unfold(X_projected, mode=0) + X_projected = X_projected[:, self.idx_order_] if isinstance(n_components, int): - n_features = int(np.prod(self.shape_out)) + n_features = int(np.prod(self.output_shape_)) if n_components > n_features: warn_msg = ( "n_components %d exceeds the maximum number, all features will be returned." % n_components @@ -265,9 +385,9 @@ def transform(self, X): logging.warning(warn_msg) warnings.warn(warn_msg) n_components = n_features - x_projected = x_projected[:, :n_components] + X_projected = X_projected[:, :n_components] - return x_projected + return X_projected def inverse_transform(self, X): """Reconstruct original-space tensors from projected data. @@ -279,28 +399,28 @@ def inverse_transform(self, X): Returns ------- - x_rec : ndarray of shape (n_samples, I_1, ..., I_N) + X_reconstructed : ndarray of shape (n_samples, I_1, ..., I_N) Reconstructed tensor data in the original shape. """ - # reshape X to tensor in shape (n_samples, self.shape_out) if X has been unfolded + # reshape X to tensor in shape (n_samples, self.output_shape_) if X has been unfolded if X.ndim <= 2: if X.ndim == 1: # reshape X to a 2D matrix (1, n_components) if X in shape (n_components,) X = X.reshape((1, -1)) n_samples = X.shape[0] n_features = X.shape[1] - if n_features <= np.prod(self.shape_out): - x_ = np.zeros((n_samples, np.prod(self.shape_out))) - x_[:, self.idx_order[:n_features]] = X[:] + if n_features <= np.prod(self.output_shape_): + x_ = np.zeros((n_samples, np.prod(self.output_shape_))) + x_[:, self.idx_order_[:n_features]] = X[:] else: msg = "Feature dimension exceeds the shape upper limit." logging.error(msg) raise ValueError(msg) - X = fold(x_, mode=0, shape=((n_samples,) + self.shape_out)) + X = fold(x_, mode=0, shape=((n_samples,) + self.output_shape_)) - x_rec = multi_mode_dot(X, self.proj_mats, modes=[m for m in range(1, self.n_dims)], transpose=True) + X_reconstructed = multi_mode_dot(X, self.proj_mats_, modes=[m for m in range(1, self.n_dims_)], transpose=True) - x_rec = x_rec + self.mean_ + X_reconstructed = X_reconstructed + self.mean_ - return x_rec + return X_reconstructed diff --git a/kalelinear/transformer/_tca.py b/kalelinear/transformer/_tca.py index 4f36413..5e946a6 100644 --- a/kalelinear/transformer/_tca.py +++ b/kalelinear/transformer/_tca.py @@ -73,7 +73,7 @@ def _make_eigenproblem(self, x_kernel_matrix, context): return obj, st y_kernel_matrix = self.gamma_ * np.dot(context.y_encoded, context.y_encoded.T) + (1 - self.gamma_) * identity - lap_mat = lap_norm(context.X_fit, n_neighbour=self.k, mode="connectivity") + lap_mat = lap_norm(context.X_fit, n_neighbors=self.k, mode="connectivity") obj += multi_dot([x_kernel_matrix, (mmd_matrix + self.mu * lap_mat), x_kernel_matrix]) st += multi_dot([x_kernel_matrix, h, y_kernel_matrix, h, x_kernel_matrix]) diff --git a/kalelinear/utils/_base.py b/kalelinear/utils/_base.py index 8069dab..a376195 100644 --- a/kalelinear/utils/_base.py +++ b/kalelinear/utils/_base.py @@ -8,14 +8,14 @@ from kalelinear.utils._backend import to_numpy -def lap_norm(X, n_neighbour=3, metric="cosine", mode="distance", normalise=True): +def lap_norm(X, n_neighbors=3, metric="cosine", mode="distance", normalize=True): """[summary] Parameters ---------- X : [type] [description] - n_neighbour : int, optional + n_neighbors : int, optional [description], by default 3 metric : str, optional [description], by default 'cosine' @@ -24,7 +24,7 @@ def lap_norm(X, n_neighbour=3, metric="cosine", mode="distance", normalise=True) returned matrix: 'connectivity' will return the connectivity matrix with ones and zeros, and 'distance' will return the distances between neighbors according to the given metric. - normalise : bool, optional + normalize : bool, optional [description], by default True Returns @@ -34,7 +34,7 @@ def lap_norm(X, n_neighbour=3, metric="cosine", mode="distance", normalise=True) """ x_np = to_numpy(X) n = x_np.shape[0] - knn_graph = kneighbors_graph(x_np, n_neighbour, metric=metric, mode=mode).toarray() + knn_graph = kneighbors_graph(x_np, n_neighbors, metric=metric, mode=mode).toarray() W = np.zeros((n, n)) knn_idx = np.logical_or(knn_graph, knn_graph.T) if mode == "distance": @@ -44,7 +44,7 @@ def lap_norm(X, n_neighbour=3, metric="cosine", mode="distance", normalise=True) W[knn_idx] = 1 D = np.diag(np.sum(W, axis=1)) - if normalise: + if normalize: D_ = inv(sqrtm(D)) lap_mat = np.eye(n) - multi_dot([D_, W, D_]) else: diff --git a/setup.py b/setup.py index 1c42d1b..fb32d8d 100644 --- a/setup.py +++ b/setup.py @@ -11,15 +11,7 @@ # Core dependencies frequently used in the kalelinear API -install_requires = [ - "cvxopt", - "numpy", - "osqp", - "pandas", - "scikit-learn>=1.6.0", - "scipy", - "tensorly", -] +install_requires = ["cvxopt", "numpy", "osqp", "pandas", "scikit-learn>=1.6.0", "scipy", "tensorly"] # Dependencies for all examples and tutorials example_requires = [ diff --git a/tests/conftest.py b/tests/conftest.py index 3a6dc8d..a03745f 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -2,6 +2,7 @@ import os import pytest +from scipy.io import loadmat @pytest.fixture(scope="session") @@ -9,3 +10,11 @@ def download_path(): path = os.path.join("tests", "test_data") os.makedirs(path, exist_ok=True) return path + + +@pytest.fixture(scope="module") +def gait(download_path): + gait_data_path = os.path.join(download_path, "gait.mat") + if not os.path.exists(gait_data_path): + pytest.skip("Gait data not found, skipping tests.") + return loadmat(gait_data_path) diff --git a/tests/estimator/test_estimator.py b/tests/estimator/test_estimator.py index cd994dd..342667a 100644 --- a/tests/estimator/test_estimator.py +++ b/tests/estimator/test_estimator.py @@ -7,7 +7,7 @@ @pytest.fixture def office_test_data(): - x = np.array( + X = np.array( [ [-2.2, -1.9], [-1.9, -2.1], @@ -24,62 +24,62 @@ def office_test_data(): y = np.array([0, 0, 0, 1, 1, 1, 0, 0, 1, 1]) z = np.array([1, 1, 1, 1, 1, 1, 0, 0, 0, 0]) covariate_mat = np.eye(2)[z] - return x, y, z, covariate_mat + return X, y, z, covariate_mat def _split_source_target(office_test_data, target_domain=0): - x, y, z, covariate_mat = office_test_data + X, y, z, covariate_mat = office_test_data tgt_idx = np.where(z == target_domain) src_idx = np.where(z != target_domain) - x_train = np.concatenate((x[src_idx], x[tgt_idx])) + X_train = np.concatenate((X[src_idx], X[tgt_idx])) c_train = np.concatenate((covariate_mat[src_idx], covariate_mat[tgt_idx])) y_train = y[src_idx] - return x, y, tgt_idx, src_idx, x_train, c_train, y_train + return X, y, tgt_idx, src_idx, X_train, c_train, y_train def test_coir_svm_solvers_fit_consistently(office_test_data): - x, y, tgt_idx, src_idx, x_train, c_train, y_train = _split_source_target(office_test_data) + X, y, tgt_idx, src_idx, X_train, c_train, y_train = _split_source_target(office_test_data) clf1 = estimator.CoIRSVM() clf2 = estimator.CoIRSVM(solver="cvxopt") - clf1.fit(x_train, y_train, c_train) - clf2.fit(x_train, y_train, c_train) + clf1.fit(X_train, y_train, c_train) + clf2.fit(X_train, y_train, c_train) for clf in (clf1, clf2): - y_pred = clf.predict(x[tgt_idx]) - assert clf.coef_.shape[0] == x_train.shape[0] + y_pred = clf.predict(X[tgt_idx]) + assert clf.coef_.shape[0] == X_train.shape[0] assert np.isfinite(clf.coef_).all() assert set(np.unique(y_pred)).issubset(set(np.unique(y))) def test_coir_svm_none_covariates_matches_zero_covariates(office_test_data): - x, _, tgt_idx, _, x_train, _, y_train = _split_source_target(office_test_data) - zero_covariates = np.zeros((x_train.shape[0], 1)) + X, _, tgt_idx, _, X_train, _, y_train = _split_source_target(office_test_data) + zero_covariates = np.zeros((X_train.shape[0], 1)) - clf_none = estimator.CoIRSVM().fit(x_train, y_train, covariates=None) - clf_zero = estimator.CoIRSVM().fit(x_train, y_train, covariates=zero_covariates) + clf_none = estimator.CoIRSVM().fit(X_train, y_train, covariates=None) + clf_zero = estimator.CoIRSVM().fit(X_train, y_train, covariates=zero_covariates) - dec_none = clf_none.decision_function(x[tgt_idx]) - dec_zero = clf_zero.decision_function(x[tgt_idx]) - pred_none = clf_none.predict(x[tgt_idx]) - pred_zero = clf_zero.predict(x[tgt_idx]) + dec_none = clf_none.decision_function(X[tgt_idx]) + dec_zero = clf_zero.decision_function(X[tgt_idx]) + pred_none = clf_none.predict(X[tgt_idx]) + pred_zero = clf_zero.predict(X[tgt_idx]) assert np.allclose(dec_none, dec_zero) assert np.array_equal(pred_none, pred_zero) def test_coir_ls_none_covariates_matches_zero_covariates(office_test_data): - x, _, tgt_idx, _, x_train, _, y_train = _split_source_target(office_test_data) - zero_covariates = np.zeros((x_train.shape[0], 1)) + X, _, tgt_idx, _, X_train, _, y_train = _split_source_target(office_test_data) + zero_covariates = np.zeros((X_train.shape[0], 1)) - clf_none = estimator.CoIRLS().fit(x_train, y_train, covariates=None) - clf_zero = estimator.CoIRLS().fit(x_train, y_train, covariates=zero_covariates) + clf_none = estimator.CoIRLS().fit(X_train, y_train, covariates=None) + clf_zero = estimator.CoIRLS().fit(X_train, y_train, covariates=zero_covariates) - dec_none = clf_none.decision_function(x[tgt_idx]) - dec_zero = clf_zero.decision_function(x[tgt_idx]) - pred_none = clf_none.predict(x[tgt_idx]) - pred_zero = clf_zero.predict(x[tgt_idx]) + dec_none = clf_none.decision_function(X[tgt_idx]) + dec_zero = clf_zero.decision_function(X[tgt_idx]) + pred_none = clf_none.predict(X[tgt_idx]) + pred_zero = clf_zero.predict(X[tgt_idx]) assert np.allclose(dec_none, dec_zero) assert np.array_equal(pred_none, pred_zero) @@ -87,13 +87,13 @@ def test_coir_ls_none_covariates_matches_zero_covariates(office_test_data): @pytest.mark.parametrize("estimator_cls", [estimator.CoIRSVM, estimator.CoIRLS]) def test_coir_estimators_predict_labels(estimator_cls, office_test_data): - x, y, tgt_idx, src_idx, x_train, c_train, y_train = _split_source_target(office_test_data) + X, y, tgt_idx, src_idx, X_train, c_train, y_train = _split_source_target(office_test_data) clf = estimator_cls() - clf.fit(x_train, y_train, c_train) + clf.fit(X_train, y_train, c_train) - decision = clf.decision_function(x[tgt_idx]) - y_pred = clf.predict(x[tgt_idx]) + decision = clf.decision_function(X[tgt_idx]) + y_pred = clf.predict(X[tgt_idx]) acc = accuracy_score(y[tgt_idx], y_pred) assert decision.shape[0] == len(tgt_idx[0]) @@ -103,33 +103,33 @@ def test_coir_estimators_predict_labels(estimator_cls, office_test_data): @pytest.mark.parametrize("estimator_cls", [estimator.CoIRSVM, estimator.CoIRLS]) def test_coir_covariate_encoder_accepts_string_covariates(estimator_cls, office_test_data): - x, y, z, _ = office_test_data + X, y, z, _ = office_test_data tgt_idx = np.where(z == 0) src_idx = np.where(z != 0) - x_train = np.concatenate((x[src_idx], x[tgt_idx])) + X_train = np.concatenate((X[src_idx], X[tgt_idx])) z_train = np.concatenate((z[src_idx], z[tgt_idx])) y_train = y[src_idx] string_covariates = np.where(z_train == 0, "target", "source") numeric_covariates = np.eye(2)[z_train] - string_clf = estimator_cls(covariate_encoder="onehot").fit(x_train, y_train, string_covariates) - numeric_clf = estimator_cls().fit(x_train, y_train, numeric_covariates) + string_clf = estimator_cls(covariate_encoder="onehot").fit(X_train, y_train, string_covariates) + numeric_clf = estimator_cls().fit(X_train, y_train, numeric_covariates) assert string_clf.covariate_encoder_ is not None - assert np.allclose(string_clf.decision_function(x[tgt_idx]), numeric_clf.decision_function(x[tgt_idx])) + assert np.allclose(string_clf.decision_function(X[tgt_idx]), numeric_clf.decision_function(X[tgt_idx])) @pytest.mark.parametrize("estimator_cls", [estimator.ARSVM, estimator.ARRLS]) def test_artl_estimators_predict_labels(estimator_cls, office_test_data): - x, y, z, _ = office_test_data + X, y, z, _ = office_test_data tgt_idx = np.where(z == 0) src_idx = np.where(z != 0) clf = estimator_cls() - clf.fit(x, y[src_idx], covariates=z, target_covariate=0) + clf.fit(X, y[src_idx], covariates=z, target_covariate=0) - decision = clf.decision_function(x[tgt_idx]) - y_pred = clf.predict(x[tgt_idx]) + decision = clf.decision_function(X[tgt_idx]) + y_pred = clf.predict(X[tgt_idx]) acc = accuracy_score(y[tgt_idx], y_pred) assert decision.shape[0] == len(tgt_idx[0]) @@ -139,36 +139,36 @@ def test_artl_estimators_predict_labels(estimator_cls, office_test_data): @pytest.mark.parametrize("estimator_cls", [estimator.ARSVM, estimator.ARRLS]) def test_artl_covariate_api_accepts_source_only_or_full_labels(estimator_cls, office_test_data): - x, y, z, _ = office_test_data + X, y, z, _ = office_test_data tgt_idx = np.where(z == 0) src_idx = np.where(z != 0) - covariate = estimator_cls().fit(x, y[src_idx], covariates=z, target_covariate=0) + covariate = estimator_cls().fit(X, y[src_idx], covariates=z, target_covariate=0) full_y = y.copy() full_y[tgt_idx] = -1 - covariate_full_y = estimator_cls().fit(x, full_y, covariates=z, target_covariate=0, unlabeled_value=-1) + covariate_full_y = estimator_cls().fit(X, full_y, covariates=z, target_covariate=0, unlabeled_value=-1) - assert np.allclose(covariate.decision_function(x[tgt_idx]), covariate_full_y.decision_function(x[tgt_idx])) - assert np.array_equal(covariate.predict(x[tgt_idx]), covariate_full_y.predict(x[tgt_idx])) + assert np.allclose(covariate.decision_function(X[tgt_idx]), covariate_full_y.decision_function(X[tgt_idx])) + assert np.array_equal(covariate.predict(X[tgt_idx]), covariate_full_y.predict(X[tgt_idx])) @pytest.mark.parametrize("estimator_cls", [estimator.ARSVM, estimator.ARRLS]) def test_artl_rejects_legacy_target_keyword(estimator_cls, office_test_data): - x, y, z, _ = office_test_data + X, y, z, _ = office_test_data tgt_idx = np.where(z == 0) src_idx = np.where(z != 0) with pytest.raises(TypeError, match="unexpected keyword argument 'Xt'"): - estimator_cls().fit(x[src_idx], y[src_idx], Xt=x[tgt_idx]) + estimator_cls().fit(X[src_idx], y[src_idx], Xt=X[tgt_idx]) @pytest.mark.parametrize("estimator_cls", [estimator.ARSVM, estimator.ARRLS]) def test_artl_covariate_fit_predict_returns_target_labels(estimator_cls, office_test_data): - x, y, z, _ = office_test_data + X, y, z, _ = office_test_data tgt_idx = np.where(z == 0) src_idx = np.where(z != 0) - y_pred = estimator_cls().fit_predict(x, y[src_idx], covariates=z, target_covariate=0) + y_pred = estimator_cls().fit_predict(X, y[src_idx], covariates=z, target_covariate=0) assert y_pred.shape == y[tgt_idx].shape assert set(np.unique(y_pred)).issubset(set(np.unique(y))) @@ -176,13 +176,13 @@ def test_artl_covariate_fit_predict_returns_target_labels(estimator_cls, office_ @pytest.mark.parametrize("estimator_cls", [estimator.LapSVM, estimator.LapRLS]) def test_manifold_estimators_predict_labels(estimator_cls, office_test_data): - x, y, tgt_idx, src_idx, x_train, _, y_train = _split_source_target(office_test_data) + X, y, tgt_idx, src_idx, X_train, _, y_train = _split_source_target(office_test_data) clf = estimator_cls() - clf.fit(x_train, y_train) + clf.fit(X_train, y_train) - decision = clf.decision_function(x[tgt_idx]) - y_pred = clf.predict(x[tgt_idx]) + decision = clf.decision_function(X[tgt_idx]) + y_pred = clf.predict(X[tgt_idx]) acc = accuracy_score(y[tgt_idx], y_pred) assert decision.shape[0] == len(tgt_idx[0]) @@ -191,18 +191,18 @@ def test_manifold_estimators_predict_labels(estimator_cls, office_test_data): def test_gsda_fit_predicts_target_labels(office_test_data): - x, y, z, covariate_mat = office_test_data + X, y, z, covariate_mat = office_test_data target_idx = np.where(z == 0)[0] clf = estimator.GSDA(max_iter=25, random_state=0) - clf.fit(x, y[target_idx], covariate_mat, target_idx=target_idx) + clf.fit(X, y[target_idx], covariate_mat, target_idx=target_idx) - y_proba = clf.predict_proba(x[target_idx]) - y_pred = clf.predict(x[target_idx]) + y_proba = clf.predict_proba(X[target_idx]) + y_pred = clf.predict(X[target_idx]) params = clf.get_fitted_params() - assert clf.coef_.shape == (x.shape[1],) - assert params["coef"].shape == (x.shape[1],) + assert clf.coef_.shape == (X.shape[1],) + assert params["coef"].shape == (X.shape[1],) assert np.isfinite(clf.coef_).all() assert np.isfinite(clf.intercept_) assert y_proba.shape == y[target_idx].shape @@ -214,14 +214,14 @@ def test_gsda_fit_predicts_target_labels(office_test_data): def test_gsda_covariate_encoder_accepts_string_groups(office_test_data): - x, y, z, _ = office_test_data + X, y, z, _ = office_test_data target_idx = np.where(z == 0)[0] string_groups = np.where(z == 0, "target", "source") clf = estimator.GSDA(max_iter=25, random_state=0, covariate_encoder="onehot") - clf.fit(x, y[target_idx], string_groups, target_idx=target_idx) + clf.fit(X, y[target_idx], string_groups, target_idx=target_idx) - y_pred = clf.predict(x[target_idx]) + y_pred = clf.predict(X[target_idx]) assert clf.covariate_encoder_ is not None assert y_pred.shape == y[target_idx].shape diff --git a/tests/pipeline/test_mpca_trainer.py b/tests/pipeline/test_mpca_trainer.py new file mode 100644 index 0000000..1433663 --- /dev/null +++ b/tests/pipeline/test_mpca_trainer.py @@ -0,0 +1,60 @@ +import numpy as np +import pytest +from numpy import testing +from sklearn.metrics import accuracy_score, roc_auc_score + +from kalelinear.pipeline.mpca_trainer import MPCATrainer + +CLASSIFIERS = ["svc", "linear_svc", "lr"] +PARAMS = [ + {"classifier_params": "auto", "mpca_params": None, "n_features": None, "search_params": None}, + { + "classifier_params": {"C": 1}, + "mpca_params": {"explained_variance_ratio": 0.9, "vectorize": True}, + "n_features": 100, + "search_params": {"cv": 3}, + }, +] + + +@pytest.mark.parametrize("classifier", CLASSIFIERS) +@pytest.mark.parametrize("params", PARAMS) +def test_mpca_trainer(classifier, params, gait): + X = gait["fea3D"].transpose((3, 0, 1, 2)) + X = X[:20, :] + y = gait["gnd"][:20].reshape(-1) + trainer = MPCATrainer(classifier=classifier, **params) + trainer.fit(X, y) + y_pred = trainer.predict(X) + testing.assert_equal(np.unique(y), np.unique(y_pred)) + assert accuracy_score(y, y_pred) >= 0.8 + + if classifier == "linear_svc": + with pytest.raises(Exception): + y_proba = trainer.predict_proba(X) + else: + y_proba = trainer.predict_proba(X) + assert np.max(y_proba) <= 1.0 + assert np.min(y_proba) >= 0.0 + y_ = np.zeros(y.shape) + y_[np.where(y == 1)] = 1 + assert roc_auc_score(y_, y_proba[:, 0]) >= 0.8 + + y_dec_score = trainer.decision_function(X) + assert roc_auc_score(y, y_dec_score) >= 0.8 + + if classifier == "svc" and trainer.clf.kernel == "rbf": + with pytest.raises(Exception): + trainer.mpca.inverse_transform(trainer.clf.coef_) + else: + # interpret utilities (select_top_weight/plot_weights) are not ported + # to kalelinear yet, so only check the inverse-transform path here. + weights = trainer.mpca.inverse_transform(trainer.clf.coef_) - trainer.mpca.mean_ + testing.assert_equal(weights.shape[1:], X.shape[1:]) + + +def test_invalid_init(): + with pytest.raises(Exception): + MPCATrainer(classifier="Ridge") + with pytest.raises(Exception): + MPCATrainer(classifier_params=False) diff --git a/tests/transformer/test_mida.py b/tests/transformer/test_mida.py index 276c429..4da2cd3 100644 --- a/tests/transformer/test_mida.py +++ b/tests/transformer/test_mida.py @@ -13,7 +13,7 @@ def sample_data(): # Test an extreme case of domain shift # yet the data's manifold is linearly separable - x, y, domains = make_domain_shifted_dataset( + X, y, domains = make_domain_shifted_dataset( num_domains=10, num_samples_per_class=2, num_features=20, @@ -23,63 +23,63 @@ def sample_data(): covariates = OneHotEncoder(handle_unknown="ignore").fit_transform(domains.reshape(-1, 1)).toarray() - return x, y, domains, covariates + return X, y, domains, covariates @pytest.mark.parametrize("num_components", [2, None]) def test_mida_shape_consistency(sample_data, num_components): - x, y, domains, covariates = sample_data + X, y, domains, covariates = sample_data mida = MIDA(n_components=num_components) mida.set_params(**mida.get_params()) - mida.fit(x, covariates=covariates) + mida.fit(X, covariates=covariates) # Transform the whole data - z = mida.transform(x, covariates=covariates) + z = mida.transform(X, covariates=covariates) # If num_components is not None, check the shape of the transformed data if num_components is not None: - testing.assert_equal(z.shape, (len(x), num_components)) + testing.assert_equal(z.shape, (len(X), num_components)) # Transform the source and target domain data separately source_mask = domains != 0 # return True/false mask instead of indices - z_src = mida.transform(x[source_mask], covariates=covariates[source_mask]) - z_tgt = mida.transform(x[~source_mask], covariates=covariates[~source_mask]) + z_src = mida.transform(X[source_mask], covariates=covariates[source_mask]) + z_tgt = mida.transform(X[~source_mask], covariates=covariates[~source_mask]) # Check if transformations are consistent with separate domains testing.assert_allclose(z_src, z[source_mask]) testing.assert_allclose(z_tgt, z[~source_mask]) orig_coef_dim = mida.orig_coef_.shape[1] - feature_dim = x.shape[1] + feature_dim = X.shape[1] assert mida.orig_coef_ is not None, "MIDA must have `orig_coef_` after fitting when kernel='linear'" assert orig_coef_dim == feature_dim, f"orig_coef_ shape mismatch: {orig_coef_dim} != {feature_dim}" def test_mida_inverse_transform(sample_data): - x, y, domains, covariates = sample_data + X, y, domains, covariates = sample_data mida = MIDA(fit_inverse_transform=True) - mida.fit(x, covariates=covariates) + mida.fit(X, covariates=covariates) # Transform the whole data - z = mida.transform(x, covariates=covariates) + z = mida.transform(X, covariates=covariates) # Inverse transform the data - x_rec = mida.inverse_transform(z) + X_rec = mida.inverse_transform(z) # We don't check whether the inverse transform is exactly equal to the original data # in terms of value since it is expected to be different due to the domain adaptation effect. # We only check the shape and dimensionality. - assert len(x_rec) == len(x), f"Inverse transform failed: {len(x_rec)} != {len(x)}" - assert x_rec.ndim == x.ndim, f"Inverse transform failed: {x_rec.ndim} != {x.ndim}" - testing.assert_equal(x_rec.shape, x.shape) + assert len(X_rec) == len(X), f"Inverse transform failed: {len(X_rec)} != {len(X)}" + assert X_rec.ndim == X.ndim, f"Inverse transform failed: {X_rec.ndim} != {X.ndim}" + testing.assert_equal(X_rec.shape, X.shape) @pytest.mark.parametrize("kernel", ["linear", "rbf", "cosine"]) def test_mida_support_kernel(sample_data, kernel): - x, y, domains, covariates = sample_data + X, y, domains, covariates = sample_data mida = MIDA(n_components=N_COMP_CONSTANT, kernel=kernel) - mida.fit(x, covariates=covariates) + mida.fit(X, covariates=covariates) is_linear = kernel == "linear" try: @@ -89,7 +89,7 @@ def test_mida_support_kernel(sample_data, kernel): assert not is_linear, "MIDA must not have `orig_coef_` after fitting when kernel!='linear'" # Transform the whole data - z = mida.transform(x, covariates=covariates) + z = mida.transform(X, covariates=covariates) # expect to allow multiple kernels supported assert mida._n_features_out == N_COMP_CONSTANT, f"Expected {N_COMP_CONSTANT} components, got {mida._n_features_out}" @@ -98,10 +98,10 @@ def test_mida_support_kernel(sample_data, kernel): @pytest.mark.parametrize("augment", ["pre", "post", None]) def test_mida_augment(sample_data, augment): - x, y, domains, covariates = sample_data + X, y, domains, covariates = sample_data mida = MIDA(n_components=N_COMP_CONSTANT, kernel="linear", augment=augment, fit_inverse_transform=True) - mida.fit(x, covariates=covariates) + mida.fit(X, covariates=covariates) # expect validator for domain covariates if augment=True has_validator = hasattr(mida, "_factor_validator") @@ -111,18 +111,18 @@ def test_mida_augment(sample_data, augment): assert not has_validator, "MIDA must not have `_factor_validator` after fitting when augment=False" if augment == "post": - z = mida.transform(x, covariates=covariates) + z = mida.transform(X, covariates=covariates) assert ( z.shape[1] == N_COMP_CONSTANT + covariates.shape[-1] - ), f"Expected {x.shape[1]} features, got {N_COMP_CONSTANT + covariates.shape[-1]}" + ), f"Expected {X.shape[1]} features, got {N_COMP_CONSTANT + covariates.shape[-1]}" @pytest.mark.parametrize("ignore_y", [True, False]) def test_mida_ignore_y(sample_data, ignore_y): - x, y, domains, covariates = sample_data + X, y, domains, covariates = sample_data mida = MIDA(n_components=N_COMP_CONSTANT, kernel="linear", ignore_y=ignore_y) - mida.fit(x, y, covariates=covariates) + mida.fit(X, y, covariates=covariates) # expect classes_ to be set if ignore_y=False has_classes = hasattr(mida, "classes_") @@ -134,48 +134,48 @@ def test_mida_ignore_y(sample_data, ignore_y): def test_mida_covariate_encoder_onehot(sample_data): - x, y, domains, _ = sample_data + X, y, domains, _ = sample_data mida = MIDA(n_components=2, covariate_encoder="onehot") - z = mida.fit_transform(x, y=y, covariates=domains) + z = mida.fit_transform(X, y=y, covariates=domains) - assert z.shape == (len(x), 2) + assert z.shape == (len(X), 2) assert hasattr(mida, "covariate_encoder_") assert mida.covariate_encoder_ is not None - testing.assert_equal(mida.covariates_fit_.shape[0], len(x)) + testing.assert_equal(mida.covariates_fit_.shape[0], len(X)) def test_mida_covariate_encoder_onehot_accepts_strings(sample_data): - x, y, domains, _ = sample_data + X, y, domains, _ = sample_data string_domains = np.asarray([f"domain-{domain}" for domain in domains]) mida = MIDA(n_components=2, covariate_encoder="onehot") - z = mida.fit_transform(x, y=y, covariates=string_domains) + z = mida.fit_transform(X, y=y, covariates=string_domains) - assert z.shape == (len(x), 2) + assert z.shape == (len(X), 2) assert np.issubdtype(mida.covariates_fit_.dtype, np.number) def test_mida_requires_numeric_covariates_without_encoder(sample_data): - x, y, domains, _ = sample_data + X, y, domains, _ = sample_data with pytest.raises(ValueError, match="covariate_encoder"): - MIDA(n_components=2).fit(x, y=y, covariates=domains.astype(str)) + MIDA(n_components=2).fit(X, y=y, covariates=domains.astype(str)) def test_mida_transform_rejects_mismatched_covariate_dimension(sample_data): - x, _, _, covariates = sample_data + X, _, _, covariates = sample_data - mida = MIDA(n_components=2, augment="pre").fit(x, covariates=covariates) + mida = MIDA(n_components=2, augment="pre").fit(X, covariates=covariates) - mismatched_covariates = np.hstack((covariates, np.zeros((len(x), 1)))) + mismatched_covariates = np.hstack((covariates, np.zeros((len(X), 1)))) with pytest.raises(ValueError, match="feature dimension"): - mida.transform(x, covariates=mismatched_covariates) + mida.transform(X, covariates=mismatched_covariates) @pytest.mark.parametrize("eigen_solver", ["auto", "dense", "arpack", "randomized"]) def test_mida_eigen_solver(sample_data, eigen_solver): - x, y, domains, covariates = sample_data + X, y, domains, covariates = sample_data mida = MIDA( n_components=N_COMP_CONSTANT, @@ -183,27 +183,27 @@ def test_mida_eigen_solver(sample_data, eigen_solver): eigen_solver=eigen_solver, max_iter=200 if eigen_solver == "arpack" else None, ) - mida.fit(x, covariates=covariates) + mida.fit(X, covariates=covariates) # Transform the whole data - z = mida.transform(x, covariates=covariates) + z = mida.transform(X, covariates=covariates) # expect the solver to have a consistent number of components assert mida._n_features_out == N_COMP_CONSTANT, f"Expected {N_COMP_CONSTANT} components, got {mida._n_features_out}" - assert z.shape[1] == mida._n_features_out, f"Expected {x.shape[1]} features, got {mida._n_features_out}" + assert z.shape[1] == mida._n_features_out, f"Expected {X.shape[1]} features, got {mida._n_features_out}" @pytest.mark.parametrize("scale_components", [True, False]) def test_mida_scale_components(sample_data, scale_components): - x, y, domains, covariates = sample_data + X, y, domains, covariates = sample_data mida = MIDA(n_components=N_COMP_CONSTANT, kernel="linear", scale_components=scale_components) - mida.fit(x, covariates=covariates) + mida.fit(X, covariates=covariates) # Transform the whole data - z = mida.transform(x, covariates=covariates) + z = mida.transform(X, covariates=covariates) # Expect the scale_components to have a consistent number of components # the behavior expected is the zero eigenvalues component is masked, not indexed assert mida._n_features_out == N_COMP_CONSTANT, f"Expected {N_COMP_CONSTANT} components, got {mida._n_features_out}" - assert z.shape[1] == mida._n_features_out, f"Expected {x.shape[1]} features, got {mida._n_features_out}" + assert z.shape[1] == mida._n_features_out, f"Expected {X.shape[1]} features, got {mida._n_features_out}" diff --git a/tests/transformer/test_mpca.py b/tests/transformer/test_mpca.py index 4b1ce1e..2a324df 100644 --- a/tests/transformer/test_mpca.py +++ b/tests/transformer/test_mpca.py @@ -10,15 +10,7 @@ N_COMPS = [1, 50, 100] VAR_RATIOS = [0.7, 0.95] -relative_tol = 0.00001 - - -@pytest.fixture(scope="module") -def gait(download_path): - gait_data_path = os.path.join(download_path, "gait.mat") - if not os.path.exists(gait_data_path): - pytest.skip("Gait data not found, skipping tests.") - return loadmat(gait_data_path) +RELATIVE_TOL = 0.00001 @pytest.fixture(scope="module") @@ -30,59 +22,90 @@ def baseline_model(download_path): @pytest.mark.parametrize("n_components", N_COMPS) -@pytest.mark.parametrize("var_ratio", VAR_RATIOS) -def test_mpca(var_ratio, n_components, gait): +@pytest.mark.parametrize("explained_variance_ratio", VAR_RATIOS) +def test_mpca(explained_variance_ratio, n_components, gait): # basic mpca test, return tensor - x = gait["fea3D"].transpose((3, 0, 1, 2)) - mpca = MPCA(var_ratio=var_ratio, vectorize=False) - x_proj = mpca.fit(x).transform(x) + X = gait["fea3D"].transpose((3, 0, 1, 2)) + mpca = MPCA(explained_variance_ratio=explained_variance_ratio, vectorize=False) + X_proj = mpca.fit(X).transform(X) - testing.assert_equal(x_proj.ndim, x.ndim) - testing.assert_equal(x_proj.shape[0], x.shape[0]) - for i in range(1, x.ndim): - assert x_proj.shape[i] <= x.shape[i] - testing.assert_equal(mpca.proj_mats[i - 1].shape[1], x.shape[i]) + testing.assert_equal(X_proj.ndim, X.ndim) + testing.assert_equal(X_proj.shape[0], X.shape[0]) + for i in range(1, X.ndim): + assert X_proj.shape[i] <= X.shape[i] + testing.assert_equal(mpca.proj_mats_[i - 1].shape[1], X.shape[i]) - x_rec = mpca.inverse_transform(x_proj) - testing.assert_equal(x_rec.shape, x.shape) + X_rec = mpca.inverse_transform(X_proj) + testing.assert_equal(X_rec.shape, X.shape) # test return vector mpca.set_params(**{"vectorize": True, "n_components": n_components}) - x_proj = mpca.transform(x) - testing.assert_equal(x_proj.ndim, 2) - testing.assert_equal(x_proj.shape[0], x.shape[0]) - testing.assert_equal(x_proj.shape[1], n_components) - x_rec = mpca.inverse_transform(x_proj) - testing.assert_equal(x_rec.shape, x.shape) + X_proj = mpca.transform(X) + testing.assert_equal(X_proj.ndim, 2) + testing.assert_equal(X_proj.shape[0], X.shape[0]) + testing.assert_equal(X_proj.shape[1], n_components) + X_rec = mpca.inverse_transform(X_proj) + testing.assert_equal(X_rec.shape, X.shape) # test n_samples = 1 - x0_proj = mpca.transform(x[0]) - testing.assert_equal(x0_proj.ndim, 2) - testing.assert_equal(x0_proj.shape[0], 1) - testing.assert_equal(x0_proj.shape[1], n_components) - x0_rec = mpca.inverse_transform(x0_proj.reshape(-1)) - testing.assert_equal(x0_rec.shape[1:], x[0].shape) + X0_proj = mpca.transform(X[0]) + testing.assert_equal(X0_proj.ndim, 2) + testing.assert_equal(X0_proj.shape[0], 1) + testing.assert_equal(X0_proj.shape[1], n_components) + X0_rec = mpca.inverse_transform(X0_proj.reshape(-1)) + testing.assert_equal(X0_rec.shape[1:], X[0].shape) # test n_components exceeds upper limit - mpca.set_params(**{"vectorize": True, "n_components": np.prod(x.shape[1:]) + 1}) - x_proj = mpca.transform(x) - testing.assert_equal(x_proj.shape[1], np.prod(mpca.shape_out)) + mpca.set_params(**{"vectorize": True, "n_components": np.prod(X.shape[1:]) + 1}) + X_proj = mpca.transform(X) + testing.assert_equal(X_proj.shape[1], np.prod(mpca.output_shape_)) def test_mpca_against_baseline(gait, baseline_model): - x = gait["fea3D"].transpose((3, 0, 1, 2)) + X = gait["fea3D"].transpose((3, 0, 1, 2)) baseline_proj_mats = [baseline_model["tUs"][i][0] for i in range(baseline_model["tUs"].size)] baseline_mean = baseline_model["TXmean"] - mpca = MPCA(var_ratio=0.97) - x_proj = mpca.fit(x).transform(x) + mpca = MPCA(explained_variance_ratio=0.97) + X_proj = mpca.fit(X).transform(X) testing.assert_allclose(baseline_mean, mpca.mean_) - baseline_proj_x = multi_mode_dot(x - baseline_mean, baseline_proj_mats, modes=[1, 2, 3]) + baseline_proj_X = multi_mode_dot(X - baseline_mean, baseline_proj_mats, modes=[1, 2, 3]) # check whether the output embeddings is close to the baseline output by keeping the same variance ratio 97% - testing.assert_allclose(x_proj**2, baseline_proj_x**2, rtol=relative_tol) - # testing.assert_equal(x_proj.shape, baseline_proj_x.shape) + testing.assert_allclose(X_proj**2, baseline_proj_X**2, rtol=RELATIVE_TOL) + # testing.assert_equal(X_proj.shape, baseline_proj_X.shape) - for i in range(x.ndim - 1): + for i in range(X.ndim - 1): # check whether each eigen-vector column is equal to/opposite of corresponding baseline eigen-vector column - # testing.assert_allclose(abs(mpca.proj_mats[i]), abs(baseline_proj_mats[i])) - testing.assert_allclose(mpca.proj_mats[i] ** 2, baseline_proj_mats[i] ** 2, rtol=relative_tol) + # testing.assert_allclose(abs(mpca.proj_mats_[i]), abs(baseline_proj_mats[i])) + testing.assert_allclose(mpca.proj_mats_[i] ** 2, baseline_proj_mats[i] ** 2, rtol=RELATIVE_TOL) + + +def test_transform_vectorize_override(gait): + X = gait["fea3D"].transpose((3, 0, 1, 2)) + n_components = 50 + + # init-level default and transform-level override to True + mpca = MPCA(vectorize=False).fit(X) + X_proj_tensor = mpca.transform(X) + X_proj_vec = mpca.transform(X, vectorize=True) + testing.assert_equal(X_proj_tensor.ndim, X.ndim) + testing.assert_equal(X_proj_vec.ndim, 2) + + # init-level vectorize=True and transform-level override to False + mpca = MPCA(vectorize=True, n_components=n_components).fit(X) + X_proj_vec = mpca.transform(X) + X_proj_tensor = mpca.transform(X, vectorize=False) + testing.assert_equal(X_proj_vec.ndim, 2) + testing.assert_equal(X_proj_vec.shape[1], n_components) + testing.assert_equal(X_proj_tensor.ndim, X.ndim) + + # explicit None keeps the init-level setting + X_proj_tensor_default = mpca.transform(X, vectorize=None) + testing.assert_equal(X_proj_tensor_default.ndim, 2) + + +def test_fit_empty_input_raises(): + X = np.empty((0, 4, 5, 6)) + mpca = MPCA() + with pytest.raises(ValueError, match="MPCA requires at least 2 samples to fit."): + mpca.fit(X) diff --git a/tests/transformer/test_transformer.py b/tests/transformer/test_transformer.py index 73accb2..4cb33be 100644 --- a/tests/transformer/test_transformer.py +++ b/tests/transformer/test_transformer.py @@ -7,7 +7,7 @@ @pytest.fixture def domain_adaptation_data(): - xs = np.array( + Xs = np.array( [ [-2.0, -1.8], [-1.8, -2.1], @@ -16,7 +16,7 @@ def domain_adaptation_data(): ] ) ys = np.array([0, 0, 1, 1]) - xt = np.array( + Xt = np.array( [ [-1.4, -1.2], [-1.2, -1.1], @@ -25,7 +25,7 @@ def domain_adaptation_data(): ] ) yt = np.array([0, 0, 1, 1]) - x = np.vstack((xs, xt)) + X = np.vstack((Xs, Xt)) y = np.concatenate((ys, yt)) binary_covariates = np.array([0, 0, 0, 0, 1, 1, 1, 1]) covariates = np.array( @@ -40,54 +40,54 @@ def domain_adaptation_data(): [0.0, 1.0], ] ) - return xs, ys, xt, yt, x, y, binary_covariates, covariates + return Xs, ys, Xt, yt, X, y, binary_covariates, covariates @pytest.mark.parametrize("transformer_cls", [TCA, JDA, BDA]) def test_domain_transformers_fit_transform_shapes(transformer_cls, domain_adaptation_data): - xs, ys, xt, yt, x, y, binary_covariates, _ = domain_adaptation_data + Xs, ys, Xt, yt, X, y, binary_covariates, _ = domain_adaptation_data transformer = transformer_cls(n_components=2) - x_transformed = transformer.fit_transform(x, y=y, covariates=binary_covariates, target_covariate=1) - xs_transformed = x_transformed[binary_covariates == 0] - xt_transformed = x_transformed[binary_covariates == 1] + X_transformed = transformer.fit_transform(X, y=y, covariates=binary_covariates, target_covariate=1) + Xs_transformed = X_transformed[binary_covariates == 0] + Xt_transformed = X_transformed[binary_covariates == 1] - assert xs_transformed.shape == (xs.shape[0], 2) - assert xt_transformed.shape == (xt.shape[0], 2) - assert np.isfinite(xs_transformed).all() - assert np.isfinite(xt_transformed).all() + assert Xs_transformed.shape == (Xs.shape[0], 2) + assert Xt_transformed.shape == (Xt.shape[0], 2) + assert np.isfinite(Xs_transformed).all() + assert np.isfinite(Xt_transformed).all() def test_mida_fit_transform_shapes_with_covariates(domain_adaptation_data): - _, _, _, _, x, y, _, covariates = domain_adaptation_data + _, _, _, _, X, y, _, covariates = domain_adaptation_data transformer = MIDA(n_components=2) - x_transformed = transformer.fit_transform(x, y=y, covariates=covariates) + X_transformed = transformer.fit_transform(X, y=y, covariates=covariates) - assert x_transformed.shape == (x.shape[0], 2) - assert np.isfinite(x_transformed).all() + assert X_transformed.shape == (X.shape[0], 2) + assert np.isfinite(X_transformed).all() def test_mida_transform_requires_fit(domain_adaptation_data): - xs, _, _, _, _, _, _, _ = domain_adaptation_data + Xs, _, _, _, _, _, _, _ = domain_adaptation_data transformer = MIDA(n_components=2) with pytest.raises(NotFittedError): - transformer.transform(xs) + transformer.transform(Xs) def test_transformers_store_training_projection(domain_adaptation_data): - xs, _, xt, _, x, y, binary_covariates, covariates = domain_adaptation_data + Xs, _, Xt, _, X, y, binary_covariates, covariates = domain_adaptation_data - tca = TCA(n_components=2).fit(x, y=y, covariates=binary_covariates, target_covariate=1) - jda = JDA(n_components=2).fit(x, y=y, covariates=binary_covariates, target_covariate=1) - bda = BDA(n_components=2, mu=0.25).fit(x, y=y, covariates=binary_covariates, target_covariate=1) - mida = MIDA(n_components=2).fit(x, y=y, covariates=covariates) + tca = TCA(n_components=2).fit(X, y=y, covariates=binary_covariates, target_covariate=1) + jda = JDA(n_components=2).fit(X, y=y, covariates=binary_covariates, target_covariate=1) + bda = BDA(n_components=2, mu=0.25).fit(X, y=y, covariates=binary_covariates, target_covariate=1) + mida = MIDA(n_components=2).fit(X, y=y, covariates=covariates) - assert tca.U.shape[0] == xs.shape[0] + xt.shape[0] - assert jda.U.shape[0] == xs.shape[0] + xt.shape[0] - assert bda.U.shape[0] == xs.shape[0] + xt.shape[0] - assert mida.U.shape[0] == xs.shape[0] + xt.shape[0] + assert tca.U.shape[0] == Xs.shape[0] + Xt.shape[0] + assert jda.U.shape[0] == Xs.shape[0] + Xt.shape[0] + assert bda.U.shape[0] == Xs.shape[0] + Xt.shape[0] + assert mida.U.shape[0] == Xs.shape[0] + Xt.shape[0] @pytest.mark.parametrize("transformer_cls", [TCA, JDA, BDA]) @@ -100,18 +100,18 @@ def test_mmd_transformers_reject_covariate_encoder(transformer_cls, domain_adapt @pytest.mark.parametrize("transformer_cls", [TCA, JDA, BDA]) def test_mmd_transformers_require_both_domains(transformer_cls, domain_adaptation_data): - _, _, _, _, x, _, _, _ = domain_adaptation_data + _, _, _, _, X, _, _, _ = domain_adaptation_data with pytest.raises(ValueError, match="both source and target"): - transformer_cls(n_components=2).fit(x, covariates=np.zeros(x.shape[0], dtype=int)) + transformer_cls(n_components=2).fit(X, covariates=np.zeros(X.shape[0], dtype=int)) @pytest.mark.parametrize("transformer_cls", [TCA, JDA, BDA]) def test_mmd_transformers_validate_target_covariate(transformer_cls, domain_adaptation_data): - _, _, _, _, x, _, binary_covariates, _ = domain_adaptation_data + _, _, _, _, X, _, binary_covariates, _ = domain_adaptation_data with pytest.raises(ValueError, match="target_covariate"): - transformer_cls(n_components=2).fit(x, covariates=binary_covariates, target_covariate=2) + transformer_cls(n_components=2).fit(X, covariates=binary_covariates, target_covariate=2) def test_jda_does_not_accept_mu(): @@ -121,7 +121,7 @@ def test_jda_does_not_accept_mu(): @pytest.mark.parametrize("mu", [-0.1, 1.1]) def test_bda_validates_mu(mu, domain_adaptation_data): - _, _, _, _, x, y, binary_covariates, _ = domain_adaptation_data + _, _, _, _, X, y, binary_covariates, _ = domain_adaptation_data with pytest.raises(ValueError, match="mu"): - BDA(n_components=2, mu=mu).fit(x, y=y, covariates=binary_covariates, target_covariate=1) + BDA(n_components=2, mu=mu).fit(X, y=y, covariates=binary_covariates, target_covariate=1) diff --git a/tests/utils/test_utils.py b/tests/utils/test_utils.py index c78318f..c1776a7 100644 --- a/tests/utils/test_utils.py +++ b/tests/utils/test_utils.py @@ -29,7 +29,7 @@ def test_base_init_returns_expected_shapes(sample_data): @pytest.mark.parametrize("mode", ["distance", "connectivity"]) def test_lap_norm_returns_square_matrix(sample_data, mode): - lap = lap_norm(sample_data, n_neighbour=2, mode=mode, normalise=False) + lap = lap_norm(sample_data, n_neighbors=2, mode=mode, normalize=False) assert lap.shape == (4, 4) assert np.allclose(lap, lap.T) @@ -80,7 +80,7 @@ def make_domain_shifted_dataset( w = random_state.randn(num_features) w = w / np.linalg.norm(w) - x_all = [] + X_all = [] y_all = [] domain_all = [] @@ -90,21 +90,21 @@ def make_domain_shifted_dataset( for label in [0, 1]: class_mean = (label - 0.5) * class_sep * w + domain_shift cov = np.eye(num_features) - x_class = random_state.multivariate_normal(class_mean, cov, num_samples_per_class) + X_class = random_state.multivariate_normal(class_mean, cov, num_samples_per_class) y_class = np.full(num_samples_per_class, label) domain_class = np.full(num_samples_per_class, i_domain) - x_all.append(x_class) + X_all.append(X_class) y_all.append(y_class) domain_all.append(domain_class) - x = np.vstack(x_all) + X = np.vstack(X_all) y = np.concatenate(y_all) domains = np.concatenate(domain_all) - idx = random_state.permutation(len(x)) - x = x[idx] + idx = random_state.permutation(len(X)) + X = X[idx] y = y[idx] domains = domains[idx] - return x, y, domains + return X, y, domains