Repository: Steenroder/steenroder (language: python)

Issue:
Create python-publish.yml



Discussion / hints:

Fix patch:
diff --git a/.gitattributes b/.gitattributes
new file mode 100644
index 0000000..409014f
--- /dev/null
+++ b/.gitattributes
@@ -0,0 +1,2 @@
+docs/_notebooks/ran/flat_Klein_bottle.ipynb filter=lfs diff=lfs merge=lfs -text
+docs/_notebooks/ran/cyclo-octane.ipynb filter=lfs diff=lfs merge=lfs -text
diff --git a/.github/ISSUE_TEMPLATE/bug_report.md b/.github/ISSUE_TEMPLATE/bug_report.md
new file mode 100644
index 0000000..dd84ea7
--- /dev/null
+++ b/.github/ISSUE_TEMPLATE/bug_report.md
@@ -0,0 +1,38 @@
+---
+name: Bug report
+about: Create a report to help us improve
+title: ''
+labels: ''
+assignees: ''
+
+---
+
+**Describe the bug**
+A clear and concise description of what the bug is.
+
+**To Reproduce**
+Steps to reproduce the behavior:
+1. Go to '...'
+2. Click on '....'
+3. Scroll down to '....'
+4. See error
+
+**Expected behavior**
+A clear and concise description of what you expected to happen.
+
+**Screenshots**
+If applicable, add screenshots to help explain your problem.
+
+**Desktop (please complete the following information):**
+ - OS: [e.g. iOS]
+ - Browser [e.g. chrome, safari]
+ - Version [e.g. 22]
+
+**Smartphone (please complete the following information):**
+ - Device: [e.g. iPhone6]
+ - OS: [e.g. iOS8.1]
+ - Browser [e.g. stock browser, safari]
+ - Version [e.g. 22]
+
+**Additional context**
+Add any other context about the problem here.
diff --git a/.github/ISSUE_TEMPLATE/feature_request.md b/.github/ISSUE_TEMPLATE/feature_request.md
new file mode 100644
index 0000000..bbcbbe7
--- /dev/null
+++ b/.github/ISSUE_TEMPLATE/feature_request.md
@@ -0,0 +1,20 @@
+---
+name: Feature request
+about: Suggest an idea for this project
+title: ''
+labels: ''
+assignees: ''
+
+---
+
+**Is your feature request related to a problem? Please describe.**
+A clear and concise description of what the problem is. Ex. I'm always frustrated when [...]
+
+**Describe the solution you'd like**
+A clear and concise description of what you want to happen.
+
+**Describe alternatives you've considered**
+A clear and concise description of any alternative solutions or features you've considered.
+
+**Additional context**
+Add any other context or screenshots about the feature request here.
diff --git a/.github/PULL_REQUEST_TEMPLATE.md b/.github/PULL_REQUEST_TEMPLATE.md
new file mode 100644
index 0000000..d056c21
--- /dev/null
+++ b/.github/PULL_REQUEST_TEMPLATE.md
@@ -0,0 +1,54 @@
+<!--
+Thanks for contributing a pull request! Please ensure you have taken a look at
+the guidelines for contributing: https://giotto-ai.github.io/gtda-docs/latest/contributing/#guidelines
+-->
+
+**Reference issues/PRs**
+<!--
+Example: Fixes #1234. See also #3456.
+Please use keywords (e.g., Fixes) to create link to the issues or pull requests
+you resolved, so that they will automatically be closed when your pull request
+is merged. See https://github.com/blog/1506-closing-issues-via-pull-requests.
+
+Pull requests which are not related to open issues may be rejected!
+If suggesting a new feature or change, please discuss it in an issue first.
+If fixing a bug, there should be an issue describing it with steps to reproduce.
+-->
+
+**Types of changes**
+<!--
+What types of changes does your code introduce? Put an `x` in all the boxes that apply.
+-->
+- [ ] Bug fix (non-breaking change which fixes an issue)
+- [ ] New feature (non-breaking change which adds functionality)
+- [ ] Breaking change (fix or feature that would cause existing functionality to change)
+
+**Description**
+<!--
+Describe your changes in detail.
+-->
+
+**Screenshots (if appropriate)**
+
+**Any other comments?**
+
+**Checklist**
+<!--
+Go over all the following points, and put an `x` in all the boxes that apply. 
+If you're unsure about any of these, don't hesitate to ask. We're here to help!
+-->
+- [ ] I have read the [guidelines for contributing](https://github.com/Steenroder/steenroder/CONTRIBUTING.md).
+- [ ] My code follows the code style of this project. I used `flake8` to check my Python changes.
+- [ ] My change requires a change to the documentation.
+- [ ] I have updated the documentation accordingly.
+- [ ] I have added tests to cover my changes.
+- [ ] All new and existing tests passed. I used `pytest` to check this on Python tests.
+
+<!--
+We value all user contributions, no matter how minor they are. If we are slow to
+review, either the pull request needs some benchmarking, tinkering,
+convincing, etc. or more likely the reviewers are simply busy. In either
+case, we ask for your understanding during the review process.
+
+Thanks for contributing!
+-->
diff --git a/.github/workflows/build-sphinx.yml b/.github/workflows/build-sphinx.yml
new file mode 100644
index 0000000..697490d
--- /dev/null
+++ b/.github/workflows/build-sphinx.yml
@@ -0,0 +1,46 @@
+name: build-sphinx-to-gh-pages
+
+env:
+  GITHUB_ACTOR: Steenroder
+  GITHUB_REPOSITORY: Steenroder/steenroder
+  GITHUB_TOKEN: ${{ secrets.CREATE_DEMO_SECRET }}
+
+on:
+  push:
+    branches: [ main ]
+
+jobs:
+  docs:
+    runs-on: ubuntu-latest
+    steps:
+    - uses: actions/checkout@v2
+      with:
+        lfs: true
+    - name: Checkout LFS files
+      run: git lfs checkout
+    - name: Set up Python
+      uses: actions/setup-python@v2
+      with:
+        python-version: "3.8"
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        python -m pip install .[tests]
+    - name: Test notebooks with pytest
+      run: |
+        notebooks=$(ls notebooks/*.ipynb)
+        long_notebooks=("notebooks/flat_Klein_bottle.ipynb" "notebooks/cyclo-octane.ipynb")
+        tested_notebooks=(`echo ${notebooks[@]} ${long_notebooks[@]} | tr ' ' '\n' | sort | uniq -u `)
+        pytest --nbmake --overwrite -n=auto ${tested_notebooks[@]}
+        cp ${tested_notebooks[@]} docs/_notebooks/
+    - name: Build documentation
+      uses: sphinx-toolbox/sphinx-action@master
+      with:
+        pre-build-command: "python -m pip install .[docs]"
+        docs-folder: "docs/"
+    - name: Deploy
+      uses: peaceiris/actions-gh-pages@v3
+      if: github.ref == 'refs/heads/main'
+      with:
+        github_token: ${{ secrets.GITHUB_TOKEN }}
+        publish_dir: docs/_build/html
diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
new file mode 100644
index 0000000..39d7e7a
--- /dev/null
+++ b/.github/workflows/python-package.yml
@@ -0,0 +1,40 @@
+# This workflow will install Python dependencies, run tests and lint with a variety of Python versions
+# For more information see: https://help.github.com/actions/language-and-framework-guides/using-python-with-github-actions
+
+name: Python package
+
+on:
+  push:
+    branches: [ main ]
+  pull_request:
+    branches: [ main ]
+
+jobs:
+  build:
+
+    runs-on: ubuntu-latest
+    strategy:
+      fail-fast: false
+      matrix:
+        python-version: ["3.8", "3.9", "3.10"]
+
+    steps:
+    - uses: actions/checkout@v2
+    - name: Set up Python ${{ matrix.python-version }}
+      uses: actions/setup-python@v2
+      with:
+        python-version: ${{ matrix.python-version }}
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        python -m pip install flake8 pytest
+        if [ -f requirements.txt ]; then pip install -r requirements.txt; fi
+    - name: Lint with flake8
+      run: |
+        # stop the build if there are Python syntax errors or undefined names
+        flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
+        # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
+        flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
+    - name: Test with pytest
+      run: |
+        pytest
diff --git a/.github/workflows/python-publish.yml b/.github/workflows/python-publish.yml
new file mode 100644
index 0000000..3bfabfc
--- /dev/null
+++ b/.github/workflows/python-publish.yml
@@ -0,0 +1,36 @@
+# This workflow will upload a Python Package using Twine when a release is created
+# For more information see: https://help.github.com/en/actions/language-and-framework-guides/using-python-with-github-actions#publishing-to-package-registries
+
+# This workflow uses actions that are not certified by GitHub.
+# They are provided by a third-party and are governed by
+# separate terms of service, privacy policy, and support
+# documentation.
+
+name: Upload Python Package
+
+on:
+  release:
+    types: [published]
+
+jobs:
+  deploy:
+
+    runs-on: ubuntu-latest
+
+    steps:
+    - uses: actions/checkout@v2
+    - name: Set up Python
+      uses: actions/setup-python@v2
+      with:
+        python-version: '3.x'
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install build
+    - name: Build package
+      run: python -m build
+    - name: Publish package
+      uses: pypa/gh-action-pypi-publish@27b31702a0e7fc50959f5ad993c78deac1bdfc29
+      with:
+        user: __token__
+        password: ${{ secrets.PYPI_API_TOKEN }}
diff --git a/.gitignore b/.gitignore
index 117da6c..0f38dad 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,114 @@
-*.pyc
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+.hypothesis/
+.pytest_cache/
+test-output.xml
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+docs/generated/
+docs/*.html
+docs/*.txt
+
+# PyBuilder
+target/
+
+# Jupyter Notebook
 .ipynb_checkpoints
-*DS_Store*
-__pycache__
-.idea
-steenroder.egg-info
+
+# pyenv
+.python-version
+
+# celery beat schedule file
+celerybeat-schedule
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+
+# Emacs
+*~
+
+# JetBrains
+.idea/
diff --git a/CODE_OF_CONDUCT.md b/CODE_OF_CONDUCT.md
new file mode 100644
index 0000000..fecaf4b
--- /dev/null
+++ b/CODE_OF_CONDUCT.md
@@ -0,0 +1,128 @@
+# Contributor Covenant Code of Conduct
+
+## Our Pledge
+
+We as members, contributors, and leaders pledge to make participation in our
+community a harassment-free experience for everyone, regardless of age, body
+size, visible or invisible disability, ethnicity, sex characteristics, gender
+identity and expression, level of experience, education, socio-economic status,
+nationality, personal appearance, race, religion, or sexual identity
+and orientation.
+
+We pledge to act and interact in ways that contribute to an open, welcoming,
+diverse, inclusive, and healthy community.
+
+## Our Standards
+
+Examples of behavior that contributes to a positive environment for our
+community include:
+
+* Demonstrating empathy and kindness toward other people
+* Being respectful of differing opinions, viewpoints, and experiences
+* Giving and gracefully accepting constructive feedback
+* Accepting responsibility and apologizing to those affected by our mistakes,
+  and learning from the experience
+* Focusing on what is best not just for us as individuals, but for the
+  overall community
+
+Examples of unacceptable behavior include:
+
+* The use of sexualized language or imagery, and sexual attention or
+  advances of any kind
+* Trolling, insulting or derogatory comments, and personal or political attacks
+* Public or private harassment
+* Publishing others' private information, such as a physical or email
+  address, without their explicit permission
+* Other conduct which could reasonably be considered inappropriate in a
+  professional setting
+
+## Enforcement Responsibilities
+
+Community leaders are responsible for clarifying and enforcing our standards of
+acceptable behavior and will take appropriate and fair corrective action in
+response to any behavior that they deem inappropriate, threatening, offensive,
+or harmful.
+
+Community leaders have the right and responsibility to remove, edit, or reject
+comments, commits, code, wiki edits, issues, and other contributions that are
+not aligned to this Code of Conduct, and will communicate reasons for moderation
+decisions when appropriate.
+
+## Scope
+
+This Code of Conduct applies within all community spaces, and also applies when
+an individual is officially representing the community in public spaces.
+Examples of representing our community include using an official e-mail address,
+posting via an official social media account, or acting as an appointed
+representative at an online or offline event.
+
+## Enforcement
+
+Instances of abusive, harassing, or otherwise unacceptable behavior may be
+reported to the community leaders responsible for enforcement at
+steenroder@gmail.com.
+All complaints will be reviewed and investigated promptly and fairly.
+
+All community leaders are obligated to respect the privacy and security of the
+reporter of any incident.
+
+## Enforcement Guidelines
+
+Community leaders will follow these Community Impact Guidelines in determining
+the consequences for any action they deem in violation of this Code of Conduct:
+
+### 1. Correction
+
+**Community Impact**: Use of inappropriate language or other behavior deemed
+unprofessional or unwelcome in the community.
+
+**Consequence**: A private, written warning from community leaders, providing
+clarity around the nature of the violation and an explanation of why the
+behavior was inappropriate. A public apology may be requested.
+
+### 2. Warning
+
+**Community Impact**: A violation through a single incident or series
+of actions.
+
+**Consequence**: A warning with consequences for continued behavior. No
+interaction with the people involved, including unsolicited interaction with
+those enforcing the Code of Conduct, for a specified period of time. This
+includes avoiding interactions in community spaces as well as external channels
+like social media. Violating these terms may lead to a temporary or
+permanent ban.
+
+### 3. Temporary Ban
+
+**Community Impact**: A serious violation of community standards, including
+sustained inappropriate behavior.
+
+**Consequence**: A temporary ban from any sort of interaction or public
+communication with the community for a specified period of time. No public or
+private interaction with the people involved, including unsolicited interaction
+with those enforcing the Code of Conduct, is allowed during this period.
+Violating these terms may lead to a permanent ban.
+
+### 4. Permanent Ban
+
+**Community Impact**: Demonstrating a pattern of violation of community
+standards, including sustained inappropriate behavior,  harassment of an
+individual, or aggression toward or disparagement of classes of individuals.
+
+**Consequence**: A permanent ban from any sort of public interaction within
+the community.
+
+## Attribution
+
+This Code of Conduct is adapted from the [Contributor Covenant][homepage],
+version 2.0, available at
+https://www.contributor-covenant.org/version/2/0/code_of_conduct.html.
+
+Community Impact Guidelines were inspired by [Mozilla's code of conduct
+enforcement ladder](https://github.com/mozilla/diversity).
+
+[homepage]: https://www.contributor-covenant.org
+
+For answers to common questions about this code of conduct, see the FAQ at
+https://www.contributor-covenant.org/faq. Translations are available at
+https://www.contributor-covenant.org/translations.
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
new file mode 100644
index 0000000..ae60d46
--- /dev/null
+++ b/CONTRIBUTING.md
@@ -0,0 +1,40 @@
+# Contributing guidelines
+
+We welcome any kind of contribution to our steenroder, from simple comment or question to a full fledged [pull request](https://help.github.com/articles/about-pull-requests/). Please read and follow our [Code of Conduct](CODE_OF_CONDUCT.md).
+
+A contribution can be one of the following cases:
+
+1. you have a question;
+1. you think you may have found a bug (including unexpected behavior);
+1. you want to make some kind of change to the code base (e.g. to fix a bug, to add a new feature, to update documentation).
+
+The sections below outline the steps in each case.
+
+## You have a question
+
+1. use the search functionality [here](https://github.com/Steenroder/steenroder/issues) to see if someone already filed the same issue;
+1. if your issue search did not yield any relevant results, make a new issue;
+1. apply the "Question" label; apply other labels when relevant.
+
+## You think you may have found a bug
+
+1. use the search functionality [here](https://github.com/Steenroder/steenroder/issues) to see if someone already filed the same issue;
+1. if your issue search did not yield any relevant results, make a new issue, making sure to provide enough information to the rest of the community to understand the cause and context of the problem. Depending on the issue, you may want to include:
+    - the [SHA hashcode](https://help.github.com/articles/autolinked-references-and-urls/#commit-shas) of the commit that is causing your problem;
+    - some identifying information (name and version number) for dependencies you're using;
+    - information about the operating system;
+1. apply relevant labels to the newly created issue.
+
+## You want to make some kind of change to the code base
+
+1. (**important**) announce your plan to the rest of the community _before you start working_. This announcement should be in the form of a (new) issue;
+1. (**important**) wait until some kind of consensus is reached about your idea being a good idea;
+1. if needed, fork the repository to your own Github profile and create your own feature branch off of the latest main commit. While working on your feature branch, make sure to stay up to date with the main branch by pulling in changes, possibly from the 'upstream' repository (follow the instructions [here](https://help.github.com/articles/configuring-a-remote-for-a-fork/) and [here](https://help.github.com/articles/syncing-a-fork/));
+1. Install dependencies with `python -m pip install -r requirements.txt`;
+1. make sure the existing tests still work by running ``pytest``. If project tests fails use ``pytest --keep-baked-projects`` to keep generated project in /tmp/pytest-* and investigate;
+1. add your own tests (if necessary);
+1. update or expand the documentation;
+1. push your feature branch to (your fork of) the Python Template repository on GitHub;
+1. create the pull request, e.g. following the instructions [here](https://help.github.com/articles/creating-a-pull-request/).
+
+In case you feel like you've made a valuable contribution, but you don't know how to write or run tests for it, or how to generate the documentation: don't let this discourage you from making the pull request; we can help you! Just go ahead and submit the pull request, but keep in mind that you might be asked to append additional commits to your pull request.
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..a8035dc
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2022 Steenroder
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/README.md b/README.md
deleted file mode 100644
index 010f29a..0000000
--- a/README.md
+++ /dev/null
@@ -1,5 +0,0 @@
-# steenroder
-
-The widespread use in applied topology of the barcode of filtered cellular complexes rests on a balance between discriminatory power and computability. It has long been envision that the strength of this invariant could be increase using cohomology operations. This package computes the recently defined *Sq*<sup>*k*</sup>-barcodes which have been shown to effectively increase the discriminatory power of barcodes on real-world data.
-For a complete presentation of these invariants please consult [Persistence Steenrod modules](https://arxiv.org/abs/1812.05031) by U. Lupo, A. Medina-Mardones and G. Tauzin.
-
diff --git a/README.rst b/README.rst
new file mode 100644
index 0000000..a4e191f
--- /dev/null
+++ b/README.rst
@@ -0,0 +1,109 @@
+steenroder: A Python package for the computation of persistence Steenrod barcodes
+=================================================================================
+
+The widespread use in applied topology of the barcode of filtered
+cellular complexes rests on a balance between discriminatory power and
+computability. It has long been envision that the strength of this
+invariant could be increase using cohomology operations. This package
+computes the recently defined *Sq*\ \ *k*\ -barcodes which have been
+shown to effectively increase the discriminatory power of barcodes on
+real-world data.
+
+For a complete presentation of these invariants please consult
+`Persistence Steenrod modules <https://arxiv.org/abs/1812.05031>`__ by
+U. Lupo, A. Medina-Mardones and G. Tauzin.
+
+License
+-------
+
+``steenroder`` is distributed under the `MIT
+license <https://github.com/Steenroder/steenroder/LICENSE>`__.
+
+Documentation
+-------------
+
+Please visit https://steenroder.github.io/steenroder/ and navigate to
+the version you are interested in.
+
+Installation
+------------
+
+Dependencies
+~~~~~~~~~~~~
+
+The latest stable version of ``steenroder`` requires:
+
+-  python (>= 3.8)
+-  numpy (>= 1.19.1)
+-  numba (>= 0.53.0)
+-  psutils (>= 5.8.0)
+-  gudhi (>= 3.5.0)
+-  plotly (>= 5.3.1)
+
+To run the examples, jupyter is required.
+
+.. _installation-1:
+
+Installation
+~~~~~~~~~~~~
+
+The simplest way to install ``steenroder`` is using ``pip`` ::
+
+   python -m pip install -U steenroder
+
+If necessary, this will also automatically install all the above
+dependencies. Note: we recommend upgrading ``pip`` to a recent version
+as the above may fail on very old versions.
+
+Contributing
+------------
+
+We welcome new contributors of all experience levels. The Steenroder
+community goals are to be helpful, welcoming, and effective. To learn
+more about making a contribution to ``steenroder``, please consult the
+`relevant
+page <https://github.com/Steenroder/steenroder/CONTRIBUTING.md>`__.
+
+Testing
+~~~~~~~
+
+After developer installation, you can launch the test suite from outside
+the source directory:
+
+::
+
+   pytest steenroder
+
+Important links
+---------------
+
+-  Official source code repo: https://github.com/Steenroder/steenroder
+-  Download releases: https://pypi.org/project/steenroder/
+-  Issue tracker: https://github.com/Steenroder/steenroder/issues
+
+Citing steenroder
+-----------------
+
+If you use ``steenroder`` in a scientific publication, we would
+appreciate citations to the following paper:
+
+`Persistence Steenrod modules <https://arxiv.org/abs/1812.05031>`__
+
+You can use the following BibTeX entry:
+
+::
+
+   @article{steenroder,
+          author = {{Lupo}, Umberto and {Medina-Mardones}, Anibal M. and {Tauzin}, Guillaume},
+           title = "{Persistence Steenrod modules}",
+         journal = {arXiv e-prints},
+   archivePrefix = {arXiv},
+          eprint = {1812.05031},
+    primaryClass = {math.AT},
+          adsurl = {https://ui.adsabs.harvard.edu/abs/2018arXiv181205031L},
+   }
+
+Contacts
+--------
+
+steenroder@gmail.com
diff --git a/docs/.githooks/clear_notebooks.py b/docs/.githooks/clear_notebooks.py
new file mode 100644
index 0000000..24c60ab
--- /dev/null
+++ b/docs/.githooks/clear_notebooks.py
@@ -0,0 +1,66 @@
+import argparse
+import nbformat
+
+from nbconvert.preprocessors import (
+    ClearOutputPreprocessor,
+    ExecutePreprocessor,
+)
+from nbconvert.preprocessors import CellExecutionError
+
+
+def process_notebook(
+        notebook_filename,
+        html_directory='notebook-html',
+        execute=True,
+):
+    """Checks if an IPython notebook runs without error from start to finish.
+    If so, writes the notebook to HTML (with outputs) and overwrites the
+    .ipynb file (without outputs).
+    """
+    with open(notebook_filename) as f:
+        nb = nbformat.read(f, as_version=4)
+
+    clear = ClearOutputPreprocessor()
+
+    try:
+        if execute:
+            # Check that the notebook runs
+            ep = ExecutePreprocessor(timeout=600, kernel_name='python3')
+            ep.preprocess(nb, {'metadata': {'path': ''}})
+
+        msg = ''
+
+    except CellExecutionError:
+        msg = f'\n  Error executing the notebook "{notebook_filename}".\n'
+        msg += f'  See notebook "{notebook_filename}" for the traceback.'
+
+    # Clear notebook outputs and save as .ipynb
+    clear.preprocess(nb, {})
+    with open(notebook_filename, mode='w', encoding='utf-8') as f:
+        nbformat.write(nb, f)
+
+    print(f"Processed {notebook_filename}{msg}")
+
+    return
+
+
+if __name__ == '__main__':
+
+    parser = argparse.ArgumentParser(description='read some notebok files')
+    parser.add_argument('notebooks',
+                        metavar='Notebooks',
+                        type=str,
+                        nargs='+',
+                        help='notebooks')
+
+    args = parser.parse_args()
+
+    notebooks = args.notebooks
+
+    for fn in notebooks:
+        if not fn.endswith('.ipynb'):
+            print(f'Error: file {fn} is not an IPython notebook.')
+            raise
+
+    for fn in notebooks:
+        process_notebook(fn, execute=False)
diff --git a/docs/.githooks/pre-commit b/docs/.githooks/pre-commit
new file mode 100644
index 0000000..e0480df
--- /dev/null
+++ b/docs/.githooks/pre-commit
@@ -0,0 +1,13 @@
+#!/bin/sh
+
+# Get all notebooks currently staged for commit
+notebooks=`git diff --cached --name-only --diff-filter=d | grep  .ipynb`
+
+# Checks if $notebooks is non-empty
+if [ ! -z "$notebooks" ]
+then
+    # Run the script to clear the notebooks
+    python .githooks/clear_notebooks.py $notebooks
+    # Re-add notebooks to staging area
+    git add $notebooks
+fi
diff --git a/docs/.nojekyll b/docs/.nojekyll
new file mode 100644
index 0000000..e69de29
diff --git a/docs/Makefile b/docs/Makefile
new file mode 100644
index 0000000..0cdb184
--- /dev/null
+++ b/docs/Makefile
@@ -0,0 +1,23 @@
+# Minimal makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line.
+SPHINXOPTS    =
+SPHINXBUILD   = sphinx-build
+SOURCEDIR     = .
+BUILDDIR      = _build
+
+# Put it first so that "make" without argument is like "make help".
+help:
+	@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
+
+.PHONY: help Makefile
+
+github:
+	@make html
+	@cp -a "$(BUILDDIR)"/html/. ./docs
+
+# Catch-all target: route all unknown targets to Sphinx using the new
+# "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).
+%: Makefile
+	@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
diff --git a/docs/_notebooks/ran/cyclo-octane.ipynb b/docs/_notebooks/ran/cyclo-octane.ipynb
new file mode 100644
index 0000000..da01460
--- /dev/null
+++ b/docs/_notebooks/ran/cyclo-octane.ipynb
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:53566a136f5593b4db51a40fd0049c120edb297ef7f54573cbd6bef1e9d35685
+size 6868515
diff --git a/docs/_notebooks/ran/flat_Klein_bottle.ipynb b/docs/_notebooks/ran/flat_Klein_bottle.ipynb
new file mode 100644
index 0000000..45dd457
--- /dev/null
+++ b/docs/_notebooks/ran/flat_Klein_bottle.ipynb
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:581560fc22a7f36963279982f56de2227aec5da9c84ccdad814ca47508afd3ef
+size 3976816
diff --git a/docs/_static/style.css b/docs/_static/style.css
new file mode 100644
index 0000000..3e44f4a
--- /dev/null
+++ b/docs/_static/style.css
@@ -0,0 +1,4 @@
+{% extends "!layout.html" %}
+{% block extrahead %}
+    <link href="{{ pathto("_static/style.css", True) }}" rel="stylesheet" type="text/css">
+{% endblock %}
diff --git a/docs/_templates/function.rst b/docs/_templates/function.rst
new file mode 100644
index 0000000..7131c79
--- /dev/null
+++ b/docs/_templates/function.rst
@@ -0,0 +1,14 @@
+:mod:`{{module}}`.{{objname}}
+{{ underline }}====================
+
+.. currentmodule:: {{ module }}
+
+.. autofunction:: {{ objname }}
+
+..
+    Exclude sphinx-gallery generated examples since we use binder for now
+    include:: {{module}}.{{objname}}.examples
+
+.. raw:: html
+
+    <div class="clearer"></div>
diff --git a/docs/api_reference.rst b/docs/api_reference.rst
new file mode 100644
index 0000000..af9157d
--- /dev/null
+++ b/docs/api_reference.rst
@@ -0,0 +1,53 @@
+API Reference
+=============
+
+.. toctree::
+   :maxdepth: 2
+
+.. automodule:: steenroder
+   :no-members:
+   :no-inherited-members:
+
+
+Core
+----
+.. currentmodule:: steenroder
+
+.. autosummary::
+   :toctree: generated/
+   :template: function.rst
+
+   compute_reduced_triangular
+   compute_barcode_and_coho_reps
+   compute_steenrod_matrix
+   compute_steenrod_barcode
+   barcodes
+   rips_barcodes
+
+
+Preprocessing
+-------------
+.. currentmodule:: steenroder
+
+.. autosummary::
+   :toctree: generated/
+   :template: function.rst
+
+    make_simplices
+    make_filtration
+    sort_filtration_by_dim
+    cone
+    suspension
+    sphere
+    wedge
+
+
+Visualization
+-------------
+.. currentmodule:: steenroder
+
+.. autosummary::
+   :toctree: generated/
+   :template: function.rst
+
+   plot_diagrams
diff --git a/docs/build_site.sh b/docs/build_site.sh
new file mode 100644
index 0000000..5616f2a
--- /dev/null
+++ b/docs/build_site.sh
@@ -0,0 +1,63 @@
+#!/bin/bash
+set -x
+
+apt-get update
+apt-get -y install git rsync python3-sphinx
+
+pwd ls -lah
+export SOURCE_DATE_EPOCH=$(git log -1 --pretty=%ct)
+
+##############
+# BUILD DOCS #
+##############
+
+# Python Sphinx, configured with source/conf.py
+# See https://www.sphinx-doc.org/
+make clean
+make html
+
+#######################
+# Update GitHub Pages #
+#######################
+
+git config --global user.name "${GITHUB_ACTOR}"
+git config --global user.email "${GITHUB_ACTOR}@users.noreply.github.com"
+
+docroot=`mktemp -d`
+rsync -av "build/html/" "${docroot}/"
+
+pushd "${docroot}"
+
+git init
+git remote add deploy "https://token:${GITHUB_TOKEN}@github.com/${GITHUB_REPOSITORY}.git"
+git checkout -b gh-pages
+
+# Adds .nojekyll file to the root to signal to GitHub that
+# directories that start with an underscore (_) can remain
+touch .nojekyll
+
+# Add README
+cat > README.md <<EOF
+# README for the GitHub Pages Branch
+This branch is simply a cache for the website served from https://Steenroder.github.io/steenroder/,
+and is  not intended to be viewed on github.com.
+
+For more information on how this site is built using Sphinx, Read the Docs, and GitHub Actions/Pages, see:
+ * https://www.docslikecode.com/articles/github-pages-python-sphinx/
+ * https://tech.michaelaltfield.net/2020/07/18/sphinx-rtd-github-pages-1
+EOF
+
+# Copy the resulting html pages built from Sphinx to the gh-pages branch
+git add .
+
+# Make a commit with changes and any new files
+msg="Updating Docs for commit ${GITHUB_SHA} made on `date -d"@${SOURCE_DATE_EPOCH}" --iso-8601=seconds` from ${GITHUB_REF} by ${GITHUB_ACTOR}"
+git commit -am "${msg}"
+
+# overwrite the contents of the gh-pages branch on our github.com repo
+git push deploy gh-pages --force
+
+popd # return to main repo sandbox root
+
+# exit cleanly
+exit 0
diff --git a/docs/conf.py b/docs/conf.py
new file mode 100644
index 0000000..63ab9e0
--- /dev/null
+++ b/docs/conf.py
@@ -0,0 +1,127 @@
+# -*- coding: utf-8 -*-
+#
+# Configuration file for the Sphinx documentation builder.
+#
+# This file does only contain a selection of the most common options. For a
+# full list see the documentation:
+# http://www.sphinx-doc.org/en/master/config
+
+# -- Path setup --------------------------------------------------------------
+
+# If extensions (or modules to document with autodoc) are in another directory,
+# add these directories to sys.path here. If the directory is relative to the
+# documentation root, use os.path.abspath to make it absolute, like shown here.
+#
+from steenroder import __version__
+
+import os
+import sys
+sys.path.insert(0, os.path.abspath('.'))
+
+
+# -- Project information -----------------------------------------------------
+
+project = 'steenroder'
+copyright = '2022, Steenroder'
+author = 'Umberto Lupo, Anibal Medina-Mardones, and Guillaume Tauzin'
+
+# The short X.Y version
+version = ''
+# The full version, including alpha/beta/rc tags
+release = __version__
+
+# -- General configuration ---------------------------------------------------
+
+# If your documentation needs a minimal Sphinx version, state it here.
+#
+# needs_sphinx = '1.0'
+
+# Add any Sphinx extension module names here, as strings. They can be
+# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
+# ones.
+extensions = [
+    'sphinx.ext.githubpages',
+    'sphinx.ext.autosummary',
+    'sphinx.ext.viewcode',
+    'sphinx.ext.mathjax',
+    'sphinx.ext.napoleon',
+]
+
+# For maths, use mathjax by default and svg if NO_MATHJAX env variable is set
+# (useful for viewing the doc offline)
+if os.environ.get('NO_MATHJAX'):
+    extensions.append('sphinx.ext.imgmath')
+    imgmath_image_format = 'svg'
+else:
+    extensions.append('sphinx.ext.mathjax')
+    mathjax_path = ('https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.0/'
+                    'MathJax.js?config=TeX-AMS_SVG')
+
+autodoc_default_options = {'members': True, 'inherited-members': True}
+
+# Add any paths that contain templates here, relative to this directory.
+templates_path = ['_templates']
+
+# generate autosummary even if no references
+autosummary_generate = True
+
+# The suffix of source filenames.
+source_suffix = '.rst'
+
+# The encoding of source files.
+# source_encoding = 'utf-8'
+
+# The master toctree document.
+master_doc = 'index'
+
+# List of patterns, relative to source directory, that match files and
+# directories to ignore when looking for source files.
+# This pattern also affects html_static_path and html_extra_path.
+exclude_patterns = ['_build', '**.ipynb_checkpoints']
+
+# If true, '()' will be appended to :func: etc. cross-reference text.
+add_function_parentheses = False
+
+# If true, the current module name will be prepended to all description
+# unit titles (such as .. function::).
+# add_module_names = True
+
+# If true, sectionauthor and moduleauthor directives will be shown in the
+# output. They are ignored by default.
+# show_authors = False
+
+# The name of the Pygments (syntax highlighting) style to use.
+pygments_style = 'sphinx'
+
+# A list of ignored prefixes for module index sorting.
+# modindex_common_prefix = []
+# -- Options for HTML output -------------------------------------------------
+
+# The theme to use for HTML and HTML Help pages.  See the documentation for
+# a list of builtin themes.
+#
+html_title = project + ' version ' + release
+html_theme = 'sphinx_rtd_theme'
+html_theme_options = {
+    'collapse_navigation': False,
+}
+# Add any paths that contain custom static files (such as style sheets) here,
+# relative to this directory. They are copied after the builtin static files,
+# so a file named "default.css" will overwrite the builtin "default.css".
+html_static_path = ['_static']
+
+# Custom sidebar templates, must be a dictionary that maps document names
+# to template names.
+#
+# The default sidebars (for documents that don't match any pattern) are
+# defined by theme itself.  Builtin themes are using these templates by
+# default: ``['localtoc.html', 'relations.html', 'sourcelink.html',
+# 'searchbox.html']``.
+#
+# html_sidebars = {}
+
+
+# -- Options for HTMLHelp output ---------------------------------------------
+
+# Output file base name for HTML help builder.
+htmlhelp_basename = 'steenroder'
diff --git a/docs/index.rst b/docs/index.rst
new file mode 100644
index 0000000..900c885
--- /dev/null
+++ b/docs/index.rst
@@ -0,0 +1,18 @@
+Welcome to steenroder's documentation!
+======================================
+
+.. toctree::
+   :maxdepth: 2
+   :caption: Contents:
+
+   introduction
+   user_guide
+   api_reference
+
+
+Indices and tables
+==================
+
+* :ref:`genindex`
+* :ref:`modindex`
+* :ref:`search`
diff --git a/docs/introduction.rst b/docs/introduction.rst
new file mode 100644
index 0000000..eac9809
--- /dev/null
+++ b/docs/introduction.rst
@@ -0,0 +1,97 @@
+Introduction
+============
+
+Context
+-------
+
+The widespread use in applied topology of the barcode of filtered
+cellular complexes rests on a balance between discriminatory power and
+computability. It has long been envision that the strength of this
+invariant could be increase using cohomology operations. This package
+computes the recently defined *Sq*\ \ *k*\ -barcodes which have been
+shown to effectively increase the discriminatory power of barcodes on
+real-world data.
+
+For a complete presentation of these invariants please consult
+`Persistence Steenrod modules <https://arxiv.org/abs/1812.05031>`__ by
+U. Lupo, A. Medina-Mardones and G. Tauzin.
+
+
+License
+-------
+
+``steenroder`` is distributed under the `MIT
+license <https://github.com/Steenroder/steenroder/LICENSE>`__.
+
+
+Installation
+------------
+
+Dependencies
+~~~~~~~~~~~~
+
+The latest stable version of ``steenroder`` requires:
+
+-  python (>= 3.8)
+-  numpy (>= 1.19.1)
+-  numba (>= 0.53.0)
+-  psutils (>= 5.8.0)
+-  gudhi (>= 3.5.0)
+-  plotly (>= 5.3.1)
+
+To run the examples, jupyter is required.
+
+.. _installation-1:
+
+Installation
+~~~~~~~~~~~~
+
+The simplest way to install ``steenroder`` is using ``pip`` ::
+
+   python -m pip install -U steenroder
+
+If necessary, this will also automatically install all the above
+dependencies. Note: we recommend upgrading ``pip`` to a recent version
+as the above may fail on very old versions.
+
+Contributing
+------------
+
+We welcome new contributors of all experience levels. The Steenroder
+community goals are to be helpful, welcoming, and effective. To learn
+more about making a contribution to ``steenroder``, please consult the
+`relevant
+page <https://github.com/Steenroder/steenroder/CONTRIBUTING.md>`__.
+
+Testing
+~~~~~~~
+
+After developer installation, you can launch the test suite from outside
+the source directory:
+
+::
+
+   pytest steenroder
+
+
+Citing steenroder
+-----------------
+
+If you use ``steenroder`` in a scientific publication, we would
+appreciate citations to the following paper:
+
+`Persistence Steenrod modules <https://arxiv.org/abs/1812.05031>`__
+
+You can use the following BibTeX entry:
+
+::
+
+   @article{steenroder,
+          author = {{Lupo}, Umberto and {Medina-Mardones}, Anibal M. and {Tauzin}, Guillaume},
+           title = "{Persistence Steenrod modules}",
+         journal = {arXiv e-prints},
+   archivePrefix = {arXiv},
+          eprint = {1812.05031},
+    primaryClass = {math.AT},
+          adsurl = {https://ui.adsabs.harvard.edu/abs/2018arXiv181205031L},
+   }
diff --git a/docs/make.bat b/docs/make.bat
new file mode 100644
index 0000000..8084272
--- /dev/null
+++ b/docs/make.bat
@@ -0,0 +1,35 @@
+@ECHO OFF
+
+pushd %~dp0
+
+REM Command file for Sphinx documentation
+
+if "%SPHINXBUILD%" == "" (
+	set SPHINXBUILD=sphinx-build
+)
+set SOURCEDIR=.
+set BUILDDIR=_build
+
+if "%1" == "" goto help
+
+%SPHINXBUILD% >NUL 2>NUL
+if errorlevel 9009 (
+	echo.
+	echo.The 'sphinx-build' command was not found. Make sure you have Sphinx
+	echo.installed, then set the SPHINXBUILD environment variable to point
+	echo.to the full path of the 'sphinx-build' executable. Alternatively you
+	echo.may add the Sphinx directory to PATH.
+	echo.
+	echo.If you don't have Sphinx installed, grab it from
+	echo.https://www.sphinx-doc.org/
+	exit /b 1
+)
+
+%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O%
+goto end
+
+:help
+%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O%
+
+:end
+popd
diff --git a/docs/requirements.txt b/docs/requirements.txt
new file mode 100644
index 0000000..681958a
--- /dev/null
+++ b/docs/requirements.txt
@@ -0,0 +1,4 @@
+sphinx>=4.5.0
+jinja2==3.0.0
+sphinx_rtd_theme
+ipykernel
diff --git a/docs/user_guide.rst b/docs/user_guide.rst
new file mode 100644
index 0000000..9936096
--- /dev/null
+++ b/docs/user_guide.rst
@@ -0,0 +1,20 @@
+
+User Guide
+==========
+
+
+Tutorials
+---------
+
+   _notebooks/barcode.ipynb
+   _notebooks/st_barcode.ipynb
+   _notebooks/squares.ipynb
+   _notebooks/cp2.ipynb
+   _notebooks/rpn.ipynb
+
+
+Examples
+--------
+
+   _notebooks/ran/flat_Klein_bottle.ipynb
+   _notebooks/ran/cyclo-octane.ipynb
diff --git a/notebooks/barcode.ipynb b/notebooks/barcode.ipynb
index 8ca728b..9d91edc 100644
--- a/notebooks/barcode.ipynb
+++ b/notebooks/barcode.ipynb
@@ -72,7 +72,7 @@
    "source": [
     "from steenroder import barcodes\n",
     "\n",
-    "barcode, st_barcode = barcodes(0, filtration)\n",
+    "barcode, steenrod_barcode = barcodes(0, filtration)\n",
     "\n",
     "print(f'Barcode of persistent relative cohomology:')\n",
     "for d in range(maxdim + 1):\n",
@@ -212,7 +212,7 @@
    "source": [
     "#### The $R = D^\\perp V$ decomposition\n",
     "\n",
-    "To compute the matrices $R$ and $V$ in the decomposition $R = D^\\perp V$ ``steenroder`` uses the method ``get_reduced_triangular``.\n",
+    "To compute the matrices $R$ and $V$ in the decomposition $R = D^\\perp V$ ``steenroder`` uses the method ``compute_reduced_triangular``.\n",
     "It actually produces $R^\\mathrm{A}$ and $V^\\mathrm{A}$.\n",
     "More precisely, the output ``reduced`` is a tuple of ``numba.typed.List`` with ``reduced[d]`` the ``d``-dimensional part of $R^\\mathrm{A}$.\n",
     "More explicitly, for an integer $j$ we have that an integer $i$ is in the tuple ``reduced[d][j]`` if and only if $R^\\mathrm{A}_{ij} = 1$. \n",
@@ -225,9 +225,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from steenroder import get_reduced_triangular\n",
+    "from steenroder import compute_reduced_triangular\n",
     "\n",
-    "spx2idx, idxs, reduced, triangular = get_reduced_triangular(filtration_by_dim)"
+    "spx2idx, idxs, reduced, triangular = compute_reduced_triangular(filtration_by_dim)"
    ]
   },
   {
@@ -296,9 +296,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from steenroder import get_barcode_and_coho_reps\n",
+    "from steenroder import compute_barcode_and_coho_reps\n",
     "\n",
-    "barcode, coho_reps = get_barcode_and_coho_reps(idxs, reduced, triangular)"
+    "barcode, coho_reps = compute_barcode_and_coho_reps(idxs, reduced, triangular)"
    ]
   },
   {
diff --git a/notebooks/cp2.ipynb b/notebooks/cp2.ipynb
index 5aeb8be..191f043 100644
--- a/notebooks/cp2.ipynb
+++ b/notebooks/cp2.ipynb
@@ -15,31 +15,31 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from steenroder import barcodes\n",
-    "from examples import cp2, wedge, sphere, suspension, cone\n",
+    "from steenroder import barcodes, wedge, sphere, suspension, cone\n",
+    "from examples import cp2\n",
     "import numpy as np\n",
     "\n",
     "def print_rel_barcodes(k, filtration):\n",
     "    '''Auxiliary function to print outputs in the\n",
     "    persistent relative cohomology context'''\n",
-    "    barcode, st_barcode = barcodes(k, filtration)\n",
+    "    barcode, steenrod_barcode = barcodes(k, filtration)\n",
     "    print(f'*) Regular (infinite only):')\n",
     "    for d, bars in enumerate(barcode):\n",
     "        print(d, bars[np.where(bars == -1)[0],:])\n",
     "    print(f'*) Sq^{k} (non-zero only):')\n",
-    "    for d, bars in enumerate(st_barcode):\n",
+    "    for d, bars in enumerate(steenrod_barcode):\n",
     "        if bars.size > 0:\n",
     "            print(d, bars)\n",
     "            \n",
     "def print_abs_barcodes(k, filtration, length=1):\n",
     "    '''Auxiliary function to print outputs in the\n",
     "    persistent absolute cohomology context'''\n",
-    "    barcode, st_barcode = barcodes(k, filtration, absolute=True)\n",
+    "    barcode, steenrod_barcode = barcodes(k, filtration, absolute=True)\n",
     "    print(f'*) Regular (length > {length} only):')\n",
     "    for d, bars in enumerate(barcode):\n",
     "        print(d, bars[np.where(bars[:,1]-bars[:,0] > length)])\n",
     "    print(f'*) Sq^{k} (non-zero only):')\n",
-    "    for d, bars in enumerate(st_barcode):\n",
+    "    for d, bars in enumerate(steenrod_barcode):\n",
     "        if bars.size > 0:\n",
     "            print(d, bars)"
    ]
diff --git a/notebooks/cyclo-octane.ipynb b/notebooks/cyclo-octane.ipynb
index 23a1420..acbed06 100644
--- a/notebooks/cyclo-octane.ipynb
+++ b/notebooks/cyclo-octane.ipynb
@@ -480,7 +480,7 @@
     "    \"absolute\": False\n",
     "}\n",
     "\n",
-    "barcode, st_barcode = rips_barcodes(X, **barcodes_kwargs, verbose=True)"
+    "barcode, steenrod_barcode = rips_barcodes(X, **barcodes_kwargs, verbose=True)"
    ]
   },
   {
@@ -561,7 +561,7 @@
     "for dim in range(n_dims):\n",
     "    segs = []\n",
     "    multiplicities = {}\n",
-    "    dgm = st_barcode[dim]\n",
+    "    dgm = steenrod_barcode[dim]\n",
     "    for p in dgm:\n",
     "        if tuple(p) in multiplicities:\n",
     "            multiplicities[tuple(p)] += 1\n",
@@ -635,7 +635,7 @@
     "barcodes_kwargs[\"absolute\"] = True\n",
     "barcodes_kwargs[\"max_edge_length\"] = max_edge_length\n",
     "\n",
-    "barcode, st_barcode = rips_barcodes(X, **barcodes_kwargs, verbose=True)"
+    "barcode, steenrod_barcode = rips_barcodes(X, **barcodes_kwargs, verbose=True)"
    ]
   },
   {
@@ -716,7 +716,7 @@
     "for dim in range(n_dims):\n",
     "    segs = []\n",
     "    multiplicities = {}\n",
-    "    dgm = st_barcode[dim]\n",
+    "    dgm = steenrod_barcode[dim]\n",
     "    dgm = dgm[dgm[:, 1] - dgm[:, 0] > lifetime_thresh]\n",
     "    for p in dgm:\n",
     "        if tuple(p) in multiplicities:\n",
@@ -766,7 +766,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "Python 3",
    "language": "python",
    "name": "python3"
   },
@@ -780,7 +780,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.10"
+   "version": "3.9.1"
   }
  },
  "nbformat": 4,
diff --git a/notebooks/data.py b/notebooks/data.py
deleted file mode 100644
index e78935a..0000000
--- a/notebooks/data.py
+++ /dev/null
@@ -1,76 +0,0 @@
-### RP2 ###
-
-top_rp2 = (
-    (1, 2, 4), (2, 3, 4), (1, 3, 5), (2, 3, 5), (1, 4, 5),
-    (1, 2, 6), (1, 3, 6), (3, 4, 6), (2, 5, 6), (4, 5, 6)
-)
-
-### RP3 ###
-
-top_rp3 = (
-    (1, 2, 3, 4), (1, 2, 3, 5), (1, 3, 4, 5), (1, 2, 4, 6), (1, 2, 5, 6),
-    (1, 4, 5, 6), (2, 3, 4, 7), (2, 3, 5, 8), (2, 3, 7, 8), (2, 4, 6, 9),
-    (2, 4, 7, 9), (2, 5, 6, 10), (2, 5, 8, 10), (2, 6, 9, 10), (2, 7, 8, 11),
-    (2, 7, 9, 11), (2, 8, 10, 11), (2, 9, 10, 11), (3, 4, 5, 12), (3, 4, 7, 12),
-    (3, 5, 8, 12), (9, 10, 11, 12), (3, 7, 8, 13), (6, 9, 10, 13), (3, 7, 12, 13),
-    (3, 8, 12, 13), (8, 9, 12, 13), (7, 10, 12, 13), (9, 10, 12, 13), (4, 5, 6, 14),
-    (4, 6, 9, 14), (5, 6, 10, 14), (7, 8, 11, 14), (4, 5, 12, 14), (7, 8, 13, 14),
-    (6, 9, 13, 14), (8, 9, 13, 14), (6, 10, 13, 14), (7, 10, 13, 14), (4, 7, 9, 15),
-    (8, 10, 11, 15), (4, 7, 12, 15), (7, 10, 12, 15), (10, 11, 12, 15), (4, 9, 14, 15),
-    (8, 9, 14, 15), (8, 11, 14, 15), (4, 12, 14, 15), (11, 12, 14, 15), (5, 8, 10, 16),
-    (7, 9, 11, 16), (5, 8, 12, 16), (8, 9, 12, 16), (9, 11, 12, 16), (5, 10, 14, 16),
-    (7, 10, 14, 16), (7, 11, 14, 16), (5, 12, 14, 16), (11, 12, 14, 16), (7, 9, 15, 16),
-    (8, 9, 15, 16), (7, 10, 15, 16), (8, 10, 15, 16)
-)
-
-### RP4 ###
-
-top_rp4 = (
-    (1, 2, 4, 5, 11), (3, 4, 6, 7, 11), (1, 3, 8, 9, 11), (2, 6, 8, 10, 11), (5, 7, 9, 10, 11),
-    (4, 5, 7, 8, 12), (2, 4, 6, 9, 12), (1, 3, 8, 9, 12), (2, 3, 5, 10, 12), (1, 6, 7, 10, 12),
-    (2, 4, 5, 11, 12), (2, 4, 6, 11, 12), (4, 5, 7, 11, 12), (4, 6, 7, 11, 12), (2, 5, 10, 11, 12),
-    (2, 6, 10, 11, 12), (5, 7, 10, 11, 12), (6, 7, 10, 11, 12), (3, 5, 6, 8, 13), (1, 2, 7, 8, 13),
-    (2, 4, 6, 9, 13), (1, 3, 4, 10, 13), (5, 7, 9, 10, 13), (1, 2, 4, 11, 13), (1, 3, 4, 11, 13),
-    (2, 4, 6, 11, 13), (3, 4, 6, 11, 13), (1, 2, 8, 11, 13), (1, 3, 8, 11, 13), (2, 6, 8, 11, 13),
-    (3, 6, 8, 11, 13), (1, 3, 8, 12, 13), (3, 5, 8, 12, 13), (1, 7, 8, 12, 13), (5, 7, 8, 12, 13),
-    (1, 3, 10, 12, 13), (3, 5, 10, 12, 13), (1, 7, 10, 12, 13), (5, 7, 10, 12, 13), (1, 2, 4, 5, 14),
-    (3, 5, 6, 8, 14), (2, 3, 7, 9, 14), (1, 6, 7, 10, 14), (4, 8, 9, 10, 14), (3, 6, 7, 11, 14),
-    (3, 6, 8, 11, 14), (3, 7, 9, 11, 14), (3, 8, 9, 11, 14), (6, 7, 10, 11, 14), (6, 8, 10, 11, 14),
-    (7, 9, 10, 11, 14), (8, 9, 10, 11, 14), (2, 3, 5, 12, 14), (2, 4, 5, 12, 14), (3, 5, 8, 12, 14),
-    (4, 5, 8, 12, 14), (2, 3, 9, 12, 14), (2, 4, 9, 12, 14), (3, 8, 9, 12, 14), (4, 8, 9, 12, 14),
-    (1, 2, 4, 13, 14), (1, 2, 7, 13, 14), (2, 4, 9, 13, 14), (2, 7, 9, 13, 14), (1, 4, 10, 13, 14),
-    (1, 7, 10, 13, 14), (4, 9, 10, 13, 14), (7, 9, 10, 13, 14), (3, 4, 6, 7, 15), (1, 2, 7, 8, 15),
-    (1, 5, 6, 9, 15), (2, 3, 5, 10, 15), (4, 8, 9, 10, 15), (1, 2, 5, 11, 15), (1, 2, 8, 11, 15),
-    (1, 5, 9, 11, 15), (1, 8, 9, 11, 15), (2, 5, 10, 11, 15), (2, 8, 10, 11, 15), (5, 9, 10, 11, 15),
-    (8, 9, 10, 11, 15), (1, 6, 7, 12, 15), (4, 6, 7, 12, 15), (1, 7, 8, 12, 15), (4, 7, 8, 12, 15),
-    (1, 6, 9, 12, 15), (4, 6, 9, 12, 15), (1, 8, 9, 12, 15), (4, 8, 9, 12, 15), (3, 4, 6, 13, 15),
-    (3, 5, 6, 13, 15), (4, 6, 9, 13, 15), (5, 6, 9, 13, 15), (3, 4, 10, 13, 15), (3, 5, 10, 13, 15),
-    (4, 9, 10, 13, 15), (5, 9, 10, 13, 15), (1, 2, 5, 14, 15), (2, 3, 5, 14, 15), (1, 5, 6, 14, 15),
-    (3, 5, 6, 14, 15), (1, 2, 7, 14, 15), (2, 3, 7, 14, 15), (1, 6, 7, 14, 15), (3, 6, 7, 14, 15),
-    (4, 5, 7, 8, 16), (1, 5, 6, 9, 16), (2, 3, 7, 9, 16), (1, 3, 4, 10, 16), (2, 6, 8, 10, 16),
-    (1, 3, 4, 11, 16), (1, 4, 5, 11, 16), (3, 4, 7, 11, 16), (4, 5, 7, 11, 16), (1, 3, 9, 11, 16),
-    (1, 5, 9, 11, 16), (3, 7, 9, 11, 16), (5, 7, 9, 11, 16), (1, 3, 9, 12, 16), (2, 3, 9, 12, 16),
-    (1, 6, 9, 12, 16), (2, 6, 9, 12, 16), (1, 3, 10, 12, 16), (2, 3, 10, 12, 16), (1, 6, 10, 12, 16),
-    (2, 6, 10, 12, 16), (2, 6, 8, 13, 16), (5, 6, 8, 13, 16), (2, 7, 8, 13, 16), (5, 7, 8, 13, 16),
-    (2, 6, 9, 13, 16), (5, 6, 9, 13, 16), (2, 7, 9, 13, 16), (5, 7, 9, 13, 16), (1, 4, 5, 14, 16),
-    (1, 5, 6, 14, 16), (4, 5, 8, 14, 16), (5, 6, 8, 14, 16), (1, 4, 10, 14, 16), (1, 6, 10, 14, 16),
-    (4, 8, 10, 14, 16), (6, 8, 10, 14, 16), (2, 3, 7, 15, 16), (3, 4, 7, 15, 16), (2, 7, 8, 15, 16),
-    (4, 7, 8, 15, 16), (2, 3, 10, 15, 16), (3, 4, 10, 15, 16), (2, 8, 10, 15, 16), (4, 8, 10, 15, 16)
-)
-
-### CP2 ###
-
-top_cp2 = (
-    (1, 2, 4, 5, 6), (2, 3, 5, 6, 4), (3, 1, 6, 4, 5),
-    (1, 2, 4, 5, 9), (2, 3, 5, 6, 7), (3, 1, 6, 4, 8),
-    (2, 3, 6, 4, 9), (3, 1, 4, 5, 7), (1, 2, 5, 6, 8),
-    (3, 1, 5, 6, 9), (1, 2, 6, 4, 7), (2, 3, 4, 5, 8),
-    (4, 5, 7, 8, 9), (5, 6, 8, 9, 7), (6, 4, 9, 7, 8),
-    (4, 5, 7, 8, 3), (5, 6, 8, 9, 1), (6, 4, 9, 7, 2),
-    (5, 6, 9, 7, 3), (6, 4, 7, 8, 1), (4, 5, 8, 9, 2),
-    (6, 4, 8, 9, 3), (4, 5, 9, 7, 1), (5, 6, 7, 8, 2),
-    (7, 8, 1, 2, 3), (8, 9, 2, 3, 1), (9, 7, 3, 1, 2),
-    (7, 8, 1, 2, 6), (8, 9, 2, 3, 4), (9, 7, 3, 1, 5),
-    (8, 9, 3, 1, 6), (9, 7, 1, 2, 4), (7, 8, 2, 3, 5),
-    (9, 7, 2, 3, 6), (7, 8, 3, 1, 4), (8, 9, 1, 2, 5)
-)
diff --git a/notebooks/examples.py b/notebooks/examples.py
index 6d75703..6d98bdf 100644
--- a/notebooks/examples.py
+++ b/notebooks/examples.py
@@ -1,145 +1,80 @@
-from itertools import combinations, chain
-from random import shuffle
-from gudhi import SimplexTree
-from data import top_rp2, top_rp3, top_rp4, top_cp2
-
-
-def make_simplices(top_faces):
-    """It returns the simplices in a complex given its top dimensional
-    simplices only.
-
-    Parameters
-    ----------
-    top_faces : list of tuple of int
-        Represents the top faces of a simplicial complex.
-
-    Returns
-    -------
-    simplices : list of tuple of int
-        Represent all simplices in the minimal complex with the given top
-        dimensional simplices.
-
-    """
-    dim = len(top_faces[0]) - 1
-    simplices = set()
-    for top_face in top_faces:
-        for d in range(dim + 1):
-            simplices |= set(tuple(sorted(s)) for s in combinations(top_face, d + 1))
-    return list(simplices)
-
-
-def make_filtration(simplices):
-    """It returns an ordering of the simplices of a simplicial complex defining a
-    filtration of it.
-
-    Parameters
-    ----------
-    simplices : list of tuple of int
-        Represents the simplices of a simplicial complex.
-
-    Returns
-    -------
-    filtration : list of tuple of int
-        Represents a filtration for the given simplicial complex.
-
-    """
-    st = SimplexTree()
-    shuffle(simplices)
-    for spx in simplices:
-        st.insert(spx)
-    st.make_filtration_non_decreasing()
-    filtration = [tuple(pair[0]) for pair in st.get_filtration()]
-    return filtration
-
-
-def cone(simplices):
-    """Returns a filtration of the cone of a given simplicial complex.
-
-    Parameters
-    ----------
-    simplices : list of tuple of int
-        Represents the simplices of a simplicial complex.
-
-    Returns
-    -------
-    cone : list of tuple of int
-        Represents a filtration of the cone of the given simplicial complex.
-
-    """
-    vertices = set(v for v in chain.from_iterable(simplices))
-    v = max(vertices) + 1
-    new_simplices = simplices + [tuple(list(s) + [v]) for s in simplices] + [(v,)]
-    cone = make_filtration(new_simplices)
-    return cone
-
-
-def suspension(simplices):
-    """Returns a filtration of the suspension of a given simplicial complex.
-
-    Parameters
-    ----------
-    simplices : list of tuple of int
-        Represents the simplices of a simplicial complex.
-
-    Returns
-    -------
-    suspension : list of tuple of int
-        Represents a filtration of the suspension of the given simplicial complex.
-
-    """
-    vertices = set(v for v in chain.from_iterable(simplices))
-    v, w = max(vertices) + 1, max(vertices) + 2
-    new_simplices = list(simplices)
-    new_simplices += [tuple(list(s) + [v]) for s in simplices] + [(v,)] +\
-                     [tuple(list(s) + [w]) for s in simplices] + [(w,)]
-    suspension = make_filtration(new_simplices)
-    return suspension
-
-
-def sphere(n):
-    """Returns a filtration of the n-sphere.
-
-    Parameters
-    ----------
-    n : int
-        The dimension of the sphere.
-
-    Returns
-    -------
-    sphere : list of tuple of int
-        Represents a filtration of the n-sphere.
-
-    """
-    top_faces = [s for s in combinations(range(1, n + 3), n + 1)]
-    simplices = make_simplices(top_faces)
-    return make_filtration(simplices)
-
-
-def wedge(simplices_1, simplices_2):
-    """Returns a filtration wedge of two simplicial complexes.
-
-    Parameters
-    ----------
-    simplices_1 : list of tuple of int
-        Represents the first simplicial complex.
-
-    simplices_2 : list of tuple of int
-        Represents the second simplicial complex.
-
-    Returns
-    -------
-    wedge : list of tuple of int
-        Represents a filtration of the wedge of the given simplicial complexes.
-
-    """
-    vertices_1 = set(v for v in chain.from_iterable(simplices_1))
-    v = max(vertices_1)
-    simplices = list(simplices_1)
-    simplices.remove((v,))
-    for s in simplices_2:
-        simplices.append(tuple(map(lambda i: i + v - 1, s)))
-    wedge = make_filtration(simplices)
-    return wedge
+"""Example datasets"""
+
+from steenroder import make_simplices, make_filtration
+
+
+# RP2
+top_rp2 = (
+    (1, 2, 4), (2, 3, 4), (1, 3, 5), (2, 3, 5), (1, 4, 5),
+    (1, 2, 6), (1, 3, 6), (3, 4, 6), (2, 5, 6), (4, 5, 6)
+)
+
+# RP3
+top_rp3 = (
+    (1, 2, 3, 4), (1, 2, 3, 5), (1, 3, 4, 5), (1, 2, 4, 6), (1, 2, 5, 6),
+    (1, 4, 5, 6), (2, 3, 4, 7), (2, 3, 5, 8), (2, 3, 7, 8), (2, 4, 6, 9),
+    (2, 4, 7, 9), (2, 5, 6, 10), (2, 5, 8, 10), (2, 6, 9, 10), (2, 7, 8, 11),
+    (2, 7, 9, 11), (2, 8, 10, 11), (2, 9, 10, 11), (3, 4, 5, 12), (3, 4, 7, 12),
+    (3, 5, 8, 12), (9, 10, 11, 12), (3, 7, 8, 13), (6, 9, 10, 13), (3, 7, 12, 13),
+    (3, 8, 12, 13), (8, 9, 12, 13), (7, 10, 12, 13), (9, 10, 12, 13), (4, 5, 6, 14),
+    (4, 6, 9, 14), (5, 6, 10, 14), (7, 8, 11, 14), (4, 5, 12, 14), (7, 8, 13, 14),
+    (6, 9, 13, 14), (8, 9, 13, 14), (6, 10, 13, 14), (7, 10, 13, 14), (4, 7, 9, 15),
+    (8, 10, 11, 15), (4, 7, 12, 15), (7, 10, 12, 15), (10, 11, 12, 15), (4, 9, 14, 15),
+    (8, 9, 14, 15), (8, 11, 14, 15), (4, 12, 14, 15), (11, 12, 14, 15), (5, 8, 10, 16),
+    (7, 9, 11, 16), (5, 8, 12, 16), (8, 9, 12, 16), (9, 11, 12, 16), (5, 10, 14, 16),
+    (7, 10, 14, 16), (7, 11, 14, 16), (5, 12, 14, 16), (11, 12, 14, 16), (7, 9, 15, 16),
+    (8, 9, 15, 16), (7, 10, 15, 16), (8, 10, 15, 16)
+)
+
+# RP4
+top_rp4 = (
+    (1, 2, 4, 5, 11), (3, 4, 6, 7, 11), (1, 3, 8, 9, 11), (2, 6, 8, 10, 11), (5, 7, 9, 10, 11),
+    (4, 5, 7, 8, 12), (2, 4, 6, 9, 12), (1, 3, 8, 9, 12), (2, 3, 5, 10, 12), (1, 6, 7, 10, 12),
+    (2, 4, 5, 11, 12), (2, 4, 6, 11, 12), (4, 5, 7, 11, 12), (4, 6, 7, 11, 12), (2, 5, 10, 11, 12),
+    (2, 6, 10, 11, 12), (5, 7, 10, 11, 12), (6, 7, 10, 11, 12), (3, 5, 6, 8, 13), (1, 2, 7, 8, 13),
+    (2, 4, 6, 9, 13), (1, 3, 4, 10, 13), (5, 7, 9, 10, 13), (1, 2, 4, 11, 13), (1, 3, 4, 11, 13),
+    (2, 4, 6, 11, 13), (3, 4, 6, 11, 13), (1, 2, 8, 11, 13), (1, 3, 8, 11, 13), (2, 6, 8, 11, 13),
+    (3, 6, 8, 11, 13), (1, 3, 8, 12, 13), (3, 5, 8, 12, 13), (1, 7, 8, 12, 13), (5, 7, 8, 12, 13),
+    (1, 3, 10, 12, 13), (3, 5, 10, 12, 13), (1, 7, 10, 12, 13), (5, 7, 10, 12, 13), (1, 2, 4, 5, 14),
+    (3, 5, 6, 8, 14), (2, 3, 7, 9, 14), (1, 6, 7, 10, 14), (4, 8, 9, 10, 14), (3, 6, 7, 11, 14),
+    (3, 6, 8, 11, 14), (3, 7, 9, 11, 14), (3, 8, 9, 11, 14), (6, 7, 10, 11, 14), (6, 8, 10, 11, 14),
+    (7, 9, 10, 11, 14), (8, 9, 10, 11, 14), (2, 3, 5, 12, 14), (2, 4, 5, 12, 14), (3, 5, 8, 12, 14),
+    (4, 5, 8, 12, 14), (2, 3, 9, 12, 14), (2, 4, 9, 12, 14), (3, 8, 9, 12, 14), (4, 8, 9, 12, 14),
+    (1, 2, 4, 13, 14), (1, 2, 7, 13, 14), (2, 4, 9, 13, 14), (2, 7, 9, 13, 14), (1, 4, 10, 13, 14),
+    (1, 7, 10, 13, 14), (4, 9, 10, 13, 14), (7, 9, 10, 13, 14), (3, 4, 6, 7, 15), (1, 2, 7, 8, 15),
+    (1, 5, 6, 9, 15), (2, 3, 5, 10, 15), (4, 8, 9, 10, 15), (1, 2, 5, 11, 15), (1, 2, 8, 11, 15),
+    (1, 5, 9, 11, 15), (1, 8, 9, 11, 15), (2, 5, 10, 11, 15), (2, 8, 10, 11, 15), (5, 9, 10, 11, 15),
+    (8, 9, 10, 11, 15), (1, 6, 7, 12, 15), (4, 6, 7, 12, 15), (1, 7, 8, 12, 15), (4, 7, 8, 12, 15),
+    (1, 6, 9, 12, 15), (4, 6, 9, 12, 15), (1, 8, 9, 12, 15), (4, 8, 9, 12, 15), (3, 4, 6, 13, 15),
+    (3, 5, 6, 13, 15), (4, 6, 9, 13, 15), (5, 6, 9, 13, 15), (3, 4, 10, 13, 15), (3, 5, 10, 13, 15),
+    (4, 9, 10, 13, 15), (5, 9, 10, 13, 15), (1, 2, 5, 14, 15), (2, 3, 5, 14, 15), (1, 5, 6, 14, 15),
+    (3, 5, 6, 14, 15), (1, 2, 7, 14, 15), (2, 3, 7, 14, 15), (1, 6, 7, 14, 15), (3, 6, 7, 14, 15),
+    (4, 5, 7, 8, 16), (1, 5, 6, 9, 16), (2, 3, 7, 9, 16), (1, 3, 4, 10, 16), (2, 6, 8, 10, 16),
+    (1, 3, 4, 11, 16), (1, 4, 5, 11, 16), (3, 4, 7, 11, 16), (4, 5, 7, 11, 16), (1, 3, 9, 11, 16),
+    (1, 5, 9, 11, 16), (3, 7, 9, 11, 16), (5, 7, 9, 11, 16), (1, 3, 9, 12, 16), (2, 3, 9, 12, 16),
+    (1, 6, 9, 12, 16), (2, 6, 9, 12, 16), (1, 3, 10, 12, 16), (2, 3, 10, 12, 16), (1, 6, 10, 12, 16),
+    (2, 6, 10, 12, 16), (2, 6, 8, 13, 16), (5, 6, 8, 13, 16), (2, 7, 8, 13, 16), (5, 7, 8, 13, 16),
+    (2, 6, 9, 13, 16), (5, 6, 9, 13, 16), (2, 7, 9, 13, 16), (5, 7, 9, 13, 16), (1, 4, 5, 14, 16),
+    (1, 5, 6, 14, 16), (4, 5, 8, 14, 16), (5, 6, 8, 14, 16), (1, 4, 10, 14, 16), (1, 6, 10, 14, 16),
+    (4, 8, 10, 14, 16), (6, 8, 10, 14, 16), (2, 3, 7, 15, 16), (3, 4, 7, 15, 16), (2, 7, 8, 15, 16),
+    (4, 7, 8, 15, 16), (2, 3, 10, 15, 16), (3, 4, 10, 15, 16), (2, 8, 10, 15, 16), (4, 8, 10, 15, 16)
+)
+
+# CP2
+top_cp2 = (
+    (1, 2, 4, 5, 6), (2, 3, 5, 6, 4), (3, 1, 6, 4, 5),
+    (1, 2, 4, 5, 9), (2, 3, 5, 6, 7), (3, 1, 6, 4, 8),
+    (2, 3, 6, 4, 9), (3, 1, 4, 5, 7), (1, 2, 5, 6, 8),
+    (3, 1, 5, 6, 9), (1, 2, 6, 4, 7), (2, 3, 4, 5, 8),
+    (4, 5, 7, 8, 9), (5, 6, 8, 9, 7), (6, 4, 9, 7, 8),
+    (4, 5, 7, 8, 3), (5, 6, 8, 9, 1), (6, 4, 9, 7, 2),
+    (5, 6, 9, 7, 3), (6, 4, 7, 8, 1), (4, 5, 8, 9, 2),
+    (6, 4, 8, 9, 3), (4, 5, 9, 7, 1), (5, 6, 7, 8, 2),
+    (7, 8, 1, 2, 3), (8, 9, 2, 3, 1), (9, 7, 3, 1, 2),
+    (7, 8, 1, 2, 6), (8, 9, 2, 3, 4), (9, 7, 3, 1, 5),
+    (8, 9, 3, 1, 6), (9, 7, 1, 2, 4), (7, 8, 2, 3, 5),
+    (9, 7, 2, 3, 6), (7, 8, 3, 1, 4), (8, 9, 1, 2, 5)
+)
 
 
 rp2 = make_filtration(make_simplices(top_rp2))
diff --git a/notebooks/flat_Klein_bottle.ipynb b/notebooks/flat_Klein_bottle.ipynb
index 64eb6af..382ba14 100644
--- a/notebooks/flat_Klein_bottle.ipynb
+++ b/notebooks/flat_Klein_bottle.ipynb
@@ -184,7 +184,7 @@
     "    \"absolute\": False\n",
     "}\n",
     "\n",
-    "barcode, st_barcode = rips_barcodes(dm, **barcodes_kwargs, verbose=True)"
+    "barcode, steenrod_barcode = rips_barcodes(dm, **barcodes_kwargs, verbose=True)"
    ]
   },
   {
@@ -204,7 +204,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plot_diagrams(barcode, st_barcode, k=1, kind=\"R\")  # \"R\" for relative cohomology"
+    "plot_diagrams(barcode, steenrod_barcode, k=1, kind=\"R\")  # \"R\" for relative cohomology"
    ]
   },
   {
@@ -293,7 +293,7 @@
     "for dim in range(n_dims):\n",
     "    segs = []\n",
     "    multiplicities = {}\n",
-    "    dgm = st_barcode[dim]\n",
+    "    dgm = steenrod_barcode[dim]\n",
     "    dgm = dgm[dgm[:, 1] - dgm[:, 0] > lifetime_thresh]\n",
     "    for p in dgm:\n",
     "        if tuple(p) in multiplicities:\n",
@@ -359,7 +359,7 @@
     "barcodes_kwargs[\"absolute\"] = True\n",
     "barcodes_kwargs[\"max_edge_length\"] = max_edge_length\n",
     "\n",
-    "barcode, st_barcode = rips_barcodes(dm, **barcodes_kwargs, verbose=True)"
+    "barcode, steenrod_barcode = rips_barcodes(dm, **barcodes_kwargs, verbose=True)"
    ]
   },
   {
@@ -377,7 +377,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plot_diagrams(barcode[:3], st_barcode[:3], k=1, kind=\"A\")  # \"A\" for absolute cohomology"
+    "plot_diagrams(barcode[:3], steenrod_barcode[:3], k=1, kind=\"A\")  # \"A\" for absolute cohomology"
    ]
   },
   {
@@ -466,7 +466,7 @@
     "for dim in range(n_dims):\n",
     "    segs = []\n",
     "    multiplicities = {}\n",
-    "    dgm = st_barcode[dim]\n",
+    "    dgm = steenrod_barcode[dim]\n",
     "    dgm = dgm[dgm[:, 1] - dgm[:, 0] > lifetime_thresh]\n",
     "    for p in dgm:\n",
     "        if tuple(p) in multiplicities:\n",
@@ -516,7 +516,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "Python 3",
    "language": "python",
    "name": "python3"
   },
@@ -530,7 +530,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.10"
+   "version": "3.9.1"
   }
  },
  "nbformat": 4,
diff --git a/notebooks/rpn.ipynb b/notebooks/rpn.ipynb
index 44ee887..8b1c27f 100644
--- a/notebooks/rpn.ipynb
+++ b/notebooks/rpn.ipynb
@@ -15,31 +15,31 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from steenroder import barcodes\n",
-    "from examples import rp2, rp3, rp4, cone, suspension\n",
+    "from steenroder import barcodes, cone, suspension\n",
+    "from examples import rp2, rp3, rp4\n",
     "import numpy as np\n",
     "\n",
     "def print_rel_barcodes(k, filtrations, names):\n",
     "    for filtration, name in zip(filtrations, names):\n",
     "        print(f'\\n**) {name} ({len(filtration)} simplices):')\n",
-    "        barcode, st_barcode = barcodes(k, filtration)\n",
+    "        barcode, steenrod_barcode = barcodes(k, filtration)\n",
     "        print(f'*) Regular (infinite only):')\n",
     "        for d, bars in enumerate(barcode):\n",
     "            print(d, bars[np.where(bars == -1)[0],:])\n",
     "        print(f'*) Sq^{k} (non-zero only):')\n",
-    "        for d, bars in enumerate(st_barcode):\n",
+    "        for d, bars in enumerate(steenrod_barcode):\n",
     "            if bars.size > 0:\n",
     "                print(d, bars)\n",
     "                \n",
     "def print_abs_barcodes(k, filtrations, names, length=1):\n",
     "    for filtration, name in zip(filtrations, names):\n",
     "        print(f'\\n**) {name} ({len(filtration)} simplices):')\n",
-    "        barcode, st_barcode = barcodes(k, filtration, absolute=True)\n",
+    "        barcode, steenrod_barcode = barcodes(k, filtration, absolute=True)\n",
     "        print(f'*) Regular (length > {length} only):')\n",
     "        for d, bars in enumerate(barcode):\n",
     "            print(d, bars[np.where(bars[:,1]-bars[:,0] > length)])\n",
     "        print(f'*) Sq^{k} (non-zero only):')\n",
-    "        for d, bars in enumerate(st_barcode):\n",
+    "        for d, bars in enumerate(steenrod_barcode):\n",
     "            if bars.size > 0:\n",
     "                print(d, bars)"
    ]
diff --git a/notebooks/st_barcode.ipynb b/notebooks/st_barcode.ipynb
index d275b46..3ecbe20 100644
--- a/notebooks/st_barcode.ipynb
+++ b/notebooks/st_barcode.ipynb
@@ -249,11 +249,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from steenroder import sort_filtration_by_dim, get_reduced_triangular, get_barcode_and_coho_reps\n",
+    "from steenroder import sort_filtration_by_dim, compute_reduced_triangular, compute_barcode_and_coho_reps\n",
     "\n",
     "filtration_by_dim = sort_filtration_by_dim(filtration)\n",
-    "spx2idx, idxs, reduced, triangular = get_reduced_triangular(filtration_by_dim)\n",
-    "barcode, coho_reps = get_barcode_and_coho_reps(idxs, reduced, triangular)"
+    "spx2idx, idxs, reduced, triangular = compute_reduced_triangular(filtration_by_dim)\n",
+    "barcode, coho_reps = compute_barcode_and_coho_reps(idxs, reduced, triangular)"
    ]
   },
   {
@@ -262,7 +262,7 @@
    "source": [
     "#### Algorithm 1\n",
     "\n",
-    "Using an implementation of Algorithm 1, ``steenroder`` constructs the Steenrod matrix using the method ``get_steenrod_matrix``.\n",
+    "Using an implementation of Algorithm 1, ``steenroder`` constructs the Steenrod matrix using the method ``compute_steenrod_matrix``.\n",
     "Its output is a list of ``numba.typed.List``, one list per simplex dimension.\n",
     "Explicitly, ``steenrod_matrix[d][j]`` entry is the result of computing the Steenrod square of ``coho_reps[d][j]``.\n",
     "\n",
@@ -277,10 +277,10 @@
    },
    "outputs": [],
    "source": [
-    "from steenroder import get_steenrod_matrix\n",
+    "from steenroder import compute_steenrod_matrix\n",
     "\n",
     "k = 1\n",
-    "steenrod_matrix = get_steenrod_matrix(k, coho_reps, filtration_by_dim, spx2idx)\n",
+    "steenrod_matrix = compute_steenrod_matrix(k, coho_reps, filtration_by_dim, spx2idx)\n",
     "\n",
     "d = 1\n",
     "print(f'Bar     : Cocycle --> Sq^{k}-cocycle:')\n",
@@ -360,7 +360,7 @@
    "source": [
     "#### Algorithm 2\n",
     "\n",
-    "It is carried through in `get_steenrod_barcode`, whose output is the $Sq^k$-barcode of the filtration."
+    "It is carried through in `compute_steenrod_barcode`, whose output is the $Sq^k$-barcode of the filtration."
    ]
   },
   {
@@ -371,10 +371,10 @@
    },
    "outputs": [],
    "source": [
-    "from steenroder import get_steenrod_barcode\n",
+    "from steenroder import compute_steenrod_barcode\n",
     "\n",
     "dd = d + k\n",
-    "steenrod_barcode = get_steenrod_barcode(k, steenrod_matrix, idxs, reduced, barcode)\n",
+    "steenrod_barcode = compute_steenrod_barcode(k, steenrod_matrix, idxs, reduced, barcode)\n",
     "steenrod_barcode[dd]"
    ]
   },
@@ -470,7 +470,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def get_st_barcode(R, Q):\n",
+    "def get_steenrod_barcode(R, Q):\n",
     "    \"\"\"Returns the Steenrod barcodes implementing Algorithm 2\"\"\"\n",
     "    m = R.shape[1]-1\n",
     "    alive = {i: True for i in range(m)}\n",
@@ -485,7 +485,7 @@
     "    steenrod_barcode += [[-1,m-i] for i in alive if alive[i]]\n",
     "    return steenrod_barcode\n",
     "\n",
-    "get_st_barcode(R, Q)"
+    "get_steenrod_barcode(R, Q)"
    ]
   },
   {
diff --git a/requirements.txt b/requirements.txt
index 84855b7..573e731 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,5 +1,5 @@
-numpy
-numba
-psutil
-gudhi
-plotly
+numpy>=1.20.3
+numba>=0.53.0
+psutil>=5.8.0
+gudhi>=3.5.0
+plotly>=5.3.1
diff --git a/setup.py b/setup.py
index f4df3b8..af992a8 100644
--- a/setup.py
+++ b/setup.py
@@ -1,11 +1,74 @@
-from setuptools import setup
-
-setup(name="steenroder",
-      version="0.1.0",
-      description="Persistence Steenrod modules: the even prime",
-      url="http://github.com/Steenroder/steenroder",
-      author="Anibal M. Medina-Mardones and Umberto Lupo",
-      author_email="anibal.medinamardones@epfl.ch, umberto.lupo@epfl.ch",
-      license="MIT",
-      packages=["steenroder"],
-      zip_safe=False)
+#! /usr/bin/env python
+"""Computation of persistence Steenrod barcodes"""
+
+import os
+import codecs
+
+
+from pkg_resources.extern.packaging import version
+from setuptools import setup, find_packages
+
+version_file = os.path.join('steenroder', '_version.py')
+with open(version_file) as f:
+    exec(f.read())
+
+with open('requirements.txt') as f:
+    requirements = f.read().splitlines()
+
+with open('docs/requirements.txt') as f:
+    doc_requirements = f.read().splitlines()
+
+
+DISTNAME = 'steenroder'
+DESCRIPTION = 'Computation of persistence Steenrod barcodes '
+with codecs.open('README.rst', encoding='utf-8-sig') as f:
+    LONG_DESCRIPTION = f.read()
+LONG_DESCRIPTION_TYPE = 'text/x-rst'
+MAINTAINER = 'Anibal Medina Mardones, Umberto Lupo'
+MAINTAINER_EMAIL = 'steenroder@gmail.com'
+URL = 'https://github.com/Steenroder/steenroder'
+LICENSE = 'MIT'
+DOWNLOAD_URL = 'https://github.com/Steenroder/steenroder/tarball/v0.1.0'
+VERSION = __version__ # noqa
+CLASSIFIERS = ['Intended Audience :: Science/Research',
+               'Intended Audience :: Developers',
+               'License :: OSI Approved',
+               'Programming Language :: Python',
+               'Topic :: Software Development',
+               'Topic :: Scientific/Engineering',
+               'Operating System :: Microsoft :: Windows',
+               'Operating System :: POSIX',
+               'Operating System :: Unix',
+               'Operating System :: MacOS',
+               'Programming Language :: Python :: 3.8',
+               'Programming Language :: Python :: 3.9',
+               'Programming Language :: Python :: 3.10']
+KEYWORDS = 'topological data analysis, persistent ' + \
+    'homology, persistence steenrod modules, persistence ' + \
+    'steenrod barcodes'
+INSTALL_REQUIRES = requirements
+EXTRAS_REQUIRE = {
+      'tests': ['pytest',
+                'pytest-xdist',
+                'nbmake',
+                'pytest-cov',
+                'flake8'],
+      'docs': doc_requirements,
+}
+
+setup(name=DISTNAME,
+      maintainer=MAINTAINER,
+      maintainer_email=MAINTAINER_EMAIL,
+      description=DESCRIPTION,
+      license=LICENSE,
+      url=URL,
+      version=VERSION,
+      download_url=DOWNLOAD_URL,
+      long_description=LONG_DESCRIPTION,
+      long_description_content_type=LONG_DESCRIPTION_TYPE,
+      zip_safe=False,
+      classifiers=CLASSIFIERS,
+      packages=find_packages(),
+      keywords=KEYWORDS,
+      install_requires=INSTALL_REQUIRES,
+      extras_require=EXTRAS_REQUIRE)
diff --git a/steenroder/__init__.py b/steenroder/__init__.py
index b791913..719ab11 100644
--- a/steenroder/__init__.py
+++ b/steenroder/__init__.py
@@ -1,952 +1,38 @@
-import time
-from functools import lru_cache
-import psutil
-
-import numba as nb
-import numpy as np
-from numba.cpython.unsafe.tuple import tuple_setitem
-from numba.np.unsafe.ndarray import to_fixed_tuple
-
-import gudhi
-
-list_of_int64_typ = nb.types.List(nb.int64)
-int64_2d_array_typ = nb.types.Array(nb.int64, 2, "C")
-
-# Determine the number of available physical cores
-N_PHYSICAL_CORES = psutil.cpu_count(logical=False)
-
-
-def sort_filtration_by_dim(filtration, maxdim=None):
-    """Organize an input simplex-wise filtration by dimension.
-
-    Parameters
-    ----------
-    filtration : sequence of list-like of int
-        Represents a simplex-wise filtration. Entry ``i`` is a list/tuple/set
-        containing the integer indices of the vertices defining the ``i``th
-        simplex in the filtration.
-
-    maxdim : int or None, optional, default: None
-        Maximum simplex dimension to be included. ``None`` means that all
-        simplices are included.
-
-    Returns
-    -------
-    filtration_by_dim : list of list of ndarray
-        For each dimension ``d``, a list of 2 aligned int arrays: the first is
-        a 1D array containing the (ordered) positional indices of all
-        ``d``-dimensional simplices in `filtration`; the second is a 2D array
-        whose ``i``-th row is the (sorted) collection of vertices defining the
-        ``i``-th ``d``-dimensional simplex.
-
-    """
-    if maxdim is None:
-        maxdim = max(map(len, filtration)) - 1
-
-    filtration_by_dim = [[] for _ in range(maxdim + 1)]
-    for i, spx in enumerate(filtration):
-        spx_tup = tuple(sorted(spx))
-        dim = len(spx_tup) - 1
-        if dim <= maxdim:
-            filtration_by_dim[dim].append([i, spx_tup])
-
-    for dim, filtr in enumerate(filtration_by_dim):
-        filtration_by_dim[dim] = [np.asarray(x, dtype=np.int64)
-                                  for x in zip(*filtr)]
-
-    return filtration_by_dim
-
-
-@nb.njit
-def _twist_reduction(coboundary, triangular, pivots_lookup):
-    """Core of the persistent relative cohomology reduction algorithm using the
-    clearing optimization."""
-    n = len(coboundary)
-
-    rel_idxs_to_clear = []
-    for j in range(n - 1, -1, -1):
-        highest_one = coboundary[j][0] if coboundary[j] else -1
-        pivot_col = pivots_lookup[highest_one]
-        while (highest_one != -1) and (pivot_col != -1):
-            coboundary[j] = _symm_diff(coboundary[j][1:],
-                                       coboundary[pivot_col][1:])
-            triangular[j] = _symm_diff(triangular[j],
-                                       triangular[pivot_col])
-            highest_one = coboundary[j][0] if coboundary[j] else -1
-            pivot_col = pivots_lookup[highest_one]
-        if highest_one != -1:
-            pivots_lookup[highest_one] = j
-            rel_idxs_to_clear.append(highest_one)
-
-    return np.asarray(rel_idxs_to_clear, dtype=np.int64)
-
-
-@lru_cache
-def _reduce_single_dim(dim):
-    len_tups_dim = dim + 1
-    tuple_typ_dim = nb.types.UniTuple(nb.int64, len_tups_dim)
-    len_tups_next_dim = dim + 2
-
-    @nb.njit
-    def _inner_reduce_single_dim(idxs_dim, tups_dim, rel_idxs_to_clear,
-                                 idxs_next_dim=None, tups_next_dim=None):
-        """R = MV"""
-        # 1) Construct sp2idx_dim as a dict simplex: relative (i.e.
-        # in-dimension) index
-        # 2) Initialize type of reduced_dim (needed for type inference)
-        # 3) Construct triangular_dim, with entries denoting relative (i.e.
-        # in-dimension) indices
-        spx2idx_dim = nb.typed.Dict.empty(tuple_typ_dim, nb.int64)
-        reduced_dim = nb.typed.List.empty_list(list_of_int64_typ)
-        triangular_dim = nb.typed.List.empty_list(list_of_int64_typ)
-        for i in range(len(idxs_dim)):
-            spx = to_fixed_tuple(tups_dim[i], len_tups_dim)
-            spx2idx_dim[spx] = i
-            reduced_dim.append([nb.int64(x) for x in range(0)])
-            triangular_dim.append([i])
-
-        # Populate reduced_dim as the coboundary matrix and apply clearing
-        # WARNING: Column entries denote relative (i.e. in-dimension) indices!
-        if idxs_next_dim is not None:
-            for j in range(len(idxs_next_dim)):
-                spx = to_fixed_tuple(tups_next_dim[j], len_tups_next_dim)
-                for face in _drop_elements(spx):
-                    reduced_dim[spx2idx_dim[face]].append(j)
-
-            for rel_idx in rel_idxs_to_clear:
-                reduced_dim[rel_idx] = [nb.int64(x) for x in range(0)]
-
-            pivots_lookup = np.full(len(idxs_next_dim), -1, dtype=np.int64)
-
-            rel_idxs_to_clear = _twist_reduction(reduced_dim, triangular_dim,
-                                                 pivots_lookup)
-        else:
-            pivots_lookup = np.empty(0, dtype=np.int64)
-
-        return (spx2idx_dim, reduced_dim, triangular_dim,
-                rel_idxs_to_clear, pivots_lookup)
-
-    return _inner_reduce_single_dim
-
-
-@nb.njit
-def _fix_triangular_after_clearing(triangular, reduced_prev_dim,
-                                   rel_idxs_to_clear, pivots_lookup_prev_dim):
-    """Massage the V matrix to maintain the R = DV decomposition after clearing,
-    as described in https://arxiv.org/abs/1908.02518, Sec. 3.2."""
-    for rel_idx in rel_idxs_to_clear:
-        triangular[rel_idx] = reduced_prev_dim[pivots_lookup_prev_dim[rel_idx]]
-
-
-def get_reduced_triangular(filtration_by_dim):
-    """Find a full-rank upper-triangular matrix V such that R = DV is reduced,
-    where D is the anti-transpose of the filtration boundary matrix. Return both
-    R and V.
-
-    Parameters
-    ----------
-    filtration_by_dim : list of list of ndarray
-        For each dimension ``d``, a list of 2 aligned int arrays: the first is
-        a 1D array containing the (ordered) positional indices of all
-        ``d``-dimensional simplices in `filtration`; the second is a 2D array
-        whose ``i``-th row is the (sorted) collection of vertices defining the
-        ``i``-th ``d``-dimensional simplex.
-
-    Returns
-    -------
-    spx2idx : tuple of ``numba.typed.Dict``
-        One dictionary per simplex dimension. The dimension-``d`` dictionary has
-        the filtration ``d``-simplices (tuples of ints) as keys; the
-        corresponding values are the positional indices of those simplices
-        relative to the ``d``-dimensional portion of the filtration.
-
-    idxs : tuple of ndarray
-        For each dimension ``d``, this is the same as
-        ``filtration_by_dim[d][0]`` and is returned for convenience.
-
-    reduced : tuple of ``numba.typed.List``
-        One list of int per simplex dimension. ``reduced[d]`` is the
-        ``d``-dimensional part of the "R" matrix in R = DV. In the computation,
-        ``reduced[d][i]`` is initialized as the coboundary of the ``i``th input
-        simplex in ``filtration_by_dim[d][1]``, i.e. as the sorted list of
-        positional indices (relative to the ``(d+1)``-dimensional portion of the
-        filtration) of that simplex's cofacets.
-
-    triangular : tuple of ``numba.typed.List``
-        On list of int per simplex dimension. ``triangular[d]`` is the
-        ``d``-dimensional part of the "V" matrix in R = DV. In the computation,
-        ``triangular[d][i]`` is initialized as the singleton list ``[i]``.
-
-    """
-    maxdim = len(filtration_by_dim) - 1
-    spx2idx_idxs_reduced_triangular = []
-    # Initialize relative (i.e. in-dimension) indices to clear, as an empty
-    # int array in dim 0
-    rel_idxs_to_clear = np.empty(0, dtype=np.int64)
-    reduced_prev_dim = nb.typed.List.empty_list(list_of_int64_typ)
-    pivots_lookup_prev_dim = np.empty(0, dtype=np.int64)
-    for dim in range(maxdim):
-        reduction_dim = _reduce_single_dim(dim)
-        idxs_dim, tups_dim = filtration_by_dim[dim]
-        idxs_next_dim, tups_next_dim = filtration_by_dim[dim + 1]
-        (spx2idx_dim, reduced_dim, triangular_dim,
-         rel_idxs_to_clear_next_dim, pivots_lookup) = \
-            reduction_dim(idxs_dim,
-                          tups_dim,
-                          rel_idxs_to_clear,
-                          idxs_next_dim=idxs_next_dim,
-                          tups_next_dim=tups_next_dim)
-        _fix_triangular_after_clearing(triangular_dim,
-                                       reduced_prev_dim,
-                                       rel_idxs_to_clear,
-                                       pivots_lookup_prev_dim)
-        spx2idx_idxs_reduced_triangular.append((spx2idx_dim,
-                                                idxs_dim,
-                                                reduced_dim,
-                                                triangular_dim))
-        rel_idxs_to_clear = rel_idxs_to_clear_next_dim
-        reduced_prev_dim = reduced_dim
-        pivots_lookup_prev_dim = pivots_lookup
-
-    reduction_dim = _reduce_single_dim(maxdim)
-    idxs_dim, tups_dim = filtration_by_dim[maxdim]
-    spx2idx_dim, reduced_dim, triangular_dim, _, _ = \
-        reduction_dim(idxs_dim, tups_dim, rel_idxs_to_clear)
-    _fix_triangular_after_clearing(triangular_dim,
-                                   reduced_prev_dim,
-                                   rel_idxs_to_clear,
-                                   pivots_lookup_prev_dim)
-    spx2idx_idxs_reduced_triangular.append((spx2idx_dim,
-                                            idxs_dim,
-                                            reduced_dim,
-                                            triangular_dim))
-
-    return tuple(zip(*spx2idx_idxs_reduced_triangular))
-
-
-@nb.njit
-def get_barcode_and_coho_reps(idxs, reduced, triangular,
-                              filtration_values=None):
-    """Extract the ordinary persistent relative cohomology barcode as well as
-    one persistent relative cohomology representative per bar.
-
-    Parameters
-    ----------
-    idxs : tuple of ndarray
-        For each dimension ``d``, a 1D int array containing the (ordered)
-        positional indices of all ``d``-dimensional simplices in the filtration.
-
-    reduced : tuple of ``numba.typed.List``
-        One list of int per simplex dimension, representing the
-        ``d``-dimensional part of the "R" matrix in R = DV. In the same format
-        as returned by `get_reduced_triangular`.
-
-    triangular : tuple of ``numba.typed.List``
-        One list of int per simplex dimension, representing the
-        ``d``-dimensional part of the "V" matrix in R = DV. In the same format
-        as returned by `get_reduced_triangular`.
-
-    filtration_values : ndarray or None, optional, default: None
-        Optionally, a single 1D array of filtration values for each simplex in
-        the filtration and in all dimensions contained in `idxs`, `reduced` and
-        `triangular`. Bars with equal birth and death filtration values are
-        discarded.
-
-    Returns
-    -------
-    barcode : list of ndarray
-        For each dimension ``d``, a 2D int array of shape ``(n_bars, 2)``
-        containing the birth (entry 1) and death (entry 0) indices of persistent
-        relative homology classes in degree ``d``. See `filtration_values`.
-        Essential bars are represented by pairs with death equal to ``-1``. Bars
-        are sorted in order of decreasing birth indices.
-
-    coho_reps : list of ``numba.typed.List``
-        For each dimension ``d``, a list of representatives of persistent
-        relative cohomology classes in degree ``d``. Each such representative is
-        represented as a list of positional indices relative to the
-        ``d``-dimensional portion of the filtration. ``coho_reps[d][j]``
-        corresponds to ``barcode[d][j]``.
-
-    """
-    barcode = []
-    coho_reps = []
-
-    if filtration_values is None:
-        pairs_0 = []
-        coho_reps_0 = []
-        for i in range(len(idxs[0])):
-            if not reduced[0][i]:
-                pairs_0.append([-1, idxs[0][i]])
-                coho_reps_0.append(triangular[0][i])
-        pairs_0 = np.asarray(pairs_0)
-        lexsrt = _lexsort_barcode(pairs_0)
-        barcode.append(pairs_0[lexsrt])
-        coho_reps.append(nb.typed.List([coho_reps_0[k] for k in lexsrt]))
-
-        for dim in range(1, len(idxs)):
-            all_birth_indices = set()
-            pairs_dim = []
-            coho_reps_dim = []
-            for i in range(len(idxs[dim - 1])):
-                if reduced[dim - 1][i]:
-                    b = idxs[dim][reduced[dim - 1][i][0]]
-                    d = idxs[dim - 1][i]
-                    pairs_dim.append([d, b])
-                    coho_reps_dim.append(reduced[dim - 1][i])
-                    all_birth_indices.add(b)
-
-            for i in range(len(idxs[dim])):
-                if idxs[dim][i] not in all_birth_indices:
-                    if not reduced[dim][i]:
-                        pairs_dim.append([-1, idxs[dim][i]])
-                        coho_reps_dim.append(triangular[dim][i])
-
-            if not len(pairs_dim):
-                pairs_dim = np.empty((0, 2), dtype=np.int64)
-            else:
-                pairs_dim = np.asarray(pairs_dim)
-            lexsrt = _lexsort_barcode(pairs_dim)
-            barcode.append(pairs_dim[lexsrt])
-            coho_reps.append(nb.typed.List([coho_reps_dim[k] for k in lexsrt]))
-
-    else:
-        pairs_0 = []
-        coho_reps_0 = []
-        for i in range(len(idxs[0])):
-            if not reduced[0][i]:
-                pairs_0.append([-1, idxs[0][i]])
-                coho_reps_0.append(triangular[0][i])
-        pairs_0 = np.asarray(pairs_0)
-        lexsrt = _lexsort_barcode(pairs_0)
-        barcode.append(pairs_0[lexsrt])
-        coho_reps.append(nb.typed.List([coho_reps_0[k] for k in lexsrt]))
-
-        for dim in range(1, len(idxs)):
-            all_birth_indices = set()
-            pairs_dim = []
-            coho_reps_dim = []
-            for i in range(len(idxs[dim - 1])):
-                if reduced[dim - 1][i]:
-                    b = idxs[dim][reduced[dim - 1][i][0]]
-                    d = idxs[dim - 1][i]
-                    if filtration_values[b] != filtration_values[d]:
-                        pairs_dim.append([d, b])
-                        coho_reps_dim.append(reduced[dim - 1][i])
-                    all_birth_indices.add(b)
-
-            for i in range(len(idxs[dim])):
-                if idxs[dim][i] not in all_birth_indices:
-                    if not reduced[dim][i]:
-                        pairs_dim.append([-1, idxs[dim][i]])
-                        coho_reps_dim.append(triangular[dim][i])
-
-            if not len(pairs_dim):
-                pairs_dim = np.empty((0, 2), dtype=np.int64)
-            else:
-                pairs_dim = np.asarray(pairs_dim)
-            lexsrt = _lexsort_barcode(pairs_dim)
-            barcode.append(pairs_dim[lexsrt])
-            coho_reps.append(nb.typed.List([coho_reps_dim[k] for k in lexsrt]))
-
-    return barcode, coho_reps
-
-
-@nb.njit
-def _initialize_steenrod_matrix(num_dimensions):
-    return [nb.typed.List.empty_list(list_of_int64_typ)
-            for _ in range(num_dimensions)]
-
-
-@lru_cache
-def _populate_steenrod_matrix_single_dim(dim_plus_k):
-    length = dim_plus_k + 1
-
-    @nb.njit(parallel=True)
-    def _inner(coho_reps_dim, tups_dim, spx2idx_dim_plus_k, n_jobs=-1):
-        steenrod_matrix_dim_plus_k = \
-            nb.typed.List([[nb.int64(0) for _ in range(0)]
-                           for _ in coho_reps_dim])
-
-        if n_jobs == -1:
-            n_jobs = N_PHYSICAL_CORES
-
-        for job_idx in nb.prange(n_jobs):
-            for coho_reps_dim_idx in range(job_idx, len(coho_reps_dim), n_jobs):
-                rep = coho_reps_dim[coho_reps_dim_idx]
-                cocycle = tups_dim[np.asarray(rep)]
-
-                # STSQ
-                cochain = set(
-                    [to_fixed_tuple(np.empty(length, dtype=np.int64), length)
-                     for _ in range(0)]
-                )
-                for i in range(len(cocycle)):
-                    for j in range(i + 1, len(cocycle)):
-                        a, b = set(cocycle[i]), set(cocycle[j])
-                        u = a.union(b)
-                        if len(u) == length:
-                            u_tuple = to_fixed_tuple(np.asarray(sorted(u)),
-                                                     length)
-                            if u_tuple in spx2idx_dim_plus_k:
-                                a_bar, b_bar = a.difference(b), b.difference(a)
-                                u_bar = sorted(a_bar.union(b_bar))
-                                index = {}
-                                for v in a_bar.union(b_bar):
-                                    pos = u_tuple.index(v)
-                                    pos_bar = u_bar.index(v)
-                                    index[v] = (pos + pos_bar) % 2
-                                index_a = set()
-                                index_b = set()
-                                for v in a_bar:
-                                    index_a.add(index[v])
-                                for w in b_bar:
-                                    index_b.add(index[w])
-                                if (index_a == set([0])
-                                    and index_b == set([1])) \
-                                        or (index_a == set([1])
-                                            and index_b == set([0])):
-                                    cochain ^= {u_tuple}
-
-                steenrod_matrix_dim_plus_k[coho_reps_dim_idx] = \
-                    sorted([spx2idx_dim_plus_k[spx] for spx in cochain])
-
-        return steenrod_matrix_dim_plus_k
-
-    return _inner
-
-
-def get_steenrod_matrix(k, coho_reps, filtration_by_dim, spx2idx, n_jobs=-1):
-    """Compute the Steenrod matrices in each dimension.
-
-    Parameters
-    ----------
-    k : int
-        Positive integer defining the cohomology operation Sq^k to be performed.
-
-    coho_reps : list of ``numba.typed.List``
-        For each dimension ``d``, a list of representatives of persistent
-        relative cohomology classes in degree ``d``. In the same format as
-        returned by `get_barcode_and_coho_reps`.
-
-    filtration_by_dim : list of list of ndarray
-        For each dimension ``d``, a list of 2 aligned int arrays: the first is
-        a 1D array containing the (ordered) positional indices of all
-        ``d``-dimensional simplices in `filtration`; the second is a 2D array
-        whose ``i``-th row is the (sorted) collection of vertices defining the
-        ``i``-th ``d``-dimensional simplex.
-
-    spx2idx : tuple of ``numba.typed.Dict``
-        One dictionary per simplex dimension. The dimension-``d`` dictionary has
-        the filtration ``d``-simplices (tuples of ints) as keys; the
-        corresponding values are the positional indices of those simplices
-        relative to the ``d``-dimensional portion of the filtration.
-
-    n_jobs : int, optional, default: ``-1``
-        [Experimental] Controls the number of threads to be used during parallel
-        computation of the Steenrod squares. ``-1`` means using all available
-        physical cores.
-
-    Returns
-    -------
-    steenrod_matrix : list of ``numba.typed.List``
-        One list per simplex dimension. ``steenrod_matrix[d][j]`` is the result
-        of computing the Steenrod square of ``coho_reps[d - k][j]``.
-
-    """
-    steenrod_matrix = _initialize_steenrod_matrix(k)
-
-    for dim, coho_reps_dim in enumerate(coho_reps[:-k]):
-        dim_plus_k = dim + k
-        tups_dim = filtration_by_dim[dim][1]
-        spx2idx_dim_plus_k = spx2idx[dim + k]
-        populate_steenrod_matrix_single_dim = \
-            _populate_steenrod_matrix_single_dim(dim_plus_k)
-        steenrod_matrix_dim_plus_k = populate_steenrod_matrix_single_dim(
-            coho_reps_dim, tups_dim, spx2idx_dim_plus_k, n_jobs=n_jobs
-        )
-        steenrod_matrix.append(steenrod_matrix_dim_plus_k)
-
-    return steenrod_matrix
-
-
-@nb.njit
-def _steenrod_barcode_single_dim(steenrod_matrix_dim, n_idxs_dim, idxs_prev_dim,
-                                 reduced_prev_dim, births_dim_minus_k):
-    # Construct augmented matrix
-    augmented = []
-    for i in range(len(reduced_prev_dim)):
-        augmented.append([nb.int64(x) for x in reduced_prev_dim[i]])
-    for i in range(len(steenrod_matrix_dim)):
-        augmented.append([nb.int64(x) for x in steenrod_matrix_dim[i]])
-
-    n_cols_red = len(reduced_prev_dim)
-    n_cols_st = len(steenrod_matrix_dim)
-    n_births_dim_minus_k = len(births_dim_minus_k)
-
-    pivots_lookup = np.full(n_idxs_dim, -1, dtype=np.int64)
-    alive = np.ones(n_births_dim_minus_k, dtype=np.bool_)
-    st_barcode_dim = []
-
-    for n_cols_st_curr, birth in enumerate(births_dim_minus_k):
-        if idxs_prev_dim[-1] >= birth:
-            break
-    for i in np.arange(n_cols_red)[::-1]:
-        idx = idxs_prev_dim[i]
-        next_birth = births_dim_minus_k[n_cols_st_curr]
-        is_idx_geq = idx >= next_birth
-        if is_idx_geq:
-            if i:
-                next_idx = idxs_prev_dim[i - 1]
-                for jj in range(n_cols_st_curr, n_cols_st):
-                    if next_idx >= births_dim_minus_k[jj]:
-                        break
-                    n_cols_st_curr += 1
-            else:
-                n_cols_st_curr = n_cols_st
-
-        if augmented[i]:
-            pivots_lookup[augmented[i][0]] = i
-
-        pivot_col_idxs_from_st = []
-        for ii in range(n_cols_red, n_cols_red + n_cols_st_curr):
-            highest_one = augmented[ii][0] if augmented[ii] else -1
-            pivot_col = pivots_lookup[highest_one]
-            while (highest_one != -1) and (pivot_col != -1):
-                augmented[ii] = _symm_diff(augmented[ii][1:],
-                                           augmented[pivot_col][1:])
-                highest_one = augmented[ii][0] if augmented[ii] else -1
-                pivot_col = pivots_lookup[highest_one]
-            if highest_one != -1:
-                pivots_lookup[highest_one] = ii
-                # Record pivot indices coming from Steenrod part of augmented
-                pivot_col_idxs_from_st.append(highest_one)
-            elif alive[ii - n_cols_red]:
-                alive[ii - n_cols_red] = False
-                birth = births_dim_minus_k[ii - n_cols_red]
-                if idx < birth:
-                    st_barcode_dim.append([idx, birth])
-
-        # Reset pivots_lookup for next iteration
-        for col_idx in pivot_col_idxs_from_st:
-            pivots_lookup[col_idx] = -1
-
-    for i in range(n_births_dim_minus_k):
-        if alive[i]:
-            st_barcode_dim.append([-1, births_dim_minus_k[i]])
-
-    return st_barcode_dim
-
-
-def get_steenrod_barcode(k, steenrod_matrix, idxs, reduced, barcode,
-                         filtration_values=None):
-    """Compute the (relative) Steenrod barcodes.
-
-    Parameters
-    ----------
-    k : int
-        Positive integer defining the cohomology operation Sq^k to be performed.
-
-    steenrod_matrix : list of ``numba.typed.List``
-        One list per simplex dimension. ``steenrod_matrix[d][j]`` is the result
-        of computing the Steenrod square of the ``j``th latest (by birth)
-        persistent relative cohomology representative in degree ``d`` (and this
-        representative must represent bar ``barcode[d][j]``). See
-        `get_steenrod_matrix`.
-
-    idxs : tuple of ndarray
-        For each dimension ``d``, a 1D int array containing the (ordered)
-        positional indices of all ``d``-dimensional simplices in the filtration.
-
-    reduced : tuple of ``numba.typed.List``
-        One list of int per simplex dimension, representing the
-        ``d``-dimensional part of the "R" matrix in R = DV. In the same format
-        as returned by `get_reduced_triangular`.
-
-    barcode : list of ndarray
-        For each dimension ``d``, a 2D int array of shape ``(n_bars, 2)``
-        containing the birth (entry 1) and death (entry 0) indices of persistent
-        relative homology classes in degree ``d``. Essential bars must be
-        represented by pairs with death equal to ``-1``. Bars must be sorted in
-        order of decreasing birth indices.
-
-    filtration_values : ndarray or None, optional, default: None
-        Optionally, a single 1D array of filtration values for each simplex in
-        the filtration and in all dimensions contained in `idxs` and `reduced`.
-        Steenrod bars with equal birth and death filtration values are
-        discarded.
-
-    Returns
-    -------
-    st_barcode : list of ndarray
-        The (relative) Sq^k-barcode. For each dimension ``d``, a 2D int array
-        of shape ``(n_bars, 2)`` containing the birth (entry 1) and death (entry
-        0) indices of Steenrod bars. Essential Steenrod bars are represented by
-        pairs with death equal to ``-1``.
-
-    """
-    def nontrivial_bars(barcode_dim):
-        infinite_bars = barcode_dim[:, 0] == -1
-        return np.logical_or(
-            infinite_bars,
-            np.logical_and(np.logical_not(infinite_bars),
-                           (filtration_values[barcode_dim[:, 0]] !=
-                            filtration_values[barcode_dim[:, 1]]))
-        )
-
-    st_barcode = [np.empty((0, 2), dtype=np.int64) for _ in range(k)]
-    for dim in range(k, len(steenrod_matrix)):
-        births_dim_minus_k = barcode[dim - k][:, 1]
-        idxs_dim = idxs[dim]
-        idxs_prev_dim = idxs[dim - 1]
-        reduced_prev_dim = reduced[dim - 1]
-        st_barcode_dim = _steenrod_barcode_single_dim(steenrod_matrix[dim],
-                                                      len(idxs_dim),
-                                                      idxs_prev_dim,
-                                                      reduced_prev_dim,
-                                                      births_dim_minus_k)
-        # NB: Conversion to array must happen outside jitted code due to
-        # https://github.com/numba/numba/issues/3579
-        st_barcode_dim = \
-            np.asarray(st_barcode_dim, dtype=np.int64).reshape((-1, 2))
-        if filtration_values is not None:
-            nontrivial_mask = nontrivial_bars(st_barcode_dim)
-            st_barcode_dim = st_barcode_dim[nontrivial_mask]
-        st_barcode.append(st_barcode_dim)
-
-    return st_barcode
-
-
-def barcodes(
-        k, filtration, absolute=False, filtration_values=None,
-        return_filtration_values=False, maxdim=None, n_jobs=1, verbose=False
-):
-    """Given a filtration, compute ordinary persistent (relative or absolute)
-    cohomology barcodes and relative Steenrod barcodes.
-
-    Parameters
-    ----------
-    k : int
-        Positive integer defining the cohomology operation Sq^k to be performed.
-
-    filtration : sequence of list-like of int
-        Represents a simplex-wise filtration. Entry ``i`` is a list/tuple/set
-        containing the integer indices of the vertices defining the ``i``th
-        simplex in the filtration.
-
-    absolute : bool, optional, default: ``False``
-        If ``True``, return the ordinary persistent absolute homology barcode,
-        and move inessential relative Steenrod bars to one degree lower while
-        keeping essential bars in their degree. If ``False``, return the
-        ordinary persistent relative cohomology barcode and relative Steenrod
-        barcode.
-
-    filtration_values : ndarray or None, optional, default: None
-        Optionally, a single 1D array of filtration values for each simplex in
-        the filtration. Ordinary and Steenrod bars with equal birth and death
-        filtration values are discarded by the computation.
-
-    return_filtration_values : bool, optional, default: ``False``
-        If ``True``, birth and deaths will be expressed as filtration values
-        instead of filtration indices. Ignored if `filtration_values` is
-        ``None``.
-
-    maxdim : int or None, optional, default: None
-        Maximum simplex dimension to be included. ``None`` means that all
-        simplices are included.
-
-    n_jobs : int, optional, default: ``1``
-        [Experimental] Controls the number of threads to be used during parallel
-        computation of the Steenrod squares. ``-1`` means using all available
-        physical cores.
-
-    verbose : bool, optional, default: ``False``
-        Whether to print timings for the intermediate steps in the computation.
-
-    Returns
-    -------
-    barcode : list of ndarray
-        For each dimension ``d``, a 2D int or float array of shape ``(n_bars,
-        2)`` containing the births and deaths of persistent relative cohomology
-        classes in degree ``d``. If `absolute` is ``False``, the birth of a bar
-        is in entry 1 and the death in entry 0; otherwise, the positions are
-        reversed. Births and death are expressed either as global filtration
-        indices or as filtration values depending on `filtration_values` and
-        `return_filtration_values`. If they are expressed as indices, essential
-        bars have death equal to ``-1``; otherwise, essential bars have death
-        equal to ``numpy.inf``.
-
-    st_barcode : list of ndarray
-        The (relative) Sq^k-barcode. For each dimension ``d``, a 2D int or float
-        array of shape ``(n_bars, 2)`` containing the birth (entry 1) and death
-        (entry 0) indices of Steenrod bars. The same conventions as for
-        `barcode` are used for birth and death values.
-
-    """
-    if verbose:
-        tic = time.time()
-    filtration_by_dim = sort_filtration_by_dim(filtration, maxdim=maxdim)
-    spx2idx, idxs, reduced, triangular = \
-        get_reduced_triangular(filtration_by_dim)
-    barcode, coho_reps = \
-        get_barcode_and_coho_reps(idxs, reduced, triangular,
-                                  filtration_values=filtration_values)
-
-    if k == 0:  # Sq^0 is the identity
-        return barcode, barcode
-
-    if k < 0:  # Sq^k is the zero map
-        return barcode, list()
-
-    if verbose:
-        toc = time.time()
-        print(f"Usual barcode computed, time taken: {(toc - tic):.3f} s")
-        tic = time.time()
-    steenrod_matrix = get_steenrod_matrix(k, coho_reps, filtration_by_dim,
-                                          spx2idx, n_jobs=n_jobs)
-    if verbose:
-        toc = time.time()
-        print(f"Steenrod matrix computed, time taken: {(toc - tic):.3f} s")
-        tic = time.time()
-    st_barcode = get_steenrod_barcode(k, steenrod_matrix, idxs, reduced,
-                                      barcode,
-                                      filtration_values=filtration_values)
-    if verbose:
-        toc = time.time()
-        print(f"Steenrod barcode computed, time taken: {(toc - tic):.3f} s")
-
-    if absolute:
-        barcode = _to_absolute_barcode(
-            barcode, filtration_values=filtration_values,
-            return_filtration_values=return_filtration_values
-        )
-        st_barcode = _to_absolute_barcode(
-            st_barcode, filtration_values=filtration_values,
-            return_filtration_values=return_filtration_values
-        )
-
-        return barcode, st_barcode
-
-    elif return_filtration_values and (filtration_values is not None):
-        barcode = _to_values_barcode(barcode, filtration_values)
-        st_barcode = _to_values_barcode(st_barcode, filtration_values)
-
-    return barcode, st_barcode
-
-
-def _to_absolute_barcode(rel_barcode, filtration_values=None,
-                         return_filtration_values=True):
-    abs_barcode = []
-
-    if (not return_filtration_values) or (filtration_values is None):
-        dtype = np.int64
-        for dim, rel_barcode_dim in enumerate(rel_barcode):
-            abs_barcode_dim = []
-            for pair in rel_barcode_dim:
-                if pair[0] == -1:
-                    abs_barcode_dim.append((pair[1], -1))
-                else:
-                    abs_barcode[dim - 1].append((pair[0], pair[1]))
-            abs_barcode.append(abs_barcode_dim)
-    else:
-        dtype = filtration_values.dtype
-        for dim, rel_barcode_dim in enumerate(rel_barcode):
-            abs_barcode_dim = []
-            for pair in rel_barcode_dim:
-                if pair[0] == -1:
-                    abs_barcode_dim.append(
-                        (filtration_values[pair[1]], np.inf)
-                    )
-                else:
-                    abs_barcode[dim - 1].append(
-                        (filtration_values[pair[0]], filtration_values[pair[1]])
-                    )
-            abs_barcode.append(abs_barcode_dim)
-
-    return [np.array(abs_barcode_dim, dtype=dtype).reshape(-1, 2)
-            for abs_barcode_dim in abs_barcode]
-
-
-def _to_values_barcode(barcode, filtration_values):
-    values_barcode = []
-    for dim, barcode_dim in enumerate(barcode):
-        values_barcode_dim = []
-        for pair in barcode_dim:
-            if pair[0] == -1:
-                values_barcode_dim.append(
-                    (-np.inf, filtration_values[pair[1]])
-                )
-            else:
-                values_barcode_dim.append(
-                    (filtration_values[pair[0]], filtration_values[pair[1]])
-                )
-        values_barcode.append(
-            np.array(values_barcode_dim,
-                     dtype=filtration_values.dtype).reshape(-1, 2)
-        )
-
-    return values_barcode
-
-
-@nb.njit
-def _symm_diff(x, y):
-    n = len(x)
-    m = len(y)
-    result = []
-    i = 0
-    j = 0
-    while (i < n) and (j < m):
-        if x[i] < y[j]:
-            result.append(x[i])
-            i += 1
-        elif y[j] < x[i]:
-            result.append(y[j])
-            j += 1
-        else:
-            i += 1
-            j += 1
-
-    while i < n:
-        result.append(x[i])
-        i += 1
-
-    while j < m:
-        result.append(y[j])
-        j += 1
-
-    return result
-
-
-@nb.njit
-def _lexsort_barcode(arr):
-    return np.argsort(arr[:, 1])[::-1]
-
-
-@nb.njit
-def _drop_elements(tup: tuple):
-    for x in range(len(tup)):
-        empty = tup[:-1]  # Not empty, but the right size and will be mutated
-        idx = 0
-        for i in range(len(tup)):
-            if i != x:
-                empty = tuple_setitem(empty, idx, tup[i])
-                idx += 1
-        yield empty
-
-
-def rips_barcodes(X, max_edge_length=np.inf, distance_matrix=False,
-                  k=1, max_simplex_dimension=2, absolute=False, n_jobs=1,
-                  verbose=False):
-    """Construct a Rips filtration and compute ordinary persistent (relative or
-    absolute) cohomology barcodes and relative Steenrod barcodes.
-
-    Parameters
-    ----------
-    X : ndarray of shape (n_samples, n_features) or (n_samples, n_samples)
-        Input point cloud (if `distance_matrix` is ``False``) or distance matrix
-        (if `distance_matrix` is ``True``).
-
-    max_edge_length : float, optional, default: ``np.inf``
-        Vietoris–Rips filtration threshold.
-
-    distance_matrix : bool, optional, default: ``False``
-        Whether the input `X` is a distance matrix or a point cloud.
-
-    k : int
-        Positive integer defining the cohomology operation Sq^k to be performed.
-
-    max_simplex_dimension : int, optional, default: ``2``
-        Maximum simplex dimension to be included in the Rips filtration.
-
-    absolute : bool, optional, default: ``False``
-        If ``True``, return the ordinary persistent absolute homology barcode,
-        and move inessential relative Steenrod bars to one degree lower while
-        keeping essential bars in their degree. If ``False``, return the
-        ordinary persistent relative cohomology barcode and relative Steenrod
-        barcode.
-
-    n_jobs : int, optional, default: ``1``
-        [Experimental] Controls the number of threads to be used during parallel
-        computation of the Steenrod squares. ``-1`` means using all available
-        physical cores.
-
-    verbose : bool, optional, default: ``False``
-        Whether to print information about the intermediate steps in the
-        computation.
-
-    Returns
-    -------
-    barcode : list of ndarray
-        For each dimension ``d``, a 2D int or float array of shape ``(n_bars,
-        2)`` containing the births and deaths of persistent relative cohomology
-        classes in degree ``d``. If `absolute` is ``False``, the birth of a bar
-        is in entry 1 and the death in entry 0; otherwise, the positions are
-        reversed. Births and death are expressed either as global filtration
-        indices or as filtration values depending on `filtration_values` and
-        `return_filtration_values`. If they are expressed as indices, essential
-        bars have death equal to ``-1``; otherwise, essential bars have death
-        equal to ``numpy.inf``.
-
-    st_barcode : list of ndarray
-        The (relative) Sq^k-barcode. For each dimension ``d``, a 2D int or float
-        array of shape ``(n_bars, 2)`` containing the birth (entry 1) and death
-        (entry 0) indices of Steenrod bars. The same conventions as for
-        `barcode` are used for birth and death values.
-
-    """
-    spx_tree_kwargs = {
-        "distance_matrix" if distance_matrix else "points": X,
-        "max_edge_length": max_edge_length
-    }
-    spx_tree = gudhi.RipsComplex(**spx_tree_kwargs).\
-        create_simplex_tree(max_dimension=1)
-
-    len_filtration = None
-    if verbose:
-        for i, _ in enumerate(spx_tree.get_filtration()):
-            pass
-        len_filtration = i + 1
-        print(f"Filtration with {len_filtration} simplices up to dimension 1 "
-              "initially")
-
-        print("Running Edge Collapse...")
-
-    it = 0
-    while True:
-        spx_tree.collapse_edges()
-
-        for i, _ in enumerate(spx_tree.get_filtration()):
-            pass
-
-        if i + 1 == len_filtration:
-            break
-        else:
-            len_filtration = i + 1
-            it += 1
-            if verbose:
-                print(f"  EC iteration {it}: {len_filtration} simplices up to "
-                      "dimension 1")
-
-    # Construct the Rips filtration and run Steenroder
-    spx_tree.expansion(max_simplex_dimension)
-    if verbose:
-        for i, _ in enumerate(spx_tree.get_filtration()):
-            pass
-        print(f"There are {i} simplices up to dimension "
-              f"{max_simplex_dimension} after EC.")
-
-    filtration, filtration_values = \
-        zip(*((tuple(t[0]), t[1]) for t in spx_tree.get_filtration()))
-    filtration_values = np.asarray(filtration_values, dtype=np.float32)
-
-    barcode, st_barcode = barcodes(
-        k,
-        filtration,
-        filtration_values=filtration_values,
-        absolute=absolute,
-        return_filtration_values=True,
-        n_jobs=n_jobs,
-        verbose=verbose
-    )
-
-    return barcode, st_barcode
+from ._version import __version__
+
+from .core import (
+    compute_reduced_triangular,
+    compute_barcode_and_coho_reps,
+    compute_steenrod_matrix,
+    compute_steenrod_barcode,
+    barcodes,
+    rips_barcodes,
+)
+from .preprocessing import (
+    make_simplices,
+    make_filtration,
+    sort_filtration_by_dim,
+    cone,
+    suspension,
+    sphere,
+    wedge,
+)
+from .plotting import plot_diagrams
+
+
+__all__ = [
+    "compute_reduced_triangular",
+    "compute_barcode_and_coho_reps",
+    "compute_steenrod_matrix",
+    "compute_steenrod_barcode",
+    "barcodes",
+    "rips_barcodes",
+    "make_simplices",
+    "make_filtration",
+    "sort_filtration_by_dim",
+    "cone",
+    "suspension",
+    "sphere",
+    "wedge",
+    "plot_diagrams",
+]
diff --git a/steenroder/_version.py b/steenroder/_version.py
new file mode 100644
index 0000000..f7ac8c9
--- /dev/null
+++ b/steenroder/_version.py
@@ -0,0 +1,21 @@
+"""``steenroder`` is a python API for the computation of 
+persistent steenrod barcodes."""
+
+# PEP0440 compatible formatted version, see:
+# https://www.python.org/dev/peps/pep-0440/
+#
+# Generic release markers:
+# X.Y
+# X.Y.Z # For bugfix releases
+#
+# Admissible pre-release markers:
+# X.YaN # Alpha release
+# X.YbN # Beta release
+# X.YrcN # Release Candidate
+# X.Y # Final release
+#
+# Dev branch marker is: 'X.Y.dev' or 'X.Y.devN' where N is an integer.
+# 'X.Y.dev0' is the canonical version of 'X.Y.dev'
+#
+
+__version__ = '0.1.0'
diff --git a/steenroder/core.py b/steenroder/core.py
new file mode 100644
index 0000000..e1ff616
--- /dev/null
+++ b/steenroder/core.py
@@ -0,0 +1,912 @@
+import gudhi
+import psutil
+import time
+import numba as nb
+import numpy as np
+from functools import lru_cache
+from numba.cpython.unsafe.tuple import tuple_setitem
+from numba.np.unsafe.ndarray import to_fixed_tuple
+
+from .preprocessing import sort_filtration_by_dim
+
+
+list_of_int64_typ = nb.types.List(nb.int64)
+int64_2d_array_typ = nb.types.Array(nb.int64, 2, "C")
+
+# Determine the number of available physical cores
+N_PHYSICAL_CORES = psutil.cpu_count(logical=False)
+
+
+@nb.njit
+def _twist_reduction(coboundary, triangular, pivots_lookup):
+    """Core of the persistent relative cohomology reduction algorithm using the
+    clearing optimization."""
+    n = len(coboundary)
+
+    rel_idxs_to_clear = []
+    for j in range(n - 1, -1, -1):
+        highest_one = coboundary[j][0] if coboundary[j] else -1
+        pivot_col = pivots_lookup[highest_one]
+        while (highest_one != -1) and (pivot_col != -1):
+            coboundary[j] = _symm_diff(coboundary[j][1:],
+                                       coboundary[pivot_col][1:])
+            triangular[j] = _symm_diff(triangular[j],
+                                       triangular[pivot_col])
+            highest_one = coboundary[j][0] if coboundary[j] else -1
+            pivot_col = pivots_lookup[highest_one]
+        if highest_one != -1:
+            pivots_lookup[highest_one] = j
+            rel_idxs_to_clear.append(highest_one)
+
+    return np.asarray(rel_idxs_to_clear, dtype=np.int64)
+
+
+@lru_cache
+def _reduce_single_dim(dim):
+    len_tups_dim = dim + 1
+    tuple_typ_dim = nb.types.UniTuple(nb.int64, len_tups_dim)
+    len_tups_next_dim = dim + 2
+
+    @nb.njit
+    def _inner_reduce_single_dim(idxs_dim, tups_dim, rel_idxs_to_clear,
+                                 idxs_next_dim=None, tups_next_dim=None):
+        """R = MV"""
+        # 1) Construct sp2idx_dim as a dict simplex: relative (i.e.
+        # in-dimension) index
+        # 2) Initialize type of reduced_dim (needed for type inference)
+        # 3) Construct triangular_dim, with entries denoting relative (i.e.
+        # in-dimension) indices
+        spx2idx_dim = nb.typed.Dict.empty(tuple_typ_dim, nb.int64)
+        reduced_dim = nb.typed.List.empty_list(list_of_int64_typ)
+        triangular_dim = nb.typed.List.empty_list(list_of_int64_typ)
+        for i in range(len(idxs_dim)):
+            spx = to_fixed_tuple(tups_dim[i], len_tups_dim)
+            spx2idx_dim[spx] = i
+            reduced_dim.append([nb.int64(x) for x in range(0)])
+            triangular_dim.append([i])
+
+        # Populate reduced_dim as the coboundary matrix and apply clearing
+        # WARNING: Column entries denote relative (i.e. in-dimension) indices!
+        if idxs_next_dim is not None:
+            for j in range(len(idxs_next_dim)):
+                spx = to_fixed_tuple(tups_next_dim[j], len_tups_next_dim)
+                for face in _drop_elements(spx):
+                    reduced_dim[spx2idx_dim[face]].append(j)
+
+            for rel_idx in rel_idxs_to_clear:
+                reduced_dim[rel_idx] = [nb.int64(x) for x in range(0)]
+
+            pivots_lookup = np.full(len(idxs_next_dim), -1, dtype=np.int64)
+
+            rel_idxs_to_clear = _twist_reduction(reduced_dim, triangular_dim,
+                                                 pivots_lookup)
+        else:
+            pivots_lookup = np.empty(0, dtype=np.int64)
+
+        return (spx2idx_dim, reduced_dim, triangular_dim,
+                rel_idxs_to_clear, pivots_lookup)
+
+    return _inner_reduce_single_dim
+
+
+@nb.njit
+def _fix_triangular_after_clearing(triangular, reduced_prev_dim,
+                                   rel_idxs_to_clear, pivots_lookup_prev_dim):
+    """Massage the V matrix to maintain the R = DV decomposition after clearing,
+    as described in https://arxiv.org/abs/1908.02518, Sec. 3.2."""
+    for rel_idx in rel_idxs_to_clear:
+        triangular[rel_idx] = reduced_prev_dim[pivots_lookup_prev_dim[rel_idx]]
+
+
+def compute_reduced_triangular(filtration_by_dim):
+    """Compute the matrices R and V in the decomposition R = DV.
+
+    Explicitly, D is the coboundary matrix of the filtration, R is a reduced
+    matrix, and V is a full-rank upper-triangular one.
+
+    Parameters
+    ----------
+    filtration_by_dim : list of list of ndarray
+        For each dimension ``d``, a list of 2 aligned int arrays: the first is
+        a 1D array containing the (ordered) positional indices of all
+        ``d``-dimensional simplices in `filtration`; the second is a 2D array
+        whose ``i``-th row is the (sorted) collection of vertices defining the
+        ``i``-th ``d``-dimensional simplex.
+
+    Returns
+    -------
+    spx2idx : tuple of ``numba.typed.Dict``
+        One dictionary per simplex dimension. The dimension-``d`` dictionary has
+        the filtration ``d``-simplices (tuples of ints) as keys; the
+        corresponding values are the positional indices of those simplices
+        relative to the ``d``-dimensional portion of the filtration.
+
+    idxs : tuple of ndarray
+        For each dimension ``d``, this is the same as
+        ``filtration_by_dim[d][0]`` and is returned for convenience.
+
+    reduced : tuple of ``numba.typed.List``
+        One list of int per simplex dimension. ``reduced[d]`` is the
+        ``d``-dimensional part of the "R" matrix in R = DV. In the computation,
+        ``reduced[d][i]`` is initialized as the coboundary of the ``i``th input
+        simplex in ``filtration_by_dim[d][1]``, i.e. as the sorted list of
+        positional indices (relative to the ``(d+1)``-dimensional portion of the
+        filtration) of that simplex's cofacets.
+
+    triangular : tuple of ``numba.typed.List``
+        On list of int per simplex dimension. ``triangular[d]`` is the
+        ``d``-dimensional part of the "V" matrix in R = DV. In the computation,
+        ``triangular[d][i]`` is initialized as the singleton list ``[i]``.
+
+    """
+    maxdim = len(filtration_by_dim) - 1
+    spx2idx_idxs_reduced_triangular = []
+    # Initialize relative (i.e. in-dimension) indices to clear, as an empty
+    # int array in dim 0
+    rel_idxs_to_clear = np.empty(0, dtype=np.int64)
+    reduced_prev_dim = nb.typed.List.empty_list(list_of_int64_typ)
+    pivots_lookup_prev_dim = np.empty(0, dtype=np.int64)
+    for dim in range(maxdim):
+        reduction_dim = _reduce_single_dim(dim)
+        idxs_dim, tups_dim = filtration_by_dim[dim]
+        idxs_next_dim, tups_next_dim = filtration_by_dim[dim + 1]
+        (spx2idx_dim, reduced_dim, triangular_dim,
+         rel_idxs_to_clear_next_dim, pivots_lookup) = \
+            reduction_dim(idxs_dim,
+                          tups_dim,
+                          rel_idxs_to_clear,
+                          idxs_next_dim=idxs_next_dim,
+                          tups_next_dim=tups_next_dim)
+        _fix_triangular_after_clearing(triangular_dim,
+                                       reduced_prev_dim,
+                                       rel_idxs_to_clear,
+                                       pivots_lookup_prev_dim)
+        spx2idx_idxs_reduced_triangular.append((spx2idx_dim,
+                                                idxs_dim,
+                                                reduced_dim,
+                                                triangular_dim))
+        rel_idxs_to_clear = rel_idxs_to_clear_next_dim
+        reduced_prev_dim = reduced_dim
+        pivots_lookup_prev_dim = pivots_lookup
+
+    reduction_dim = _reduce_single_dim(maxdim)
+    idxs_dim, tups_dim = filtration_by_dim[maxdim]
+    spx2idx_dim, reduced_dim, triangular_dim, _, _ = \
+        reduction_dim(idxs_dim, tups_dim, rel_idxs_to_clear)
+    _fix_triangular_after_clearing(triangular_dim,
+                                   reduced_prev_dim,
+                                   rel_idxs_to_clear,
+                                   pivots_lookup_prev_dim)
+    spx2idx_idxs_reduced_triangular.append((spx2idx_dim,
+                                            idxs_dim,
+                                            reduced_dim,
+                                            triangular_dim))
+
+    return tuple(zip(*spx2idx_idxs_reduced_triangular))
+
+
+@nb.njit
+def compute_barcode_and_coho_reps(idxs, reduced, triangular,
+                              filtration_values=None):
+    """Compute the barcode of relative cohomology and cocycle representatives.
+
+    Parameters
+    ----------
+    idxs : tuple of ndarray
+        For each dimension ``d``, a 1D int array containing the (ordered)
+        positional indices of all ``d``-dimensional simplices in the filtration.
+
+    reduced : tuple of ``numba.typed.List``
+        One list of int per simplex dimension, representing the
+        ``d``-dimensional part of the "R" matrix in R = DV. In the same format
+        as returned by `compute_reduced_triangular`.
+
+    triangular : tuple of ``numba.typed.List``
+        One list of int per simplex dimension, representing the
+        ``d``-dimensional part of the "V" matrix in R = DV. In the same format
+        as returned by `compute_reduced_triangular`.
+
+    filtration_values : ndarray or None, optional, default: None
+        Optionally, a single 1D array of filtration values for each simplex in
+        the filtration and in all dimensions contained in `idxs`, `reduced` and
+        `triangular`. Bars with equal birth and death filtration values are
+        discarded.
+
+    Returns
+    -------
+    barcode : list of ndarray
+        For each dimension ``d``, a 2D int array of shape ``(n_bars, 2)``
+        containing the birth (entry 1) and death (entry 0) indices of persistent
+        relative homology classes in degree ``d``. See `filtration_values`.
+        Essential bars are represented by pairs with death equal to ``-1``. Bars
+        are sorted in order of decreasing birth indices.
+
+    coho_reps : list of ``numba.typed.List``
+        For each dimension ``d``, a list of representatives of persistent
+        relative cohomology classes in degree ``d``. Each such representative is
+        represented as a list of positional indices relative to the
+        ``d``-dimensional portion of the filtration. ``coho_reps[d][j]``
+        corresponds to ``barcode[d][j]``.
+
+    """
+    barcode = []
+    coho_reps = []
+
+    if filtration_values is None:
+        pairs_0 = []
+        coho_reps_0 = []
+        for i in range(len(idxs[0])):
+            if not reduced[0][i]:
+                pairs_0.append([-1, idxs[0][i]])
+                coho_reps_0.append(triangular[0][i])
+        pairs_0 = np.asarray(pairs_0)
+        lexsrt = _lexsort_barcode(pairs_0)
+        barcode.append(pairs_0[lexsrt])
+        coho_reps.append(nb.typed.List([coho_reps_0[k] for k in lexsrt]))
+
+        for dim in range(1, len(idxs)):
+            all_birth_indices = set()
+            pairs_dim = []
+            coho_reps_dim = []
+            for i in range(len(idxs[dim - 1])):
+                if reduced[dim - 1][i]:
+                    b = idxs[dim][reduced[dim - 1][i][0]]
+                    d = idxs[dim - 1][i]
+                    pairs_dim.append([d, b])
+                    coho_reps_dim.append(reduced[dim - 1][i])
+                    all_birth_indices.add(b)
+
+            for i in range(len(idxs[dim])):
+                if idxs[dim][i] not in all_birth_indices:
+                    if not reduced[dim][i]:
+                        pairs_dim.append([-1, idxs[dim][i]])
+                        coho_reps_dim.append(triangular[dim][i])
+
+            if not len(pairs_dim):
+                pairs_dim = np.empty((0, 2), dtype=np.int64)
+            else:
+                pairs_dim = np.asarray(pairs_dim)
+            lexsrt = _lexsort_barcode(pairs_dim)
+            barcode.append(pairs_dim[lexsrt])
+            coho_reps.append(nb.typed.List([coho_reps_dim[k] for k in lexsrt]))
+
+    else:
+        pairs_0 = []
+        coho_reps_0 = []
+        for i in range(len(idxs[0])):
+            if not reduced[0][i]:
+                pairs_0.append([-1, idxs[0][i]])
+                coho_reps_0.append(triangular[0][i])
+        pairs_0 = np.asarray(pairs_0)
+        lexsrt = _lexsort_barcode(pairs_0)
+        barcode.append(pairs_0[lexsrt])
+        coho_reps.append(nb.typed.List([coho_reps_0[k] for k in lexsrt]))
+
+        for dim in range(1, len(idxs)):
+            all_birth_indices = set()
+            pairs_dim = []
+            coho_reps_dim = []
+            for i in range(len(idxs[dim - 1])):
+                if reduced[dim - 1][i]:
+                    b = idxs[dim][reduced[dim - 1][i][0]]
+                    d = idxs[dim - 1][i]
+                    if filtration_values[b] != filtration_values[d]:
+                        pairs_dim.append([d, b])
+                        coho_reps_dim.append(reduced[dim - 1][i])
+                    all_birth_indices.add(b)
+
+            for i in range(len(idxs[dim])):
+                if idxs[dim][i] not in all_birth_indices:
+                    if not reduced[dim][i]:
+                        pairs_dim.append([-1, idxs[dim][i]])
+                        coho_reps_dim.append(triangular[dim][i])
+
+            if not len(pairs_dim):
+                pairs_dim = np.empty((0, 2), dtype=np.int64)
+            else:
+                pairs_dim = np.asarray(pairs_dim)
+            lexsrt = _lexsort_barcode(pairs_dim)
+            barcode.append(pairs_dim[lexsrt])
+            coho_reps.append(nb.typed.List([coho_reps_dim[k] for k in lexsrt]))
+
+    return barcode, coho_reps
+
+
+@nb.njit
+def _initialize_steenrod_matrix(num_dimensions):
+    return [nb.typed.List.empty_list(list_of_int64_typ)
+            for _ in range(num_dimensions)]
+
+
+@lru_cache
+def _populate_steenrod_matrix_single_dim(dim_plus_k):
+    length = dim_plus_k + 1
+
+    @nb.njit(parallel=True)
+    def _inner(coho_reps_dim, tups_dim, spx2idx_dim_plus_k, n_jobs=-1):
+        steenrod_matrix_dim_plus_k = \
+            nb.typed.List([[nb.int64(0) for _ in range(0)]
+                           for _ in coho_reps_dim])
+
+        if n_jobs == -1:
+            n_jobs = N_PHYSICAL_CORES
+
+        for job_idx in nb.prange(n_jobs):
+            for coho_reps_dim_idx in range(job_idx, len(coho_reps_dim), n_jobs):
+                rep = coho_reps_dim[coho_reps_dim_idx]
+                cocycle = tups_dim[np.asarray(rep)]
+
+                # STSQ
+                cochain = set(
+                    [to_fixed_tuple(np.empty(length, dtype=np.int64), length)
+                     for _ in range(0)]
+                )
+                for i in range(len(cocycle)):
+                    for j in range(i + 1, len(cocycle)):
+                        a, b = set(cocycle[i]), set(cocycle[j])
+                        u = a.union(b)
+                        if len(u) == length:
+                            u_tuple = to_fixed_tuple(np.asarray(sorted(u)),
+                                                     length)
+                            if u_tuple in spx2idx_dim_plus_k:
+                                a_bar, b_bar = a.difference(b), b.difference(a)
+                                u_bar = sorted(a_bar.union(b_bar))
+                                index = {}
+                                for v in a_bar.union(b_bar):
+                                    pos = u_tuple.index(v)
+                                    pos_bar = u_bar.index(v)
+                                    index[v] = (pos + pos_bar) % 2
+                                index_a = set()
+                                index_b = set()
+                                for v in a_bar:
+                                    index_a.add(index[v])
+                                for w in b_bar:
+                                    index_b.add(index[w])
+                                if (index_a == set([0])
+                                    and index_b == set([1])) \
+                                        or (index_a == set([1])
+                                            and index_b == set([0])):
+                                    cochain ^= {u_tuple}
+
+                steenrod_matrix_dim_plus_k[coho_reps_dim_idx] = \
+                    sorted([spx2idx_dim_plus_k[spx] for spx in cochain])
+
+        return steenrod_matrix_dim_plus_k
+
+    return _inner
+
+
+def compute_steenrod_matrix(k, coho_reps, filtration_by_dim, spx2idx, n_jobs=-1):
+    """Compute the k-th Steenrod matrices in each dimension.
+
+    Parameters
+    ----------
+    k : int
+        Positive integer defining the cohomology operation Sq^k to be performed.
+
+    coho_reps : list of ``numba.typed.List``
+        For each dimension ``d``, a list of representatives of persistent
+        relative cohomology classes in degree ``d``. In the same format as
+        returned by `compute_barcode_and_coho_reps`.
+
+    filtration_by_dim : list of list of ndarray
+        For each dimension ``d``, a list of 2 aligned int arrays: the first is
+        a 1D array containing the (ordered) positional indices of all
+        ``d``-dimensional simplices in `filtration`; the second is a 2D array
+        whose ``i``-th row is the (sorted) collection of vertices defining the
+        ``i``-th ``d``-dimensional simplex.
+
+    spx2idx : tuple of ``numba.typed.Dict``
+        One dictionary per simplex dimension. The dimension-``d`` dictionary has
+        the filtration ``d``-simplices (tuples of ints) as keys; the
+        corresponding values are the positional indices of those simplices
+        relative to the ``d``-dimensional portion of the filtration.
+
+    n_jobs : int, optional, default: ``-1``
+        [Experimental] Controls the number of threads to be used during parallel
+        computation of the Steenrod squares. ``-1`` means using all available
+        physical cores.
+
+    Returns
+    -------
+    steenrod_matrix : list of ``numba.typed.List``
+        One list per simplex dimension. ``steenrod_matrix[d][j]`` is the result
+        of computing the Steenrod square of ``coho_reps[d - k][j]``.
+
+    """
+    steenrod_matrix = _initialize_steenrod_matrix(k)
+
+    for dim, coho_reps_dim in enumerate(coho_reps[:-k]):
+        dim_plus_k = dim + k
+        tups_dim = filtration_by_dim[dim][1]
+        spx2idx_dim_plus_k = spx2idx[dim + k]
+        populate_steenrod_matrix_single_dim = \
+            _populate_steenrod_matrix_single_dim(dim_plus_k)
+        steenrod_matrix_dim_plus_k = populate_steenrod_matrix_single_dim(
+            coho_reps_dim, tups_dim, spx2idx_dim_plus_k, n_jobs=n_jobs
+        )
+        steenrod_matrix.append(steenrod_matrix_dim_plus_k)
+
+    return steenrod_matrix
+
+
+@nb.njit
+def _steenrod_barcode_single_dim(steenrod_matrix_dim, n_idxs_dim, idxs_prev_dim,
+                                 reduced_prev_dim, births_dim_minus_k):
+    # Construct augmented matrix
+    augmented = []
+    for i in range(len(reduced_prev_dim)):
+        augmented.append([nb.int64(x) for x in reduced_prev_dim[i]])
+    for i in range(len(steenrod_matrix_dim)):
+        augmented.append([nb.int64(x) for x in steenrod_matrix_dim[i]])
+
+    n_cols_red = len(reduced_prev_dim)
+    n_cols_st = len(steenrod_matrix_dim)
+    n_births_dim_minus_k = len(births_dim_minus_k)
+
+    pivots_lookup = np.full(n_idxs_dim, -1, dtype=np.int64)
+    alive = np.ones(n_births_dim_minus_k, dtype=np.bool_)
+    steenrod_barcode_dim = []
+
+    for n_cols_st_curr, birth in enumerate(births_dim_minus_k):
+        if idxs_prev_dim[-1] >= birth:
+            break
+    for i in np.arange(n_cols_red)[::-1]:
+        idx = idxs_prev_dim[i]
+        next_birth = births_dim_minus_k[n_cols_st_curr]
+        is_idx_geq = idx >= next_birth
+        if is_idx_geq:
+            if i:
+                next_idx = idxs_prev_dim[i - 1]
+                for jj in range(n_cols_st_curr, n_cols_st):
+                    if next_idx >= births_dim_minus_k[jj]:
+                        break
+                    n_cols_st_curr += 1
+            else:
+                n_cols_st_curr = n_cols_st
+
+        if augmented[i]:
+            pivots_lookup[augmented[i][0]] = i
+
+        pivot_col_idxs_from_st = []
+        for ii in range(n_cols_red, n_cols_red + n_cols_st_curr):
+            highest_one = augmented[ii][0] if augmented[ii] else -1
+            pivot_col = pivots_lookup[highest_one]
+            while (highest_one != -1) and (pivot_col != -1):
+                augmented[ii] = _symm_diff(augmented[ii][1:],
+                                           augmented[pivot_col][1:])
+                highest_one = augmented[ii][0] if augmented[ii] else -1
+                pivot_col = pivots_lookup[highest_one]
+            if highest_one != -1:
+                pivots_lookup[highest_one] = ii
+                # Record pivot indices coming from Steenrod part of augmented
+                pivot_col_idxs_from_st.append(highest_one)
+            elif alive[ii - n_cols_red]:
+                alive[ii - n_cols_red] = False
+                birth = births_dim_minus_k[ii - n_cols_red]
+                if idx < birth:
+                    steenrod_barcode_dim.append([idx, birth])
+
+        # Reset pivots_lookup for next iteration
+        for col_idx in pivot_col_idxs_from_st:
+            pivots_lookup[col_idx] = -1
+
+    for i in range(n_births_dim_minus_k):
+        if alive[i]:
+            steenrod_barcode_dim.append([-1, births_dim_minus_k[i]])
+
+    return steenrod_barcode_dim
+
+
+def compute_steenrod_barcode(k, steenrod_matrix, idxs, reduced, barcode,
+                         filtration_values=None):
+    """Compute the Sq^k-barcode of relative cohomology.
+
+    Parameters
+    ----------
+    k : int
+        Positive integer defining the cohomology operation Sq^k to be performed.
+
+    steenrod_matrix : list of ``numba.typed.List``
+        One list per simplex dimension. ``steenrod_matrix[d][j]`` is the result
+        of computing the Steenrod square of the ``j``th latest (by birth)
+        persistent relative cohomology representative in degree ``d`` (and this
+        representative must represent bar ``barcode[d][j]``). See
+        `compute_steenrod_matrix`.
+
+    idxs : tuple of ndarray
+        For each dimension ``d``, a 1D int array containing the (ordered)
+        positional indices of all ``d``-dimensional simplices in the filtration.
+
+    reduced : tuple of ``numba.typed.List``
+        One list of int per simplex dimension, representing the
+        ``d``-dimensional part of the "R" matrix in R = DV. In the same format
+        as returned by `compute_reduced_triangular`.
+
+    barcode : list of ndarray
+        For each dimension ``d``, a 2D int array of shape ``(n_bars, 2)``
+        containing the birth (entry 1) and death (entry 0) indices of persistent
+        relative homology classes in degree ``d``. Essential bars must be
+        represented by pairs with death equal to ``-1``. Bars must be sorted in
+        order of decreasing birth indices.
+
+    filtration_values : ndarray or None, optional, default: None
+        Optionally, a single 1D array of filtration values for each simplex in
+        the filtration and in all dimensions contained in `idxs` and `reduced`.
+        Steenrod bars with equal birth and death filtration values are
+        discarded.
+
+    Returns
+    -------
+    steenrod_barcode : list of ndarray
+        The (relative) Sq^k-barcode. For each dimension ``d``, a 2D int array
+        of shape ``(n_bars, 2)`` containing the birth (entry 1) and death (entry
+        0) indices of Steenrod bars. Essential Steenrod bars are represented by
+        pairs with death equal to ``-1``.
+
+    """
+    def nontrivial_bars(barcode_dim):
+        infinite_bars = barcode_dim[:, 0] == -1
+        return np.logical_or(
+            infinite_bars,
+            np.logical_and(np.logical_not(infinite_bars),
+                           (filtration_values[barcode_dim[:, 0]] !=
+                            filtration_values[barcode_dim[:, 1]]))
+        )
+
+    steenrod_barcode = [np.empty((0, 2), dtype=np.int64) for _ in range(k)]
+    for dim in range(k, len(steenrod_matrix)):
+        births_dim_minus_k = barcode[dim - k][:, 1]
+        idxs_dim = idxs[dim]
+        idxs_prev_dim = idxs[dim - 1]
+        reduced_prev_dim = reduced[dim - 1]
+        steenrod_barcode_dim = _steenrod_barcode_single_dim(steenrod_matrix[dim],
+                                                      len(idxs_dim),
+                                                      idxs_prev_dim,
+                                                      reduced_prev_dim,
+                                                      births_dim_minus_k)
+        # NB: Conversion to array must happen outside jitted code due to
+        # https://github.com/numba/numba/issues/3579
+        steenrod_barcode_dim = \
+            np.asarray(steenrod_barcode_dim, dtype=np.int64).reshape((-1, 2))
+        if filtration_values is not None:
+            nontrivial_mask = nontrivial_bars(steenrod_barcode_dim)
+            steenrod_barcode_dim = steenrod_barcode_dim[nontrivial_mask]
+        steenrod_barcode.append(steenrod_barcode_dim)
+
+    return steenrod_barcode
+
+
+def barcodes(
+        k, filtration, absolute=False, filtration_values=None,
+        return_filtration_values=False, maxdim=None, n_jobs=1, verbose=False
+):
+    """Given a filtration, compute ordinary persistent (relative or absolute)
+    cohomology barcodes and relative Steenrod barcodes.
+
+    Parameters
+    ----------
+    k : int
+        Positive integer defining the cohomology operation Sq^k to be performed.
+
+    filtration : sequence of list-like of int
+        Represents a simplex-wise filtration. Entry ``i`` is a list/tuple/set
+        containing the integer indices of the vertices defining the ``i``-th
+        simplex in the filtration.
+
+    absolute : bool, optional, default: ``False``
+        If ``True``, return the ordinary persistent absolute homology barcode,
+        and move inessential relative Steenrod bars to one degree lower while
+        keeping essential bars in their degree. If ``False``, return the
+        ordinary persistent relative cohomology barcode and relative Steenrod
+        barcode.
+
+    filtration_values : ndarray or None, optional, default: None
+        Optionally, a single 1D array of filtration values for each simplex in
+        the filtration. Ordinary and Steenrod bars with equal birth and death
+        filtration values are discarded by the computation.
+
+    return_filtration_values : bool, optional, default: ``False``
+        If ``True``, birth and deaths will be expressed as filtration values
+        instead of filtration indices. Ignored if `filtration_values` is
+        ``None``.
+
+    maxdim : int or None, optional, default: None
+        Maximum simplex dimension to be included. ``None`` means that all
+        simplices are included.
+
+    n_jobs : int, optional, default: ``1``
+        [Experimental] Controls the number of threads to be used during parallel
+        computation of the Steenrod squares. ``-1`` means using all available
+        physical cores.
+
+    verbose : bool, optional, default: ``False``
+        Whether to print timings for the intermediate steps in the computation.
+
+    Returns
+    -------
+    barcode : list of ndarray
+        For each dimension ``d``, a 2D int or float array of shape ``(n_bars,
+        2)`` containing the births and deaths of persistent relative cohomology
+        classes in degree ``d``. If `absolute` is ``False``, the birth of a bar
+        is in entry 1 and the death in entry 0; otherwise, the positions are
+        reversed. Births and death are expressed either as global filtration
+        indices or as filtration values depending on `filtration_values` and
+        `return_filtration_values`. If they are expressed as indices, essential
+        bars have death equal to ``-1``; otherwise, essential bars have death
+        equal to ``numpy.inf``.
+
+    steenrod_barcode : list of ndarray
+        The (relative) Sq^k-barcode. For each dimension ``d``, a 2D int or float
+        array of shape ``(n_bars, 2)`` containing the birth (entry 1) and death
+        (entry 0) indices of Steenrod bars. The same conventions as for
+        `barcode` are used for birth and death values.
+
+    """
+    if verbose:
+        tic = time.time()
+    filtration_by_dim = sort_filtration_by_dim(filtration, maxdim=maxdim)
+    spx2idx, idxs, reduced, triangular = \
+        compute_reduced_triangular(filtration_by_dim)
+    barcode, coho_reps = \
+        compute_barcode_and_coho_reps(idxs, reduced, triangular,
+                                  filtration_values=filtration_values)
+
+    if k == 0:  # Sq^0 is the identity
+        return barcode, barcode
+
+    if k < 0:  # Sq^k is the zero map
+        return barcode, list()
+
+    if verbose:
+        toc = time.time()
+        print(f"Usual barcode computed, time taken: {(toc - tic):.3f} s")
+        tic = time.time()
+    steenrod_matrix = compute_steenrod_matrix(k, coho_reps, filtration_by_dim,
+                                          spx2idx, n_jobs=n_jobs)
+    if verbose:
+        toc = time.time()
+        print(f"Steenrod matrix computed, time taken: {(toc - tic):.3f} s")
+        tic = time.time()
+    steenrod_barcode = compute_steenrod_barcode(k, steenrod_matrix, idxs, reduced,
+                                      barcode,
+                                      filtration_values=filtration_values)
+    if verbose:
+        toc = time.time()
+        print(f"Steenrod barcode computed, time taken: {(toc - tic):.3f} s")
+
+    if absolute:
+        barcode = _to_absolute_barcode(
+            barcode, filtration_values=filtration_values,
+            return_filtration_values=return_filtration_values
+        )
+        steenrod_barcode = _to_absolute_barcode(
+            steenrod_barcode, filtration_values=filtration_values,
+            return_filtration_values=return_filtration_values
+        )
+
+        return barcode, steenrod_barcode
+
+    elif return_filtration_values and (filtration_values is not None):
+        barcode = _to_values_barcode(barcode, filtration_values)
+        steenrod_barcode = _to_values_barcode(steenrod_barcode, filtration_values)
+
+    return barcode, steenrod_barcode
+
+
+def _to_absolute_barcode(rel_barcode, filtration_values=None,
+                         return_filtration_values=True):
+    abs_barcode = []
+
+    if (not return_filtration_values) or (filtration_values is None):
+        dtype = np.int64
+        for dim, rel_barcode_dim in enumerate(rel_barcode):
+            abs_barcode_dim = []
+            for pair in rel_barcode_dim:
+                if pair[0] == -1:
+                    abs_barcode_dim.append((pair[1], -1))
+                else:
+                    abs_barcode[dim - 1].append((pair[0], pair[1]))
+            abs_barcode.append(abs_barcode_dim)
+    else:
+        dtype = filtration_values.dtype
+        for dim, rel_barcode_dim in enumerate(rel_barcode):
+            abs_barcode_dim = []
+            for pair in rel_barcode_dim:
+                if pair[0] == -1:
+                    abs_barcode_dim.append(
+                        (filtration_values[pair[1]], np.inf)
+                    )
+                else:
+                    abs_barcode[dim - 1].append(
+                        (filtration_values[pair[0]], filtration_values[pair[1]])
+                    )
+            abs_barcode.append(abs_barcode_dim)
+
+    return [np.array(abs_barcode_dim, dtype=dtype).reshape(-1, 2)
+            for abs_barcode_dim in abs_barcode]
+
+
+def _to_values_barcode(barcode, filtration_values):
+    values_barcode = []
+    for dim, barcode_dim in enumerate(barcode):
+        values_barcode_dim = []
+        for pair in barcode_dim:
+            if pair[0] == -1:
+                values_barcode_dim.append(
+                    (-np.inf, filtration_values[pair[1]])
+                )
+            else:
+                values_barcode_dim.append(
+                    (filtration_values[pair[0]], filtration_values[pair[1]])
+                )
+        values_barcode.append(
+            np.array(values_barcode_dim,
+                     dtype=filtration_values.dtype).reshape(-1, 2)
+        )
+
+    return values_barcode
+
+
+@nb.njit
+def _symm_diff(x, y):
+    n = len(x)
+    m = len(y)
+    result = []
+    i = 0
+    j = 0
+    while (i < n) and (j < m):
+        if x[i] < y[j]:
+            result.append(x[i])
+            i += 1
+        elif y[j] < x[i]:
+            result.append(y[j])
+            j += 1
+        else:
+            i += 1
+            j += 1
+
+    while i < n:
+        result.append(x[i])
+        i += 1
+
+    while j < m:
+        result.append(y[j])
+        j += 1
+
+    return result
+
+
+@nb.njit
+def _lexsort_barcode(arr):
+    return np.argsort(arr[:, 1])[::-1]
+
+
+@nb.njit
+def _drop_elements(tup: tuple):
+    for x in range(len(tup)):
+        empty = tup[:-1]  # Not empty, but the right size and will be mutated
+        idx = 0
+        for i in range(len(tup)):
+            if i != x:
+                empty = tuple_setitem(empty, idx, tup[i])
+                idx += 1
+        yield empty
+
+
+def rips_barcodes(X, max_edge_length=np.inf, distance_matrix=False,
+                  k=1, max_simplex_dimension=2, absolute=False, n_jobs=1,
+                  verbose=False):
+    """Construct a Rips filtration and compute ordinary persistent (relative or
+    absolute) cohomology barcodes and relative Steenrod barcodes.
+
+    Parameters
+    ----------
+    X : ndarray of shape (n_samples, n_features) or (n_samples, n_samples)
+        Input point cloud (if `distance_matrix` is ``False``) or distance matrix
+        (if `distance_matrix` is ``True``).
+
+    max_edge_length : float, optional, default: ``np.inf``
+        Vietoris–Rips filtration threshold.
+
+    distance_matrix : bool, optional, default: ``False``
+        Whether the input `X` is a distance matrix or a point cloud.
+
+    k : int
+        Positive integer defining the cohomology operation Sq^k to be performed.
+
+    max_simplex_dimension : int, optional, default: ``2``
+        Maximum simplex dimension to be included in the Rips filtration.
+
+    absolute : bool, optional, default: ``False``
+        If ``True``, return the ordinary persistent absolute homology barcode,
+        and move inessential relative Steenrod bars to one degree lower while
+        keeping essential bars in their degree. If ``False``, return the
+        ordinary persistent relative cohomology barcode and relative Steenrod
+        barcode.
+
+    n_jobs : int, optional, default: ``1``
+        [Experimental] Controls the number of threads to be used during parallel
+        computation of the Steenrod squares. ``-1`` means using all available
+        physical cores.
+
+    verbose : bool, optional, default: ``False``
+        Whether to print information about the intermediate steps in the
+        computation.
+
+    Returns
+    -------
+    barcode : list of ndarray
+        For each dimension ``d``, a 2D int or float array of shape ``(n_bars,
+        2)`` containing the births and deaths of persistent relative cohomology
+        classes in degree ``d``. If `absolute` is ``False``, the birth of a bar
+        is in entry 1 and the death in entry 0; otherwise, the positions are
+        reversed. Births and death are expressed either as global filtration
+        indices or as filtration values depending on `filtration_values` and
+        `return_filtration_values`. If they are expressed as indices, essential
+        bars have death equal to ``-1``; otherwise, essential bars have death
+        equal to ``numpy.inf``.
+
+    steenrod_barcode : list of ndarray
+        The (relative) Sq^k-barcode. For each dimension ``d``, a 2D int or float
+        array of shape ``(n_bars, 2)`` containing the birth (entry 1) and death
+        (entry 0) indices of Steenrod bars. The same conventions as for
+        `barcode` are used for birth and death values.
+
+    """
+    spx_tree_kwargs = {
+        "distance_matrix" if distance_matrix else "points": X,
+        "max_edge_length": max_edge_length
+    }
+    spx_tree = gudhi.RipsComplex(**spx_tree_kwargs).\
+        create_simplex_tree(max_dimension=1)
+
+    len_filtration = None
+    if verbose:
+        for i, _ in enumerate(spx_tree.get_filtration()):
+            pass
+        len_filtration = i + 1
+        print(f"Filtration with {len_filtration} simplices up to dimension 1 "
+              "initially")
+
+        print("Running Edge Collapse...")
+
+    it = 0
+    while True:
+        spx_tree.collapse_edges()
+
+        for i, _ in enumerate(spx_tree.get_filtration()):
+            pass
+
+        if i + 1 == len_filtration:
+            break
+        else:
+            len_filtration = i + 1
+            it += 1
+            if verbose:
+                print(f"  EC iteration {it}: {len_filtration} simplices up to "
+                      "dimension 1")
+
+    # Construct the Rips filtration and run Steenroder
+    spx_tree.expansion(max_simplex_dimension)
+    if verbose:
+        for i, _ in enumerate(spx_tree.get_filtration()):
+            pass
+        print(f"There are {i} simplices up to dimension "
+              f"{max_simplex_dimension} after EC.")
+
+    filtration, filtration_values = \
+        zip(*((tuple(t[0]), t[1]) for t in spx_tree.get_filtration()))
+    filtration_values = np.asarray(filtration_values, dtype=np.float32)
+
+    barcode, steenrod_barcode = barcodes(
+        k,
+        filtration,
+        filtration_values=filtration_values,
+        absolute=absolute,
+        return_filtration_values=True,
+        n_jobs=n_jobs,
+        verbose=verbose
+    )
+
+    return barcode, steenrod_barcode
diff --git a/steenroder/plotting.py b/steenroder/plotting.py
index 9ee3a23..a64f703 100644
--- a/steenroder/plotting.py
+++ b/steenroder/plotting.py
@@ -30,7 +30,7 @@ def _compute_bounds(barcode):
     return min_val_display, max_val_display, posinfinity_val, neginfinity_val
 
 
-def plot_diagrams(barcode, st_barcode, k=None, kind=None, tex=False,
+def plot_diagrams(barcode, steenrod_barcode, k=None, kind=None, tex=False,
                   plotly_params=None):
     """Plot a regular persistence barcode and a Steenrod barcode as diagrams on
     a common birth-death plane.
@@ -42,14 +42,14 @@ def plot_diagrams(barcode, st_barcode, k=None, kind=None, tex=False,
         with 2 columns, containing birth-death pairs in degree ``d``, to be used
         as coordinates in the two-dimensional plot.
 
-    st_barcode : list of ndarray
+    steenrod_barcode : list of ndarray
         The (relative) Sq^k-barcode to plot. For each dimension ``d``, a 2D
         array with 2 columns, containing the birth-death pairs of Steenrod bars
         in degree ``d``, to be used as coordinates in the two-dimensional plot.
 
     k : int or None, optional, default: ``None``
         Positive integer defining the cohomology operation Sq^k that was
-        performed to obtain `st_barcode`. Only used for labelling.
+        performed to obtain `steenrod_barcode`. Only used for labelling.
 
     kind : ``"R"`` | ``"A"`` or None, optional, default: ``None``
         Whether the barcodes to be plotted come from absolute or relative
@@ -96,10 +96,10 @@ def _connect_st_label(st_label, h_label, tex):
     else:
         legend_title = "Cohomology"
         h_subscript = r"" if tex else ""
-    homology_dimensions = range(max(len(barcode), len(st_barcode)))
+    homology_dimensions = range(max(len(barcode), len(steenrod_barcode)))
 
     min_val_display, max_val_display, posinfinity_val, neginfinity_val = \
-        _compute_bounds(barcode + st_barcode)
+        _compute_bounds(barcode + steenrod_barcode)
 
     fig = gobj.Figure()
     fig.add_trace(gobj.Scatter(
@@ -115,7 +115,7 @@ def _connect_st_label(st_label, h_label, tex):
         h_label = (rf"\mathcal{{H}}^{{{dim}}}" if tex else f"H^{dim}") + \
                   h_subscript
         h_label, st_h_label = _connect_st_label(st_label, h_label, tex)
-        for label, symbol, ms, bc in ([st_h_label, "diamond", 10, st_barcode],
+        for label, symbol, ms, bc in ([st_h_label, "diamond", 10, steenrod_barcode],
                                       [h_label, "circle", 8, barcode]):
             subbc = bc[dim].copy()
             unique, inverse, counts = np.unique(
diff --git a/steenroder/preprocessing.py b/steenroder/preprocessing.py
new file mode 100644
index 0000000..4b5c6b3
--- /dev/null
+++ b/steenroder/preprocessing.py
@@ -0,0 +1,181 @@
+import numpy as np
+from itertools import combinations, chain
+from random import shuffle
+from gudhi import SimplexTree
+
+
+def make_simplices(top_faces):
+    """Return all simplices of a complex given its top ones.
+
+    Parameters
+    ----------
+    top_faces : list of tuple of int
+        Represents the top faces of a simplicial complex.
+
+    Returns
+    -------
+    simplices : list of tuple of int
+        Represent all simplices in the minimal complex with the given top
+        dimensional simplices.
+
+    """
+    dim = len(top_faces[0]) - 1
+    simplices = set()
+    for top_face in top_faces:
+        for d in range(dim + 1):
+            simplices |= set(tuple(sorted(s)) for s in combinations(top_face, d + 1))
+    return list(simplices)
+
+
+def make_filtration(simplices):
+    """Return a filtration of a given complex.
+
+    Parameters
+    ----------
+    simplices : list of tuple of int
+        Represents the simplices of a simplicial complex.
+
+    Returns
+    -------
+    filtration : list of tuple of int
+        Represents a filtration for the given simplicial complex.
+
+    """
+    st = SimplexTree()
+    shuffle(simplices)
+    for spx in simplices:
+        st.insert(spx)
+    st.make_filtration_non_decreasing()
+    filtration = [tuple(pair[0]) for pair in st.get_filtration()]
+    return filtration
+
+
+def sort_filtration_by_dim(filtration, maxdim=None):
+    """Organize an input filtration by dimension.
+
+    Parameters
+    ----------
+    filtration : sequence of list-like of int
+        Represents a simplex-wise filtration. Entry ``i`` is a list/tuple/set
+        containing the integer indices of the vertices defining the ``i`` th
+        simplex in the filtration.
+
+    maxdim : int or None, optional, default: None
+        Maximum simplex dimension to be included. ``None`` means that all
+        simplices are included.
+
+    Returns
+    -------
+    filtration_by_dim : list of list of ndarray
+        For each dimension ``d``, a list of 2 aligned int arrays: the first is
+        a 1D array containing the (ordered) positional indices of all
+        ``d``-dimensional simplices in `filtration`; the second is a 2D array
+        whose ``i``-th row is the (sorted) collection of vertices defining the
+        ``i``-th ``d``-dimensional simplex.
+
+    """
+    if maxdim is None:
+        maxdim = max(map(len, filtration)) - 1
+
+    filtration_by_dim = [[] for _ in range(maxdim + 1)]
+    for i, spx in enumerate(filtration):
+        spx_tup = tuple(sorted(spx))
+        dim = len(spx_tup) - 1
+        if dim <= maxdim:
+            filtration_by_dim[dim].append([i, spx_tup])
+
+    for dim, filtr in enumerate(filtration_by_dim):
+        filtration_by_dim[dim] = [np.asarray(x, dtype=np.int64)
+                                  for x in zip(*filtr)]
+
+    return filtration_by_dim
+
+
+def cone(simplices):
+    """Return a filtration of the cone of a complex.
+
+    Parameters
+    ----------
+    simplices : list of tuple of int
+        Represents the simplices of a simplicial complex.
+
+    Returns
+    -------
+    cone : list of tuple of int
+        Represents a filtration of the cone of the given simplicial complex.
+
+    """
+    vertices = set(v for v in chain.from_iterable(simplices))
+    v = max(vertices) + 1
+    new_simplices = simplices + [tuple(list(s) + [v]) for s in simplices] + [(v,)]
+    cone = make_filtration(new_simplices)
+    return cone
+
+
+def suspension(simplices):
+    """Return a filtration of the suspension of a complex.
+
+    Parameters
+    ----------
+    simplices : list of tuple of int
+        Represents the simplices of a simplicial complex.
+
+    Returns
+    -------
+    suspension : list of tuple of int
+        Represents a filtration of the suspension of the given simplicial complex.
+
+    """
+    vertices = set(v for v in chain.from_iterable(simplices))
+    v, w = max(vertices) + 1, max(vertices) + 2
+    new_simplices = list(simplices)
+    new_simplices += [tuple(list(s) + [v]) for s in simplices] + [(v,)] +\
+                     [tuple(list(s) + [w]) for s in simplices] + [(w,)]
+    suspension = make_filtration(new_simplices)
+    return suspension
+
+
+def sphere(n):
+    """Return a filtration of the n-sphere.
+
+    Parameters
+    ----------
+    n : int
+        The dimension of the sphere.
+
+    Returns
+    -------
+    sphere : list of tuple of int
+        Represents a filtration of the n-sphere.
+
+    """
+    top_faces = [s for s in combinations(range(1, n + 3), n + 1)]
+    simplices = make_simplices(top_faces)
+    return make_filtration(simplices)
+
+
+def wedge(simplices_1, simplices_2):
+    """Return a filtration of the wedge of two complexes.
+
+    Parameters
+    ----------
+    simplices_1 : list of tuple of int
+        Represents the first simplicial complex.
+
+    simplices_2 : list of tuple of int
+        Represents the second simplicial complex.
+
+    Returns
+    -------
+    wedge : list of tuple of int
+        Represents a filtration of the wedge of the given simplicial complexes.
+
+    """
+    vertices_1 = set(v for v in chain.from_iterable(simplices_1))
+    v = max(vertices_1)
+    simplices = list(simplices_1)
+    simplices.remove((v,))
+    for s in simplices_2:
+        simplices.append(tuple(map(lambda i: i + v - 1, s)))
+    wedge = make_filtration(simplices)
+    return wedge


Test patch:
diff --git a/steenroder/tests/__init__.py b/steenroder/tests/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/steenroder/tests/test_core.py b/steenroder/tests/test_core.py
new file mode 100644
index 0000000..06e8cab
--- /dev/null
+++ b/steenroder/tests/test_core.py
@@ -0,0 +1,2 @@
+def test_core():
+    pass