Repository: ArnaoutLab/diversity (language: python)

Issue:
Power mean with 0 weights
@Elliot-D-Hill: Unit testing power mean revealed the following bugs:

1. Power mean of order 0 is defined as geometric mean which equals 1 when all weights are 0. However, diversity should be 0 in that case. How and where should we fix this? Some thoughts:
    - Technically power mean is correct, so somehow Metacommunity should take care of this case
    - I wouldn't expect anyone to have a good reason to pass all-0 weights because that is not interesting, so maybe it's a matter of input validation?

2. When order < -100 or > 100 and all weights are 0, power_mean resorts to "initial" keyword arguments of numpy.amin and numpy.amax, which are inf and -inf, respectively. Technically, the power mean should be undefined. I suggest validating and raising exception.

_Originally posted by @JasperBraun in https://github.com/Elliot-D-Hill/metacommunity/issues/7#issuecomment-1008064776_
Fix metacommunity.SimilarityFromFile
Include species_order derivation in constructor.


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 5ee3fe1c..46eb8843 100644
--- a/.gitignore
+++ b/.gitignore
@@ -133,3 +133,4 @@ test_files/
 
 junk_code.txt
 examples.ipynb
+CHANGELOG
diff --git a/README.md b/README.md
index e489262b..53a8a139 100644
--- a/README.md
+++ b/README.md
@@ -1,2 +1,673 @@
-# diversity-with-similarity
-Alpha and beta diversity with species similarity
+# diversity: similarity-sensitive diversity indices
+
+![Tests](https://github.com/Elliot-D-Hill/diversity/actions/workflows/tests.yml/badge.svg)
+[![Python 3.8](https://img.shields.io/badge/python-3.8-blue.svg)](https://www.python.org/downloads/release/python-380/)
+
+* [About](#about)
+* [Installation](#installation)
+* [Usage and Examples](#usage-and-examples)
+    + [Python](#python)
+        + [Similarity from array or dataframe](#similarity-from-array-or-dataframe)
+        + [Similarity from file](#similarity-from-file)
+        + [Similarity from function](#similarity-from-function)
+    + [Command line interface](#command-line-interface)
+        + [Metagenomics dataset](#metagenomics-dataset)
+* [Background](#background)
+    + [Diversity indices](#diversity-indices)
+    + [Partitioned diversity](#partitioned-diversity)
+    + [Frequency-sensitive diversity](#frequency-sensitive-diversity)
+    + [Similarity-sensitive diversity](#similarity-sensitive-diversity)
+    + [One package to calculate them all](#one-package-to-calculate-them-all)
+* [Alternatives](#alternatives)
+
+## About
+
+For a rigorous mathematical treatment of the diversity indices calculated by `diversity`, see [Reeve et al., 2014](https://arxiv.org/abs/1404.6520). A brief informal discussion can be found in the [background section](#background).
+
+## Installation
+
+To install, execute:
+```bash
+pip install diversity
+```
+
+**Dependencies:**
+- pytest >= 6.2.5
+- numpy >= 1.22.0
+- pandas >= 1.3.5
+
+## Usage and Examples
+
+
+### Python
+
+```python
+from diversity.metacommunity import make_metacommunity
+from pandas import read_csv
+from numpy.linalg import norm
+from scipy.spatial.distance import pdist, squareform
+```
+The data for the following examples can be found at this [Gist](https://gist.github.com/Elliot-D-Hill/1200c74c83a4be1f63f9ef56b333bb99).
+
+```python
+gist_url = 'https://gist.github.com/Elliot-D-Hill/1200c74c83a4be1f63f9ef56b333bb99/raw/71e2c2f026b9a495023e6e71ee1dcc545cb53972/'
+```
+
+
+#### Similarity from array or dataframe
+
+First we load the [Palmer penguins dataset](https://github.com/allisonhorst/palmerpenguins) (see [Gorman, Williams, and Fraser, 2014](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0090081) for more details). To simplify our examples, we have preprocessed the dataset to include only the means of four features for the three penguin species (Adelie, Chinstrap, and Gentoo).
+
+
+```python
+penguin_features_filepath = gist_url + 'penguin_features.csv'
+penguin_features = read_csv(penguin_features_filepath)
+penguin_features.head()
+```
+
+
+
+
+<div>
+<table border="1" class="dataframe">
+  <thead>
+    <tr style="text-align: right;">
+      <th></th>
+      <th>species</th>
+      <th>culmen_length_mm</th>
+      <th>culmen_depth_mm</th>
+      <th>flipper_length_mm</th>
+      <th>body_mass_g</th>
+    </tr>
+  </thead>
+  <tbody>
+    <tr>
+      <th>0</th>
+      <td>Adelie</td>
+      <td>38.823973</td>
+      <td>18.347260</td>
+      <td>190.102740</td>
+      <td>3706.164384</td>
+    </tr>
+    <tr>
+      <th>1</th>
+      <td>Chinstrap</td>
+      <td>48.833824</td>
+      <td>18.420588</td>
+      <td>195.823529</td>
+      <td>3733.088235</td>
+    </tr>
+    <tr>
+      <th>2</th>
+      <td>Gentoo</td>
+      <td>47.542500</td>
+      <td>15.002500</td>
+      <td>217.233333</td>
+      <td>5090.625000</td>
+    </tr>
+  </tbody>
+</table>
+</div>
+
+
+
+Next, we load in the penguin subcommunity-species counts data.
+
+
+```python
+penguin_counts_filepath = gist_url + 'penguin_counts.csv'
+penguin_counts = read_csv(penguin_counts_filepath)
+penguin_counts.head()
+```
+
+
+
+
+<div>
+
+<table border="1" class="dataframe">
+  <thead>
+    <tr style="text-align: right;">
+      <th></th>
+      <th>subcommunity</th>
+      <th>species</th>
+      <th>count</th>
+    </tr>
+  </thead>
+  <tbody>
+    <tr>
+      <th>0</th>
+      <td>Biscoe</td>
+      <td>Gentoo</td>
+      <td>120</td>
+    </tr>
+    <tr>
+      <th>1</th>
+      <td>Dream</td>
+      <td>Chinstrap</td>
+      <td>68</td>
+    </tr>
+    <tr>
+      <th>2</th>
+      <td>Dream</td>
+      <td>Adelie</td>
+      <td>55</td>
+    </tr>
+    <tr>
+      <th>3</th>
+      <td>Torgersen</td>
+      <td>Adelie</td>
+      <td>47</td>
+    </tr>
+    <tr>
+      <th>4</th>
+      <td>Biscoe</td>
+      <td>Adelie</td>
+      <td>44</td>
+    </tr>
+  </tbody>
+</table>
+</div>
+
+
+
+First, we generate a similarity matrix from the penguin features by taking the euclidean distance between each row in the dataset. Then, we convert the distances to similarities by scaling between 0 and 1.
+
+
+```python
+species_order = penguin_features.pop('species')
+distances = pdist(penguin_features)
+penguin_distance_matrix = squareform(distances)
+penguin_similarity_matrix = 1 / (1  + penguin_distance_matrix)
+penguin_similarity_matrix
+```
+
+    array([[1.00000000e+00, 3.30156918e-02, 7.21626767e-04],
+           [3.30156918e-02, 1.00000000e+00, 7.35991958e-04],
+           [7.21626767e-04, 7.35991958e-04, 1.00000000e+00]])
+
+
+Last, we create a metacommunity object from the counts, spcecies order, and similarity matrix.
+
+```python
+penguin_metacommunity = make_metacommunity(
+    penguin_counts, 
+    species_order=species_order, 
+    similarity_matrix=penguin_similarity_matrix
+)
+```
+
+We can calculate all diviersty measures at once for a given viewpoint.
+
+
+```python
+penguin_metacommunity.subcommunities_to_dataframe(viewpoint=0)
+```
+
+<div>
+
+<table border="1" class="dataframe">
+  <thead>
+    <tr style="text-align: right;">
+      <th></th>
+      <th>community</th>
+      <th>viewpoint</th>
+      <th>alpha</th>
+      <th>rho</th>
+      <th>beta</th>
+      <th>gamma</th>
+      <th>normalized_alpha</th>
+      <th>normalized_rho</th>
+      <th>normalized_beta</th>
+    </tr>
+  </thead>
+  <tbody>
+    <tr>
+      <th>0</th>
+      <td>Biscoe</td>
+      <td>0</td>
+      <td>4.068632</td>
+      <td>1.635146</td>
+      <td>0.611566</td>
+      <td>2.638069</td>
+      <td>1.997771</td>
+      <td>0.802886</td>
+      <td>1.245507</td>
+    </tr>
+    <tr>
+      <th>1</th>
+      <td>Dream</td>
+      <td>0</td>
+      <td>5.253772</td>
+      <td>1.735989</td>
+      <td>0.576041</td>
+      <td>3.539491</td>
+      <td>1.934772</td>
+      <td>0.639301</td>
+      <td>1.564208</td>
+    </tr>
+    <tr>
+      <th>2</th>
+      <td>Torgersen</td>
+      <td>0</td>
+      <td>7.106383</td>
+      <td>3.155993</td>
+      <td>0.316858</td>
+      <td>2.251711</td>
+      <td>1.000000</td>
+      <td>0.444107</td>
+      <td>2.251711</td>
+    </tr>
+  </tbody>
+</table>
+</div>
+
+
+
+We can change the viewpoint parameter and recalculate diversity measures.
+
+
+```python
+penguin_metacommunity.subcommunities_to_dataframe(viewpoint=2)
+```
+
+<div>
+
+<table border="1" class="dataframe">
+  <thead>
+    <tr style="text-align: right;">
+      <th></th>
+      <th>community</th>
+      <th>viewpoint</th>
+      <th>alpha</th>
+      <th>rho</th>
+      <th>beta</th>
+      <th>gamma</th>
+      <th>normalized_alpha</th>
+      <th>normalized_rho</th>
+      <th>normalized_beta</th>
+    </tr>
+  </thead>
+  <tbody>
+    <tr>
+      <th>0</th>
+      <td>Biscoe</td>
+      <td>2</td>
+      <td>3.351522</td>
+      <td>1.233509</td>
+      <td>0.810695</td>
+      <td>2.615200</td>
+      <td>1.645657</td>
+      <td>0.605675</td>
+      <td>1.651050</td>
+    </tr>
+    <tr>
+      <th>1</th>
+      <td>Dream</td>
+      <td>2</td>
+      <td>5.202916</td>
+      <td>1.424577</td>
+      <td>0.701963</td>
+      <td>3.132199</td>
+      <td>1.916044</td>
+      <td>0.524620</td>
+      <td>1.906143</td>
+    </tr>
+    <tr>
+      <th>2</th>
+      <td>Torgersen</td>
+      <td>2</td>
+      <td>7.106383</td>
+      <td>3.155993</td>
+      <td>0.316858</td>
+      <td>2.251711</td>
+      <td>1.000000</td>
+      <td>0.444107</td>
+      <td>2.251711</td>
+    </tr>
+  </tbody>
+</table>
+</div>
+
+
+
+We can also compute all metacommunity diversity measures at once. 
+
+
+```python
+penguin_metacommunity.metacommunity_to_dataframe(viewpoint=0)
+```
+
+<div>
+
+<table border="1" class="dataframe">
+  <thead>
+    <tr style="text-align: right;">
+      <th></th>
+      <th>community</th>
+      <th>viewpoint</th>
+      <th>alpha</th>
+      <th>rho</th>
+      <th>beta</th>
+      <th>gamma</th>
+      <th>normalized_alpha</th>
+      <th>normalized_rho</th>
+      <th>normalized_beta</th>
+    </tr>
+  </thead>
+  <tbody>
+    <tr>
+      <th>0</th>
+      <td>metacommunity</td>
+      <td>0</td>
+      <td>4.932543</td>
+      <td>1.886294</td>
+      <td>0.557012</td>
+      <td>2.915662</td>
+      <td>1.834166</td>
+      <td>0.692157</td>
+      <td>1.504464</td>
+    </tr>
+  </tbody>
+</table>
+</div>
+
+
+
+If we only want to compute a single measure, we call them individually.
+
+
+```python
+penguin_metacommunity.subcommunity_alpha(viewpoint=0)
+```
+
+    array([4.06863173, 5.25377155, 7.10638298])
+
+
+
+#### Similarity from file
+
+For large datasets, the similarity matrix may not fit in RAM. To avoid loading the entire matrix into RAM, use the similarity_matrix_filepath argument to read a file from a hard disk drive.
+
+This example calculates diversity indices for the gut microbiomes of four cyclists. For our example, we use a subsample of this dataset. The full dataset can be found [here](https://trace.ncbi.nlm.nih.gov/Traces/sra/?study=SRP219106). See [Peterson, et al., 2017](https://microbiomejournal.biomedcentral.com/articles/10.1186/s40168-017-0320-4) for more details.
+
+
+```python
+# read the files from gist
+metagenome_similarity_filepath = gist_url + 'metagenome_similarity_matrix.csv'
+metagenome_counts_filepath = gist_url + 'metagenome_counts.csv'
+metagenome_similarity = read_csv(metagenome_similarity_filepath)
+metagenome_counts = read_csv(metagenome_counts_filepath)
+
+# save the similarity matrix file locally
+similarity_filepath = 'metagenome_similarity_matrix.csv'
+metagenome_similarity.to_csv(similarity_filepath, index=None)
+```
+
+
+```python
+metagenome = make_metacommunity(
+    metagenome_counts, 
+    similarity_matrix_filepath=similarity_filepath
+)
+```
+
+
+```python
+metagenome.subcommunities_to_dataframe(viewpoint=0)
+```
+
+
+
+
+<div>
+<table border="1" class="dataframe">
+  <thead>
+    <tr style="text-align: right;">
+      <th></th>
+      <th>community</th>
+      <th>viewpoint</th>
+      <th>alpha</th>
+      <th>rho</th>
+      <th>beta</th>
+      <th>gamma</th>
+      <th>normalized_alpha</th>
+      <th>normalized_rho</th>
+      <th>normalized_beta</th>
+    </tr>
+  </thead>
+  <tbody>
+    <tr>
+      <th>0</th>
+      <td>Breezer</td>
+      <td>0</td>
+      <td>5.620958</td>
+      <td>3.878248</td>
+      <td>0.257848</td>
+      <td>1.452272</td>
+      <td>1.404630</td>
+      <td>0.969142</td>
+      <td>1.031841</td>
+    </tr>
+    <tr>
+      <th>1</th>
+      <td>Commencal</td>
+      <td>0</td>
+      <td>5.765467</td>
+      <td>3.971145</td>
+      <td>0.251817</td>
+      <td>1.452311</td>
+      <td>1.441604</td>
+      <td>0.992950</td>
+      <td>1.007101</td>
+    </tr>
+    <tr>
+      <th>2</th>
+      <td>IronHorse</td>
+      <td>0</td>
+      <td>5.613436</td>
+      <td>3.934638</td>
+      <td>0.254153</td>
+      <td>1.426263</td>
+      <td>1.403254</td>
+      <td>0.983586</td>
+      <td>1.016688</td>
+    </tr>
+    <tr>
+      <th>3</th>
+      <td>Scott</td>
+      <td>0</td>
+      <td>5.398303</td>
+      <td>3.821598</td>
+      <td>0.261671</td>
+      <td>1.410603</td>
+      <td>1.350040</td>
+      <td>0.955728</td>
+      <td>1.046323</td>
+    </tr>
+  </tbody>
+</table>
+</div>
+
+
+
+#### Similarity from function
+
+If the similarity matrix file doesn't fit in RAM and we don't have a precomputed similarity matrix, we can define a custom similarity function to generate the similarity matrix on the fly. If the matrix is large or the similarity function is computationally expensive, then this approach can take a long time.
+
+
+```python
+def euclidean_similarity(a, b):
+    euclidean_distance = norm(a - b)
+    return 1 / (1 + euclidean_distance)
+```
+
+
+```python
+penguin_metacommunity = make_metacommunity(
+    penguin_counts, 
+    species_order=species_order, 
+    similarity_function=euclidean_similarity, 
+    features=penguin_features
+)
+```
+
+
+```python
+penguin_metacommunity.subcommunities_to_dataframe(viewpoint=0)
+```
+
+
+
+
+<div>
+<table border="1" class="dataframe">
+  <thead>
+    <tr style="text-align: right;">
+      <th></th>
+      <th>community</th>
+      <th>viewpoint</th>
+      <th>alpha</th>
+      <th>rho</th>
+      <th>beta</th>
+      <th>gamma</th>
+      <th>normalized_alpha</th>
+      <th>normalized_rho</th>
+      <th>normalized_beta</th>
+    </tr>
+  </thead>
+  <tbody>
+    <tr>
+      <th>0</th>
+      <td>Biscoe</td>
+      <td>0</td>
+      <td>4.068632</td>
+      <td>1.635146</td>
+      <td>0.611566</td>
+      <td>2.638069</td>
+      <td>1.997771</td>
+      <td>0.802886</td>
+      <td>1.245507</td>
+    </tr>
+    <tr>
+      <th>1</th>
+      <td>Dream</td>
+      <td>0</td>
+      <td>5.253772</td>
+      <td>1.735989</td>
+      <td>0.576041</td>
+      <td>3.539491</td>
+      <td>1.934772</td>
+      <td>0.639301</td>
+      <td>1.564208</td>
+    </tr>
+    <tr>
+      <th>2</th>
+      <td>Torgersen</td>
+      <td>0</td>
+      <td>7.106383</td>
+      <td>3.155993</td>
+      <td>0.316858</td>
+      <td>2.251711</td>
+      <td>1.000000</td>
+      <td>0.444107</td>
+      <td>2.251711</td>
+    </tr>
+  </tbody>
+</table>
+</div>
+
+
+
+
+```python
+
+```
+
+### Command line interface
+
+For calculating diversity indices from the command line, we will need the following files, which can be downloaded from this [Gist](https://gist.github.com/Elliot-D-Hill/1200c74c83a4be1f63f9ef56b333bb99).
+
+`metagenome_counts.csv`
+```
+subcommunity,species,counts
+subcommunity_1,Akkermansia muciniphila,4394
+subcommunity_1,Faecalibacterium prausnitzii,1807
+...
+subcommunity_4,Blautia coccoides,1
+```
+
+and
+
+`metagenome_similarity_matrix.csv`
+```
+Akkermansia muciniphila,Faecalibacterium prausnitzii,...,Shigella dysenteriae
+1.0,0.6671415709437735,...,0.6093549494839958
+...
+0.6093549494839958,0.6185898023961985,...,1.0
+```
+
+Once the files are downloaded, execute
+```bash
+python -m diversity -i metagenome_counts.csv -s metagenome_similarity_matrix.csv
+```
+
+This is what is printed to the terminal window:
+```
+2022-01-19T14:02:44-0500 (MainProcess, MainThread) INFO python3.9.4 /Users/jasper/diversity/src/diversity/__main__.py -i metagenome_counts.csv -s metagenome_similarity_matrix.csv -v 0 1 2 inf
+community viewpoint alpha rho beta gamma normalized_alpha normalized_rho normalized_beta
+subcommunity_1 0.00 5.6056 3.8792 0.2578 1.4479 1.4008 0.9694 1.0316
+subcommunity_2 0.00 5.7490 3.9715 0.2518 1.4480 1.4375 0.9930 1.0070
+subcommunity_3 0.00 5.5985 3.9352 0.2541 1.4222 1.3995 0.9837 1.0165
+subcommunity_4 0.00 5.3856 3.8230 0.2616 1.4068 1.3469 0.9561 1.0459
+metacommunity 0.00 5.5847 3.9022 0.2563 1.4312 1.3962 0.9756 1.0253
+subcommunity_1 1.00 5.6017 3.8719 0.2583 1.4468 1.3998 0.9675 1.0335
+subcommunity_2 1.00 5.7449 3.9703 0.2519 1.4470 1.4365 0.9927 1.0073
+subcommunity_3 1.00 5.5895 3.9347 0.2541 1.4206 1.3973 0.9836 1.0167
+subcommunity_4 1.00 5.3633 3.8167 0.2620 1.4052 1.3413 0.9545 1.0477
+metacommunity 1.00 5.5732 3.8980 0.2565 1.4298 1.3933 0.9745 1.0262
+subcommunity_1 2.00 5.5979 3.8646 0.2588 1.4456 1.3989 0.9657 1.0355
+subcommunity_2 2.00 5.7410 3.9692 0.2519 1.4459 1.4355 0.9925 1.0076
+subcommunity_3 2.00 5.5808 3.9342 0.2542 1.4189 1.3951 0.9835 1.0168
+subcommunity_4 2.00 5.3421 3.8106 0.2624 1.4037 1.3360 0.9530 1.0493
+metacommunity 2.00 5.5617 3.8937 0.2568 1.4283 1.3904 0.9734 1.0271
+subcommunity_1 inf 5.2934 3.5613 0.2808 1.3452 1.3228 0.8899 1.1237
+subcommunity_2 inf 5.4517 3.8544 0.2594 1.3452 1.3632 0.9638 1.0376
+subcommunity_3 inf 5.2445 3.8535 0.2595 1.3452 1.3110 0.9633 1.0381
+subcommunity_4 inf 4.9728 3.6221 0.2761 1.3452 1.2436 0.9058 1.1040
+metacommunity inf 4.9728 3.5613 0.2594 1.3452 1.2436 0.8899 1.0376
+2022-01-19T14:02:44-0500 (MainProcess, MainThread) INFO Done!
+2022-01-19T14:02:44-0500 (MainProcess, MainThread) INFO process shutting down
+```
+
+Notes:
+- you can use .csv or .tsv for input files, but output is tab-delimited
+- piping output will not include log statements
+- for further options execute `python -m diversity -h`
+
+
+## Background
+
+### Diversity indices
+
+A ***community*** is a collection of elements called ***individuals***, each of which is assigned a label called its ***species***, where multiple individuals may have the same species. A ***diversity index*** is a statistic associated with a community, which describes how much the species of its individuals vary. For example, a community of many individuals of the same species has a very low diversity whereas a community with multiple species and the same amount of individuals per species has a high diversity.
+
+### Partitioned diversity
+
+Some diversity indices compare the diversities of subsets of a community with respect to the overall community. The subsets are called ***subcommunities***, while the overall community is called a ***metacommunity***. For example, two communities with the same frequency distribution across two disjoint sets of the same number of unique species each comprise half of the combined metacommunity diversity. Examples of interpretations of some of these diversity indices include: the "representativeness" of a subcommunity of its metacommunity, and the "contribution" of a subcommunity to its metacommunity's diversity. 
+
+### Frequency-sensitive diversity
+
+[In 1973, Hill introduced a framework](https://doi.org/10.2307/1934352) which unifies commonly used diversity indices into a single parameterized family of diversity measures. The so-called ***viewpoint parameter*** can be thought of as the sensitivity to rare species. At one end of the spectrum, when the viewpoint parameter is set to 0, species frequency is ignored entirely, and only the number of distinct species matters, while at the other end of the spectrum, when the viewpoint parameter is set to ∞, only the highest frequency species in a community is considered by the corresponding diversity measure. Common diversity measures such as ***species richness***, ***Shannon entropy***, the ***Gini-Simpson index***, and the ***Berger-Parker index*** have simple and natural relationships with Hill's indices at different values for the viewpoint parameter (0, 1, 2, ∞, respectively).
+
+### Similarity-sensitive diversity
+
+In addition to being sensitive to frequency, it often makes sense to account for similarity in a diversity measure. For example, a community of 2 different types of rodents, may be considered less diverse as the same community where one of the rodent types was replaced by the same number of individuals of a bird species. [Reeve et al.](https://arxiv.org/abs/1404.6520) and [Leinster and Cobbold](https://doi.org/10.1890/10-2402.1) present a general mathematically rigorous way of incorporating similarity measures into Hill's framework. The result is a family of similarity-sensitive diversity indices parameterized by the same viewpoint parameter as well as the similarity function used for the species in the meta- or subcommunities of interest. As opposed to accounting for distinct species and their frequency, these similarity-sensitive diversity measures can be interpreted as accounting for different possibly overlapping clusters of mutually similar species and their combined frequencies.
+
+### One package to calculate them all
+
+The `diversity` package is able to calculate all of the similarity- and frequency-sensitive subcommunity and metacommunity diversity measures described in [Reeve et al.](https://arxiv.org/abs/1404.6520). See the paper for more in-depth information on their derivation and interpretation.
+
+## Alternatives
+
+To date, we know of no other python package that implements the similarity-sensitive diversity measures calculated by `diversity`. An [R package](https://github.com/boydorr/rdiversity) and a [julia package](https://github.com/EcoJulia/Diversity.jl) exist. However, both packages require the species similarities to be stored in the form of a matrix in memory. That approach does not scale to the amount of species in some applications, such as immune repertoires. `diversity` allows the user to store the similarity matrix in a file, or simply provide a python function that computes the similarities on the fly.
diff --git a/__main__.py b/__main__.py
deleted file mode 100644
index b871068d..00000000
--- a/__main__.py
+++ /dev/null
@@ -1,60 +0,0 @@
-from parameters import configure_arguments
-from log import LOG_HANDLER, LOGGER
-from diversity import Metacommunity
-from numpy import genfromtxt
-from sys import argv
-from platform import python_version
-from logging import captureWarnings, getLogger
-
-
-# Ensure warnings are handled properly.
-captureWarnings(True)
-getLogger('py.warnings').addHandler(LOG_HANDLER)
-
-########################################################################
-
-
-def main():
-    # Parse and validate arguments
-    parser = configure_arguments()
-    args = parser.parse_args()
-
-    LOGGER.setLevel(args.log_level)
-    LOGGER.info(' '.join([f'python{python_version()}', *argv]))
-    LOGGER.debug(f'args: {args}')
-
-    # FIXME equal species are listed multiple times (once for each
-    # subcommunity they are members of, which is deceiving
-    data = genfromtxt(args.filepath, delimiter=',', dtype=object)
-    LOGGER.debug(f'data: {data}')
-    counts = data[:, :3]
-    features = data[:, 3:]
-    viewpoint = args.q[0]
-    meta = Metacommunity(counts, viewpoint, args.Z, features=features)
-
-    print('\n')
-
-    print(meta.alpha)
-    print(meta.rho)
-    print(meta.beta)
-    print(meta.gamma)
-    print(meta.normalized_alpha)
-    print(meta.normalized_rho)
-    print(meta.normalized_beta)
-
-    print('\n')
-
-    print(meta.A)
-    print(meta.R)
-    print(meta.B)
-    print(meta.G)
-    print(meta.normalized_A)
-    print(meta.normalized_R)
-    print(meta.normalized_B)
-
-    LOGGER.info('Done!')
-
-
-########################################################################
-if __name__ == "__main__":
-    main()
diff --git a/diversity.py b/diversity.py
deleted file mode 100644
index 37218643..00000000
--- a/diversity.py
+++ /dev/null
@@ -1,259 +0,0 @@
-"""Module for calculating metacommunity and subcommunity diversity measures.
-
-Classes
--------
-
-Abundance
-    Represents the relative abundances or the metacommunity and its subcommunities
-
-
-Similarity
-    Represents the similarities between species weighted by their relative abundances
-
-Metacommunity
-    Represents a metacommunity made up of subcommunities and computes metacommunity
-    subcommunity diversity measures.
-"""
-from csv import reader, writer
-from pathlib import Path
-from dataclasses import dataclass, field
-from functools import cached_property
-from typing import Callable
-from pandas import DataFrame
-from numpy import (amin, dot, array, empty, unique, isclose,
-                   prod, zeros, sum as numpy_sum, broadcast_to, power,
-                   multiply, divide, float64, inf)
-
-########################################################################
-
-
-@dataclass
-class Abundance:
-
-    counts: array
-
-    @cached_property
-    def p(self):
-        """
-        Calculates the relative abundance of each species in the metacommunity
-        """
-        return self.P.sum(axis=1, keepdims=True)
-
-    @cached_property
-    def P(self):
-        """Calculates the relative abundance of each species in each subcommunity.
-
-        Returns
-        -------
-        A 2D numpy.ndarray where rows are unique species and columns are subcommunities
-        and each element is a species count in a subcommunity
-        """
-        rows, row_pos = unique(self.counts[:, 0], return_inverse=True)
-        cols, col_pos = unique(self.counts[:, 2], return_inverse=True)
-        metacommunity_counts = zeros(
-            (len(rows), len(cols)), dtype=float64)
-        metacommunity_counts[row_pos, col_pos] = self.counts[:, 1]
-        total_abundance = metacommunity_counts.sum()
-        return metacommunity_counts / total_abundance
-
-    @cached_property
-    def w(self):
-        """
-        Calculates each subcommunity's relative abundance as a proportion of the 
-        total metacommunity's relative abundance
-        """
-        return self.P.sum(axis=0)
-
-    @cached_property
-    def normalized_P(self):
-        """
-        Calculates normalized subcommunity relative abundances
-        """
-        return self.P / self.w
-
-
-@dataclass
-class Similarity:
-
-    abundance: Abundance
-    Z: array = None
-    z_filepath: str = None
-    similarity_fn: Callable = None
-    features: array = None
-
-    @cached_property
-    def Zp(self):
-        return self.calculate_zp(self.abundance.p)
-
-    @cached_property
-    def ZP(self):
-        return self.calculate_zp(self.abundance.P)
-
-    @cached_property
-    def normalized_ZP(self):
-        return self.calculate_zp(self.abundance.normalized_P)
-
-    def write_similarity_matrix(self):
-        n_species = self.features.shape[0]
-        z_i = empty(n_species, dtype=float64)
-        with open(self.z_filepath, 'w') as f:
-            csv_writer = writer(f)
-            for species_i in self.features:
-                for j, species_j in enumerate(self.features):
-                    z_i[j] = self.similarity_fn(species_i, species_j)
-                csv_writer.writerow(z_i)
-
-    def zp_from_file(self, P):
-        ZP = empty(P.shape, dtype=float64)
-        with open(self.z_filepath, 'r') as f:
-            for i, z_i in enumerate(reader(f)):
-                z_i = array(z_i, dtype=float64)
-                ZP[i, :] = dot(z_i, P)
-        return ZP
-
-    def zp_from_array(self, P):
-        return dot(self.Z, P)
-
-    def calculate_zp(self, P):
-        if self.Z is not None:
-            return self.zp_from_array(P)
-        if not self.z_filepath.is_file():
-            self.write_similarity_matrix()
-        return self.zp_from_file(P)
-
-
-@dataclass
-class Metacommunity:
-    """Class for metacommunities and calculation their diversity.
-
-    Attributes
-    ----------
-    counts: numpy.ndarray
-    _viewpoint: float
-    z_filepath: str
-    similarity_fn: Callable
-    features: np.ndarray
-    abundance: Abundance
-    similarity: Similarity
-    """
-
-    counts: array
-    _viewpoint: float
-    Z: array = None
-    z_filepath: str = None
-    similarity_fn: Callable = None
-    features: array = None
-    abundance: Abundance = field(init=False)
-    similarity: Similarity = field(init=False)
-
-    def __post_init__(self):
-        self.abundance = Abundance(self.counts)
-        if self.z_filepath:
-            self.z_filepath = Path(self.z_filepath)
-        self.similarity = Similarity(
-            self.abundance,
-            Z=self.Z,
-            z_filepath=self.z_filepath,
-            similarity_fn=self.similarity_fn,
-            features=self.features)
-
-    # FIXME validate new viewpoint
-
-    def set_viewpoint(self, viewpoint):
-        self._viewpoint = viewpoint
-
-    @property
-    def alpha(self):
-        return self.subcommunity_measure(1, self.similarity.ZP)
-
-    @property
-    def rho(self):
-        return self.subcommunity_measure(self.similarity.Zp, self.similarity.ZP)
-
-    @property
-    def beta(self):
-        return 1 / self.rho
-
-    @property
-    def gamma(self):
-        return self.subcommunity_measure(1, self.similarity.Zp)
-
-    @property
-    def normalized_alpha(self):
-        return self.subcommunity_measure(1, self.similarity.normalized_ZP)
-
-    @property
-    def normalized_rho(self):
-        return self.subcommunity_measure(self.similarity.Zp, self.similarity.normalized_ZP)
-
-    @property
-    def normalized_beta(self):
-        return 1 / self.normalized_rho
-
-    @property
-    def A(self):
-        return self.metacommunity_measure(self.alpha)
-
-    @property
-    def R(self):
-        return self.metacommunity_measure(self.rho)
-
-    @property
-    def B(self):
-        return self.metacommunity_measure(self.beta)
-
-    @property
-    def G(self):
-        return self.metacommunity_measure(self.gamma)
-
-    @property
-    def normalized_A(self):
-        return self.metacommunity_measure(self.normalized_alpha)
-
-    @property
-    def normalized_R(self):
-        return self.metacommunity_measure(self.normalized_rho)
-
-    @property
-    def normalized_B(self):
-        return self.metacommunity_measure(self.normalized_beta)
-
-    def subcommunity_measure(self, numerator, denominator):
-        similarities = divide(numerator, denominator, out=zeros(
-            denominator.shape), where=denominator != 0)
-        return self.power_mean(self.abundance.normalized_P, similarities)
-
-    def metacommunity_measure(self, subcommunity_measure):
-        return self.power_mean(self.abundance.w, subcommunity_measure)
-
-    def power_mean(self, weights, items):
-        """Calculates a weighted power mean.
-
-        Parameters
-        ----------
-        weights: numpy.ndarray
-            The weights corresponding to items.
-        items: numpy.ndarray
-            The elements for which the weighted power mean is computed.
-
-        Returns
-        -------
-        The power mean of items with exponent order, weighted by weights.
-        When order is close to 1 or less than -100, analytical formulas
-        for the limits at 1 and -infinity are used respectively.
-        """
-        order = 1 - self._viewpoint
-        mask = weights != 0
-        if isclose(order, 0):
-            return prod(power(items, weights, where=mask), axis=0, where=mask)
-        elif order < -100:
-            items = broadcast_to(items, weights.shape)
-            return amin(items, axis=0, where=mask, initial=inf)
-        items_power = power(items, order, where=mask)
-        items_product = multiply(items_power, weights, where=mask)
-        items_sum = numpy_sum(items_product, axis=0, where=mask)
-        return power(items_sum, 1 / order)
-
-    # FIXME implement me!
-    def format_results(self):
-        pass  # DataFrame() maybe?
diff --git a/parameters.py b/parameters.py
deleted file mode 100644
index c1804e79..00000000
--- a/parameters.py
+++ /dev/null
@@ -1,79 +0,0 @@
-"""Module for Chubacabra's argument parser configuration.
-
-Classes
--------
-ArgumentWarning
-    Warning category for problematic arguments.
-ValidateQ
-    Validator for -q parameter.
-
-Functions
----------
-configure_arguments
-    Creates argument parser.
-"""
-from argparse import Action, ArgumentParser
-from warnings import warn
-
-########################################################################
-
-
-class ArgumentWarning(Warning):
-    """Used for warnings related to problematic argument choices."""
-    pass
-
-########################################################################
-
-
-class ValidateQ(Action):
-    """Validator for -q parameter."""
-
-    def __call__(self, parser, args, values, option_string=None):
-        """Validates -q parameter.
-
-        Warns if arguments larger than 100 are passed, reminding that
-        they are treated as infinity in the diversity calculation.
-        """
-        if any([q > 100 for q in values]):
-            warn("q > 100.0 defaults to the analytical formula for q = inf.",
-                 category=ArgumentWarning)
-        setattr(args, self.dest, values)
-
-########################################################################
-
-
-def configure_arguments():
-    """Creates argument parser.
-
-    Returns
-    -------
-    argparse.ArgumentParser configured to handle command-line arguments
-    for Chubacabra.
-    """
-    parser = ArgumentParser()
-    parser.add_argument(
-        "filepath",
-        type=str,
-        help=("A csv file where the first 3 columns of the file are the"
-              " species name, its count, and subcommunity name, and all"
-              " following columns are features of that species that"
-              " will be used to calculate similarity between species."))
-    parser.add_argument(
-        "-l", "--log_level",
-        help=("Logging verbosity level. Must be one of DEBUG, INFO,"
-              " WARNING, ERROR, CRITICAL (listed in decreasing"
-              " verbosity)."),
-        default="INFO")
-    parser.add_argument(
-        "-q",
-        nargs='+',
-        type=float,
-        help="A list of viewpoint parameters.",
-        action=ValidateQ)
-    parser.add_argument(
-        "Z",
-        type=str,
-        help=("The filepath to a csv file containing a symmetric"
-              " similarity matrix. If the file does not exist, one will"
-              " be created with the user defined similarity function."))
-    return parser
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 00000000..374b58cb
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,6 @@
+[build-system]
+requires = [
+    "setuptools>=42",
+    "wheel"
+]
+build-backend = "setuptools.build_meta"
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 00000000..8dbcd3f4
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,2 @@
+numpy==1.22.0
+pandas==1.3.5
diff --git a/requirements_dev.txt b/requirements_dev.txt
new file mode 100644
index 00000000..fe9151a7
--- /dev/null
+++ b/requirements_dev.txt
@@ -0,0 +1,3 @@
+pytest==6.2.5
+tox==3.24.5
+tox-gh-actions==2.9.1
diff --git a/setup.cfg b/setup.cfg
new file mode 100644
index 00000000..229cc84d
--- /dev/null
+++ b/setup.cfg
@@ -0,0 +1,24 @@
+[metadata]
+name = diversity
+version = 1.0.0
+author = Elliot Hill, Jasper Braun
+author_email = elliot.douglas.hill@gmail.com
+description = Calculations of diversity with similarity
+long_description = file: README.md
+long_description_content_type = text/markdown
+url = https://github.com/Elliot-D-Hill//metacommunity
+project_urls =
+    Bug Tracker = https://github.com/Elliot-D-Hill/metacommunity/issues
+classifiers =
+    Programming Language :: Python :: 3
+    License :: OSI Approved :: MIT License
+    Operating System :: OS Independent
+
+[options]
+package_dir =
+    = src
+packages = find:
+python_requires = >=3.8
+
+[options.packages.find]
+where = src
diff --git a/__init__.py b/src/diversity/__init__.py
similarity index 100%
rename from __init__.py
rename to src/diversity/__init__.py
diff --git a/src/diversity/__main__.py b/src/diversity/__main__.py
new file mode 100644
index 00000000..4b572a1b
--- /dev/null
+++ b/src/diversity/__main__.py
@@ -0,0 +1,70 @@
+"""Main module for executing diversity on command-line.
+
+Functions
+---------
+main
+    Calculates diversities according to command-line specifications.
+"""
+from sys import argv
+from platform import python_version
+from logging import captureWarnings, getLogger
+
+from pandas import read_csv, concat
+
+from diversity.metacommunity import make_metacommunity
+from diversity.parameters import configure_arguments
+from diversity.log import LOG_HANDLER, LOGGER
+from diversity.utilities import get_file_delimiter
+
+# Ensure warnings are handled properly.
+captureWarnings(True)
+getLogger("py.warnings").addHandler(LOG_HANDLER)
+
+
+def main(args):
+    """Calculates diversity from species counts and similarities.
+
+    Parameters
+    ----------
+    args: argparse.Namespace
+        Return object of argparse.ArgumentParser object created by
+        diversity.parameters.configure_arguments and applied to command
+        line arguments.
+    """
+    LOGGER.setLevel(args.log_level)
+    LOGGER.info(" ".join([f"python{python_version()}", *argv]))
+    LOGGER.debug(f"args: {args}")
+
+    delimiter = get_file_delimiter(args.input_filepath)
+    species_counts = read_csv(
+        args.input_filepath,
+        sep=delimiter,
+        usecols=[args.subcommunity_column, args.species_column, args.count_column],
+    )[[args.subcommunity_column, args.species_column, args.count_column]]
+
+    LOGGER.debug(f"data: {species_counts}")
+
+    meta = make_metacommunity(
+        species_counts, similarity_matrix_filepath=args.similarity_matrix_filepath
+    )
+
+    community_views = []
+    for view in args.viewpoint:
+        community_views.append(meta.subcommunities_to_dataframe(view))
+        community_views.append(meta.metacommunity_to_dataframe(view))
+
+    community_views = concat(community_views, ignore_index=True)
+    community_views.viewpoint = community_views.viewpoint.map(
+        lambda v: format(v, ".2f")
+    )
+    community_views.to_csv(
+        args.output_filepath, sep="\t", float_format="%.4f", index=False
+    )
+
+    LOGGER.info("Done!")
+
+
+if __name__ == "__main__":
+    parser = configure_arguments()
+    args = parser.parse_args()
+    main(args)
diff --git a/log.py b/src/diversity/log.py
similarity index 74%
rename from log.py
rename to src/diversity/log.py
index 43a653f9..f3b1f991 100644
--- a/log.py
+++ b/src/diversity/log.py
@@ -15,14 +15,14 @@
     INFO,
     StreamHandler,
     WARNING,
-    )
+)
 from multiprocessing import get_logger
-from time import (
-    gmtime,
-    strftime)
+from time import gmtime, strftime
 
-date_format = '%Y-%m-%dT%H:%M:%S%z'
-logging_format = '%(asctime)s\t(%(processName)s, %(threadName)s)\t%(levelname)s\t%(message)s'
+date_format = "%Y-%m-%dT%H:%M:%S%z"
+logging_format = (
+    "%(asctime)s\t(%(processName)s, %(threadName)s)\t%(levelname)s\t%(message)s"
+)
 
 LOGGER = get_logger()
 formatter = Formatter(fmt=logging_format, datefmt=date_format)
diff --git a/src/diversity/metacommunity.py b/src/diversity/metacommunity.py
new file mode 100644
index 00000000..a55afbaf
--- /dev/null
+++ b/src/diversity/metacommunity.py
@@ -0,0 +1,785 @@
+"""Module for metacommunity and subcommunity diversity measures.
+
+Classes
+-------
+Abundance
+    Species abundances in metacommunity.
+Similarity
+    Abstract base class for relative abundance-weighted species
+    similarities.
+SimilarityFromFile
+    Implements Similarity by reading similarities from a file.
+SimilarityFromFunction
+    Implements Similarity by calculating similarities with a function.
+SimilarityFromMemory
+    Implements Similarity by storing similarities in memory.
+Metacommunity
+    Represents a metacommunity made up of subcommunities and computes
+    metacommunity subcommunity diversity measures.
+
+Functions
+---------
+make_similarity
+    Returns correct diversity.metacommunity.Similarity object fitting
+    parameter specification.
+make_metacommunity
+    Builds diversity.metacommunity.Metacommunity object according to
+    parameter specification.
+"""
+from abc import ABC, abstractmethod
+from csv import reader
+from functools import cached_property
+from pathlib import Path
+
+from pandas import DataFrame
+from numpy import array, empty, zeros, broadcast_to, divide, float64
+
+from diversity.utilities import (
+    InvalidArgumentError,
+    get_file_delimiter,
+    power_mean,
+    unique_correspondence,
+)
+
+
+class Abundance:
+    """Relative abundances of species in a metacommunity.
+
+    A community consists of a set of species, each of which may appear
+    any (non-negative) number of times. A metacommunity consists of one
+    or more subcommunities and can be represented by the number of
+    appearances of each species in each of the subcommunities that the
+    species appears in.
+    """
+
+    def __init__(
+        self,
+        counts,
+        species_order=None,
+        subcommunity_order=None,
+        subcommunity_column=0,
+        species_column=1,
+        count_column=2,
+    ):
+        """Initializes object.
+
+        Determines species and subcommunity orderings if needed.
+
+        Parameters
+        ----------
+        counts: numpy.ndarray
+            A 2-d numpy.ndarray with subcommunity identifiers, species
+            identifiers and number of appearances of the row's species
+            in the row's subcommunity. The column ordering is determined
+            by the subcommunity_column, species_column, and
+            counts_column parameters. Each combination of species and
+            subcommunity must appear no more than once.
+        species_order: numpy.ndarray
+            Ordered unique species identifiers. The ordering determines
+            in which order values corresponding to each species are
+            returned by methods of the object.
+        subcommunity_order: Iterable
+            Ordered unique subcommunity identifiers. The ordering
+            determines in which order values corresponding to each
+            species are returned by methods of the object.
+        subcommunity_column: int
+            Index of subcommunity identifier column in counts.
+        species_column: int
+            Index of species identifier column in counts.
+        count_column: int
+            Index of species count column in counts.
+        """
+        self.counts = counts
+        self.subcommunity_column = subcommunity_column
+        self.species_column = species_column
+        self.count_column = count_column
+        self.species_order, self.__species_unique_pos = unique_correspondence(
+            items=self.counts[:, self.species_column],
+            ordered_unique_items=species_order,
+        )
+        self.subcommunity_order, self.__subcommunity_unique_pos = unique_correspondence(
+            items=self.counts[:, self.subcommunity_column],
+            ordered_unique_items=subcommunity_order,
+        )
+
+    def __pivot_table(self):
+        """Converts long to wide formatted counts.
+
+        Returns
+        -------
+        A numpy.ndarray where rows correspond to species, columns to
+        subcommunities and each element is the count of a species in a
+        specific subcommunity.
+        """
+        table = zeros(
+            (len(self.species_order), len(self.subcommunity_order)), dtype=float64
+        )
+        table[self.__species_unique_pos, self.__subcommunity_unique_pos] = self.counts[
+            :, self.count_column
+        ].astype(float64)
+        return table
+
+    @cached_property
+    def subcommunity_abundance(self):
+        """Calculates the relative abundances in subcommunities.
+
+        Returns
+        -------
+        A numpy.ndarray of shape (n_species, n_subcommunities), where
+        rows correspond to unique species, columns correspond to
+        subcommunities and each element is the abundance of the species
+        in the subcommunity relative to the total metacommunity size.
+        The row ordering is established by the species_to_row attribute.
+        """
+        metacommunity_counts = self.__pivot_table()
+        total_abundance = metacommunity_counts.sum()
+        return metacommunity_counts / total_abundance
+
+    @cached_property
+    def metacommunity_abundance(self):
+        """Calculates the relative abundances in metacommunity.
+
+        Returns
+        -------
+        A numpy.ndarray of shape (n_species, 1), where rows correspond
+        to unique species and each row contains the relative abundance
+        of the species in the metacommunity. The row ordering is
+        established by the species_to_row attribute.
+        """
+        return self.subcommunity_abundance.sum(axis=1, keepdims=True)
+
+    @cached_property
+    def subcommunity_normalizing_constants(self):
+        """Calculates subcommunity normalizing constants.
+
+        Returns
+        -------
+        A numpy.ndarray of shape (n_subcommunities,), with the fraction
+        of each subcommunity's size of the metacommunity.
+        """
+        return self.subcommunity_abundance.sum(axis=0)
+
+    @cached_property
+    def normalized_subcommunity_abundance(self):
+        """Calculates normalized relative abundances in subcommunities.
+
+        Returns
+        -------
+        A numpy.ndarray of shape (n_species, n_subcommunities), where
+        rows correspond to unique species, columns correspond to
+        subcommunities and each element is the abundance of the species
+        in the subcommunity relative to the subcommunity size. The row
+        ordering is established by the species_to_row attribute.
+        """
+        return self.subcommunity_abundance / self.subcommunity_normalizing_constants
+
+
+class Similarity(ABC):
+    """Interface for classes computing weighted similarities."""
+
+    def __init__(self, species_order):
+        self.species_order = species_order
+
+    @abstractmethod
+    def calculate_weighted_similarities(self, relative_abundances):
+        """Calculates weighted sums of similarities to each species.
+
+        Parameters
+        ----------
+        relative_abundances: numpy.ndarray
+            Array of shape (n_species, n_communities), where rows
+            correspond to unique species, columns correspond to
+            (meta-/sub-)communities and each element is the relative
+            abundance of a species in a (meta-/sub-)community.
+
+        Returns
+        -------
+        A 2-d numpy.ndarray of shape (n_species, n_communities), where
+        rows correspond to unique species, columns correspond to
+        (meta-/sub-) communities and each element is a sum of
+        similarities to one species weighted by the similarities stored
+        in the similarities file.
+        """
+        pass
+
+
+class SimilarityFromFile(Similarity):
+    """Implements Similarity by using similarities stored in file."""
+
+    def __init__(self, similarity_matrix_filepath):
+        """Initializes object.
+
+        similarity_matrix_filepath: str
+            Path to similarities file. If similarity_function is None
+            the file must exist and contain a square matrix of
+            similarities between species, together with a header
+            containing the unique species names in the matrix's row and
+            column ordering. If similarity_function is also specified,
+            the file must not exist, but will instead be generated as
+            soon as it is needed in subsequent computations.
+        """
+        self.similarity_matrix_filepath = Path(similarity_matrix_filepath)
+        self._delimiter = get_file_delimiter(self.similarity_matrix_filepath)
+        self.species_order = self.__get_species_order()
+
+    def __get_species_order(self):
+        """The species ordering used in similarity matrix file."""
+        with open(self.similarity_matrix_filepath, "r") as file:
+            species_order = array(next(reader(file, delimiter=self._delimiter)))
+        return species_order
+
+    def calculate_weighted_similarities(self, relative_abundances):
+        """Calculates weighted sums of similarities to each species.
+
+        Similarities are read from similarities file referred to by
+        object's similarity_matrix_filepath attribute.
+
+        See diversity.metacommunity.Similarity.calculate_weighted_similarities
+        for complete specification.
+        """
+        weighted_similarities = empty(relative_abundances.shape, dtype=float64)
+        with open(self.similarity_matrix_filepath, "r") as file:
+            next(reader(file, delimiter=self._delimiter))
+            for i, row in enumerate(reader(file, delimiter=self._delimiter)):
+                similarities_row = array(row, dtype=float64)
+                weighted_similarities[i, :] = similarities_row @ relative_abundances
+        return weighted_similarities
+
+
+class SimilarityFromFunction(Similarity):
+    """Implements Similarity using a similarity function."""
+
+    def __init__(self, similarity_function, features, species_order):
+        """Initializes object.
+
+        Parameters
+        ----------
+        similarity_function: Callable
+            Callable to determine similarity between species. Must take
+            two items from the features argument and return a numeric
+            similarity value.
+        features: numpy.ndarray
+            A numpy.ndarray where each item along axis 0 comprises the
+            features of a species that are passed to
+            similarity_function. The order of features is determined
+            by species_order.
+        species_order: Iterable
+            The unique species in desired order. Row and column ordering
+            in similarity matrix calculations is determined by this
+            argument.
+        """
+
+        self.similarity_function = similarity_function
+        self.features = features
+        super().__init__(species_order)
+
+    def calculate_weighted_similarities(self, relative_abundances):
+        """Calculates weighted sums of similarities to each species.
+
+        Similarities are calculated using object's similarity_function
+        attribute.
+
+        See diversity.metacommunity.Similarity.calculate_weighted_similarities
+        for complete specification.
+        """
+        weighted_similarities = empty(relative_abundances.shape, dtype=float64)
+        similarities_row = empty(relative_abundances.shape[0], dtype=float64)
+        for i, species_i in enumerate(self.features):
+            for j, species_j in enumerate(self.features):
+                similarities_row[j] = self.similarity_function(species_i, species_j)
+            weighted_similarities[i, :] = similarities_row @ relative_abundances
+        return weighted_similarities
+
+
+class SimilarityFromMemory(Similarity):
+    """Implements Similarity using similarities stored in memory."""
+
+    def __init__(self, similarity_matrix, species_order):
+        """Initializes object.
+
+        similarity_matrix: numpy.ndarray
+            2-d array of similarities between species. Ordering of rows
+            and columns must correspond to species_order argument.
+        species_order: Iterable
+            The unique species in desired order.
+        """
+        self.similarity_matrix = similarity_matrix
+        super().__init__(species_order)
+
+    def calculate_weighted_similarities(self, relative_abundances):
+        """Calculates weighted sums of similarities to each species.
+
+        Similarities are calculated using object's similarity_matrix
+        attribute.
+
+        See diversity.metacommunity.Similarity.calculate_weighted_similarities
+        for complete specification.
+        """
+        return self.similarity_matrix @ relative_abundances
+
+
+def make_similarity(
+    similarity_matrix=None,
+    species_order=None,
+    similarity_matrix_filepath=None,
+    similarity_function=None,
+    features=None,
+):
+    """Creates a Similarity object from specified parameter combination.
+
+    Valid parameter combinations and the returned types:
+        - similarity_matrix, species_order -> SimilarityFromMemory
+        - similarity_matrix_filepath -> SimilarityFromFile
+        - species_order, similarity_matrix_filepath,
+          similarity_function, features -> SimilarityFromFunction
+
+    Parameters
+    ----------
+    similarity_matrix: numpy.ndarray
+        Used by diversity.metacommunity.SimilarityFromMemory.
+    species_order: Iterable
+        Used by diversity.metacommunity.SimilarityFromMemory, or
+        diversity.metacommunity.SimilarityFromFunction.
+    similarity_matrix_filepath: str
+        Used by diversity.metacommunity.SimilarityFromFile, or
+        diversity.metacommunity.SimilarityFromFunction.
+    similarity_function: Callable
+        Used by diversity.metacommunity.SimilarityFromFunction.
+    features: numpy.ndarray
+        Used by diversity.metacommunity.SimilarityFromFunction.
+
+    Returns
+    -------
+    An object whose type is the implementation of the abstract base
+    Similarity according to parameter specification.
+    """
+    from_memory_parameters = (similarity_matrix, species_order)
+    from_function_parameters = (
+        similarity_function,
+        features,
+        species_order,
+    )
+    from_file_parameters = (similarity_matrix_filepath,)
+
+    are_all_not_none = lambda parameters: all(p is not None for p in parameters)
+    if are_all_not_none(from_memory_parameters):
+        return SimilarityFromMemory(*from_memory_parameters)
+    elif are_all_not_none(from_function_parameters):
+        return SimilarityFromFunction(*from_function_parameters)
+    elif are_all_not_none(from_file_parameters):
+        return SimilarityFromFile(*from_file_parameters)
+    else:
+        raise InvalidArgumentError(
+            "Invalid argument combination. See the documentation for"
+            " valid argument combinations."
+        )
+
+
+class Metacommunity:
+    """Class for metacommunities and calculating their diversity.
+
+    All diversities computed by objects of this class are
+    similarity-sensitive. See https://arxiv.org/abs/1404.6520 for
+    precise definitions of the various diversity measures.
+    """
+
+    def __init__(self, similarity, abundance):
+        """Initializes object.
+
+        Parameters
+        ----------
+        similarity: diversity.metacommunity.Similarity
+            Object for calculating abundance-weighted similarities.
+        abundance: diversity.metacommunity.Abundance
+            Object whose (sub-/meta-)community species abundances are
+            used.
+        """
+        self.similarity = similarity
+        self.abundance = abundance
+
+    @cached_property
+    def __metacommunity_similarity(self):
+        """Sums of similarities weighted by metacommunity abundances."""
+        return self.similarity.calculate_weighted_similarities(
+            self.abundance.metacommunity_abundance
+        )
+
+    @cached_property
+    def __subcommunity_similarity(self):
+        """Sums of similarities weighted by subcommunity abundances."""
+        return self.similarity.calculate_weighted_similarities(
+            self.abundance.subcommunity_abundance
+        )
+
+    @cached_property
+    def __normalized_subcommunity_similarity(self):
+        """Sums of similarities weighted by the normalized subcommunity abundances."""
+        return self.similarity.calculate_weighted_similarities(
+            self.abundance.normalized_subcommunity_abundance
+        )
+
+    def subcommunity_alpha(self, viewpoint):
+        """Calculates alpha class diversities of subcommunities.
+
+        Corresponds roughly to the diversities of subcommunities
+        relative to the metacommunity.
+
+        Parameters
+        ----------
+        viewpoint: numeric
+            Non-negative number. Can be interpreted as the degree of
+            ignorance towards rare species, where 0 treats rare species
+            the same as frequent species, and infinity considers only the
+            most frequent species.
+        """
+        return self.__subcommunity_measure(viewpoint, 1, self.__subcommunity_similarity)
+
+    def subcommunity_rho(self, viewpoint):
+        """Calculates rho class diversities of subcommunities.
+
+        Corresponds roughly to how redundant each subcommunity's classes
+        are in the metacommunity.
+
+        Parameters
+        ----------
+        viewpoint: numeric
+            Non-negative number. Can be interpreted as the degree of
+            ignorance towards rare species, where 0 treats rare species
+            the same as frequent species, and infinity considers only the
+            most frequent species.
+        """
+        return self.__subcommunity_measure(
+            viewpoint,
+            self.__metacommunity_similarity,
+            self.__subcommunity_similarity,
+        )
+
+    def subcommunity_beta(self, viewpoint):
+        """Calculates beta class diversities of subcommunities.
+
+        Corresponds roughly to how distinct each subcommunity's classes
+        are from all classes in metacommunity.
+
+        Parameters
+        ----------
+        viewpoint: numeric
+            Non-negative number. Can be interpreted as the degree of
+            ignorance towards rare species, where 0 treats rare species
+            the same as frequent species, and infinity considers only
+            the most frequent species.
+        """
+        return 1 / self.subcommunity_rho(viewpoint)
+
+    def subcommunity_gamma(self, viewpoint):
+        """Calculates gamma class diversities of subcommunities.
+
+        Corresponds roughly to how much each subcommunity contributes
+        towards the metacommunity diversity.
+
+        Parameters
+        ----------
+        viewpoint: numeric
+            Non-negative number. Can be interpreted as the degree of
+            ignorance towards rare species, where 0 treats rare species
+            the same as frequent species, and infinity considers only
+            the most frequent species.
+        """
+        denominator = broadcast_to(
+            self.__metacommunity_similarity,
+            self.abundance.normalized_subcommunity_abundance.shape,
+        )
+        return self.__subcommunity_measure(viewpoint, 1, denominator)
+
+    def normalized_subcommunity_alpha(self, viewpoint):
+        """Calculates normalized alpha class diversities of subcommunities.
+
+        Corresponds roughly to the diversities of subcommunities in
+        isolation.
+
+        Parameters
+        ----------
+        viewpoint: numeric
+            Non-negative number. Can be interpreted as the degree of
+            ignorance towards rare species, where 0 treats rare species
+            the same as frequent species, and infinity considers only the
+            most frequent species.
+        """
+        return self.__subcommunity_measure(
+            viewpoint, 1, self.__normalized_subcommunity_similarity
+        )
+
+    def normalized_subcommunity_rho(self, viewpoint):
+        """Calculates normalized rho class diversities of subcommunities.
+
+        Corresponds roughly to the representativeness of subcommunities.
+
+        Parameters
+        ----------
+        viewpoint: numeric
+            Non-negative number. Can be interpreted as the degree of
+            ignorance towards rare species, where 0 treats rare species
+            the same as frequent species, and infinity considers only the
+            most frequent species.
+        """
+        return self.__subcommunity_measure(
+            viewpoint,
+            self.__metacommunity_similarity,
+            self.__normalized_subcommunity_similarity,
+        )
+
+    def normalized_subcommunity_beta(self, viewpoint):
+        """Calculates normalized rho class diversities of subcommunities.
+
+        Corresponds roughly to average diversity of subcommunities in
+        the metacommunity.
+
+        Parameters
+        ----------
+        viewpoint: numeric
+            Non-negative number. Can be interpreted as the degree of
+            ignorance towards rare species, where 0 treats rare species
+            the same as frequent species, and infinity considers only the
+            most frequent species.
+        """
+        return 1 / self.normalized_subcommunity_rho(viewpoint)
+
+    def metacommunity_alpha(self, viewpoint):
+        """Calculates alpha class diversity of metacommunity.
+
+        Corresponds roughly to the average diversity of subcommunities
+        relative to the metacommunity.
+
+        Parameters
+        ----------
+        viewpoint: numeric
+            Non-negative number. Can be interpreted as the degree of
+            ignorance towards rare species, where 0 treats rare species
+            the same as frequent species, and infinity considers only the
+            most frequent species.
+        """
+        return self.__metacommunity_measure(viewpoint, self.subcommunity_alpha)
+
+    def metacommunity_rho(self, viewpoint):
+        """Calculates rho class diversitiy of metacommunity.
+
+        Corresponds roughly to the average redundancy of subcommunities
+        in the metacommunity.
+
+        Parameters
+        ----------
+        viewpoint: numeric
+            Non-negative number. Can be interpreted as the degree of
+            ignorance towards rare species, where 0 treats rare species
+            the same as frequent species, and infinity considers only the
+            most frequent species.
+        """
+        return self.__metacommunity_measure(viewpoint, self.subcommunity_rho)
+
+    def metacommunity_beta(self, viewpoint):
+        """Calculates beta class diversity of metacommunity.
+
+        Corresponds roughly to the average distinctness of
+        subcommunities within the metacommunity.
+
+        Parameters
+        ----------
+        viewpoint: numeric
+            Non-negative number. Can be interpreted as the degree of
+            ignorance towards rare species, where 0 treats rare species
+            the same as frequent species, and infinity considers only the
+            most frequent species.
+        """
+        return self.__metacommunity_measure(viewpoint, self.subcommunity_beta)
+
+    def metacommunity_gamma(self, viewpoint):
+        """Calculates gamma class diversity of metacommunity.
+
+        Corresponds roughly to the class diversity of the unpartitioned
+        metacommunity.
+
+        Parameters
+        ----------
+        viewpoint: numeric
+            Non-negative number. Can be interpreted as the degree of
+            ignorance towards rare species, where 0 treats rare species
+            the same as frequent species, and infinity considers only the
+            most frequent species.
+        """
+        return self.__metacommunity_measure(viewpoint, self.subcommunity_gamma)
+
+    def normalized_metacommunity_alpha(self, viewpoint):
+        """Calculates alpha class diversity of metacommunity.
+
+        Corresponds roughly to the average diversity of subcommunities
+        in isolation.
+
+        Parameters
+        ----------
+        viewpoint: numeric
+            Non-negative number. Can be interpreted as the degree of
+            ignorance towards rare species, where 0 treats rare species
+            the same as frequent species, and infinity considers only the
+            most frequent species.
+        """
+        return self.__metacommunity_measure(
+            viewpoint, self.normalized_subcommunity_alpha
+        )
+
+    def normalized_metacommunity_rho(self, viewpoint):
+        """Calculates rho class diversitiy of metacommunity.
+
+        Corresponds roughly to the average representativeness of
+        subcommunities.
+
+        Parameters
+        ----------
+        viewpoint: numeric
+            Non-negative number. Can be interpreted as the degree of
+            ignorance towards rare species, where 0 treats rare species
+            the same as frequent species, and infinity considers only the
+            most frequent species.
+        """
+        return self.__metacommunity_measure(viewpoint, self.normalized_subcommunity_rho)
+
+    def normalized_metacommunity_beta(self, viewpoint):
+        """Calculates beta class diversity of metacommunity.
+
+        Corresponds roughly to the effective number of distinct
+        subcommunities.
+
+        Parameters
+        ----------
+        viewpoint: numeric
+            Non-negative number. Can be interpreted as the degree of
+            ignorance towards rare species, where 0 treats rare species
+            the same as frequent species, and infinity considers only the
+            most frequent species.
+        """
+        return self.__metacommunity_measure(
+            viewpoint, self.normalized_subcommunity_beta
+        )
+
+    def __subcommunity_measure(self, viewpoint, numerator, denominator):
+        """Calculates subcommunity diversity measures."""
+        similarities = divide(
+            numerator, denominator, out=zeros(denominator.shape), where=denominator != 0
+        )
+        return power_mean(
+            1 - viewpoint,
+            self.abundance.normalized_subcommunity_abundance,
+            similarities,
+        )
+
+    def __metacommunity_measure(self, viewpoint, subcommunity_function):
+        """Calculates metcommunity diversity measures."""
+        subcommunity_measure = subcommunity_function(viewpoint)
+        return power_mean(
+            1 - viewpoint,
+            self.abundance.subcommunity_normalizing_constants,
+            subcommunity_measure,
+        )
+
+    def subcommunities_to_dataframe(self, viewpoint):
+        """Table containing all subcommunity diversity values.
+
+        Parameters
+        ----------
+        viewpoint: numeric
+            Non-negative number. Can be interpreted as the degree of
+            ignorance towards rare species, where 0 treats rare species
+            the same as frequent species, and infinity considers only the
+            most frequent species.
+        """
+        return DataFrame(
+            {
+                "community": self.abundance.subcommunity_order,
+                "viewpoint": viewpoint,
+                "alpha": self.subcommunity_alpha(viewpoint),
+                "rho": self.subcommunity_rho(viewpoint),
+                "beta": self.subcommunity_beta(viewpoint),
+                "gamma": self.subcommunity_gamma(viewpoint),
+                "normalized_alpha": self.normalized_subcommunity_alpha(viewpoint),
+                "normalized_rho": self.normalized_subcommunity_rho(viewpoint),
+                "normalized_beta": self.normalized_subcommunity_beta(viewpoint),
+            }
+        )
+
+    def metacommunity_to_dataframe(self, viewpoint):
+        """Table containing all metacommunity diversity values.
+
+        Parameters
+        ----------
+        viewpoint: numeric
+            Non-negative number. Can be interpreted as the degree of
+            ignorance towards rare species, where 0 treats rare species
+            the same as frequent species, and infinity considers only the
+            most frequent species."""
+        return DataFrame(
+            {
+                "community": "metacommunity",
+                "viewpoint": viewpoint,
+                "alpha": self.metacommunity_alpha(viewpoint),
+                "rho": self.metacommunity_rho(viewpoint),
+                "beta": self.metacommunity_beta(viewpoint),
+                "gamma": self.metacommunity_gamma(viewpoint),
+                "normalized_alpha": self.normalized_metacommunity_alpha(viewpoint),
+                "normalized_rho": self.normalized_metacommunity_rho(viewpoint),
+                "normalized_beta": self.normalized_metacommunity_beta(viewpoint),
+            },
+            index=[0],
+        )
+
+
+def make_metacommunity(
+    counts,
+    similarity_matrix=None,
+    similarity_matrix_filepath=None,
+    similarity_function=None,
+    features=None,
+    species_order=None,
+):
+    """Builds a Metacommunity object from specified parameters.
+
+    Valid parameter combinations and the corresponding Similarity
+    implementations used are:
+        - similarity_matrix, species_order -> SimilarityFromMemory
+        - similarity_matrix_filepath -> SimilarityFromFile
+        - species_order, similarity_function, features -> SimilarityFromFunction
+
+    Parameters
+    ----------
+    counts: numpy.ndarray or pandas.DataFrame
+        See diversity.metacommunity.Abundance. If the object is a
+        pandas.DataFrame, its to_numpy method should return the expected
+        numpy.ndarray.
+    similarity_matrix: numpy.ndarray
+        See diversity.metacommunity.SimilarityFromMemory.
+    similarity_matrix_filepath: str
+        See diversity.metacommunity.SimilarityFromFile
+    similarity_function: Callable
+        See diversity.metacommunity.SimilarityFromFunction.
+    features: numpy.ndarray
+        See diversity.metacommunity.SimilarityFromFunction.
+    species_order: Iterable
+        See diversity.metacommunity.SimilarityFromMemory, or
+        diversity.metacommunity.SimilarityFromFunction.
+
+    Returns
+    -------
+    A diversity.metacommunity.Metacommunity object build according to
+    parameter specification.
+    """
+
+    if_dataframe_to_numpy = lambda x: x.to_numpy() if isinstance(x, DataFrame) else x
+    similarity_matrix = if_dataframe_to_numpy(similarity_matrix)
+    features = if_dataframe_to_numpy(features)
+    similarity = make_similarity(
+        similarity_matrix,
+        species_order,
+        similarity_matrix_filepath,
+        similarity_function,
+        features,
+    )
+    counts = if_dataframe_to_numpy(counts)
+    abundance = Abundance(counts, similarity.species_order)
+    return Metacommunity(similarity, abundance)
diff --git a/src/diversity/parameters.py b/src/diversity/parameters.py
new file mode 100644
index 00000000..dcf05f36
--- /dev/null
+++ b/src/diversity/parameters.py
@@ -0,0 +1,118 @@
+"""Module for diversity's argument parser configuration.
+
+Classes
+-------
+ValidateViewpoint
+    Validator for -viewpoint parameter.
+
+Functions
+---------
+configure_arguments
+    Creates argument parser for .
+"""
+from argparse import Action, ArgumentParser
+from sys import stdout, stdin
+from warnings import warn
+
+from numpy import inf
+
+from diversity.utilities import ArgumentWarning
+
+
+class ValidateViewpoint(Action):
+    """Validator for -viewpoint parameter."""
+
+    def __call__(self, parser, args, values, option_string=None):
+        """Validates -viewpoint parameter.
+
+        Warns if arguments larger than 100 are passed, reminding that
+        they are treated as infinity in the diversity calculation.
+        """
+        if any([viewpoint > 100 and viewpoint != inf for viewpoint in values]):
+            warn(
+                "viewpoints > 100.0 defaults to the analytical formula"
+                " for viewpoint = infinity.",
+                category=ArgumentWarning,
+            )
+        setattr(args, self.dest, values)
+
+
+def configure_arguments():
+    """Creates argument parser.
+
+    Returns
+    -------
+    argparse.ArgumentParser configured to handle command-line arguments
+    for Metacommunity.
+    """
+    parser = ArgumentParser()
+    parser.add_argument(
+        "-c",
+        "--count_column",
+        help="Header of the column containing the counts.",
+        default="count",
+    )
+    parser.add_argument(
+        "-i",
+        "--input_filepath",
+        help=(
+            "A csv or tsv file containing one column for subcommunity"
+            " names (with header 'subcommunity', or specified via -u),"
+            " one column for species names (with header 'species', or"
+            " specified via -p), and one column of counts, which"
+            " contains the number of individuals of the row's species"
+            " in the row's subcommunity (with header 'count', or"
+            " specified via -c)."
+        ),
+    )
+    parser.add_argument(
+        "-l",
+        "--log_level",
+        help=(
+            "Logging verbosity level. Must be one of DEBUG, INFO,"
+            " WARNING, ERROR, CRITICAL (listed in decreasing"
+            " verbosity)."
+        ),
+        default="INFO",
+    )
+    parser.add_argument(
+        "-o",
+        "--output_filepath",
+        default=stdout,
+        help="A filepath to where the program's output will be saved",
+    )
+    parser.add_argument(
+        "-p",
+        "--species_column",
+        help="Header of the column containing the species names.",
+        default="species",
+    )
+    parser.add_argument(
+        "-s",
+        "--similarity_matrix_filepath",
+        help=(
+            "The filepath to a csv or tsv file containing a similarity"
+            " for the species in the input file. The file must have a"
+            " header listing the species names corresponding to each"
+            " column, and column and row ordering must be the same."
+        ),
+    )
+    parser.add_argument(
+        "-u",
+        "--subcommunity_column",
+        help="Header for the column containing the subcommunity names.",
+        default="subcommunity",
+    )
+    parser.add_argument(
+        "-v",
+        "--viewpoint",
+        nargs="+",
+        type=float,
+        help=(
+            "A list of viewpoint parameters. Any non-negative number"
+            " (including inf) is valid, but viewpoint parameters"
+            " greater than 100 are treated like inf."
+        ),
+        action=ValidateViewpoint,
+    )
+    return parser
diff --git a/src/diversity/utilities.py b/src/diversity/utilities.py
new file mode 100644
index 00000000..7b5b322b
--- /dev/null
+++ b/src/diversity/utilities.py
@@ -0,0 +1,195 @@
+"""Miscellaneous helper module for the metacommunity package.
+
+Functions
+---------
+get_file_delimiter
+    Determines delimiter in datafile from file extension.
+power_mean
+    Calculates weighted power means.
+unique_correspondence
+    Describes correspondence between a sequence of items and the set of
+    their unique elements.
+
+Exceptions
+----------
+MetacommunityError
+    Base class for all custom metacommunity exceptions.
+InvalidArgumentError
+    Raised when invalid argument is passed to a function.
+ArgumentWarning
+    Used for warnings of problematic argument choices.
+"""
+from pathlib import Path
+from warnings import warn
+
+from numpy import (
+    array,
+    isclose,
+    prod,
+    amin,
+    amax,
+    multiply,
+    sum as numpy_sum,
+    inf,
+    power,
+    unique,
+)
+
+
+class MetacommunityError(Exception):
+    """Base class for all custom metacommunity exceptions."""
+
+    pass
+
+
+class InvalidArgumentError(MetacommunityError):
+    """Raised when a function receives an invalid argument."""
+
+    pass
+
+
+class ArgumentWarning(Warning):
+    """Used for warnings related to problematic argument choices."""
+
+    pass
+
+
+def get_file_delimiter(filepath):
+    """Determines delimiter of file from filepath ending.
+
+    Parameters
+    ----------
+    filepath: str
+        The filepath whose delimiter is to be determine.
+
+    Returns
+    -------
+    "," if the file has the csv extension and "\t" otherwise.
+    """
+    suffix = Path(filepath).suffix
+    if suffix == ".csv":
+        return ","
+    elif suffix == ".tsv":
+        return "\t"
+    else:
+        warn(
+            f"File extension for {filepath} not recognized. Assuming"
+            " tab-delimited file.",
+            category=ArgumentWarning,
+        )
+        return "\t"
+
+
+def __validate_power_mean_args(weights, items, atol, weight_is_nonzero):
+    """Validates arguments for power_mean.
+
+    Parameters
+    ----------
+    weights, items
+        Same as for power_mean.
+    weight_is_nonzero
+        Boolean array of same shape as weights, indicating those whose
+        absolute value meets or exceeds atol.
+    """
+    if len(weights.shape) > 2:
+        raise InvalidArgumentError(
+            f"Invalid weights shape for power_mean: {weights.shape}."
+        )
+    if weights.shape != items.shape:
+        raise InvalidArgumentError(
+            f"Shape of weights ({weights.shape}) must be the same as"
+            f" shape fo items ({items.shape})."
+        )
+    all_0_column = None
+    if len(weights.shape) == 1:
+        all_0_column = [not (weight_is_nonzero).any()]
+    elif len(weights.shape) == 2:
+        all_0_column = [
+            not (weight_is_nonzero[:, col]).any() for col in range(weights.shape[1])
+        ]
+    if any(all_0_column):
+        raise InvalidArgumentError(
+            "power_mean expects at least one zero weight. Weights are"
+            " considered 0, if absolute value does not meet or exceed"
+            f" configurable minimum threshold: {atol:.2e}."
+        )
+
+
+def power_mean(order, weights, items, atol=1e-9):
+    """Calculates weighted power means.
+
+    Parameters
+    ----------
+    order: numeric
+        Exponent used for the power mean.
+    weights: numpy.ndarray
+        The weights corresponding to items. Must be 1-d or 2-d. If 2-d,
+        each column must contain at least one wheight whose absolute
+        value meets or exceeds atol.
+    items: numpy.ndarray
+        The elements for which the weighted power mean is computed. Must
+        have same shape as weights.
+    atol: float
+        Threshold below which weights are considered to be 0.
+
+    Returns
+    -------
+    A numpy.ndarray of the power means of items along axis 0 using order
+    as exponent, weighing by weights. The array shape is the same as
+    that of weights, or items except with the 0-axis removed. In the
+    case of 1-d weights and items, the result has shape (1,). When order
+    is close to 0 (absolute value less than atol), less than -100, or
+    greater than 100 analytical formulas for the limits at 0, -infinity,
+    or infinity are used respectively. An exception is raised if all
+    weights or a column of weights (in the 2-d case) are all too close
+    to 0.
+    """
+    weight_is_nonzero = abs(weights) >= atol
+    __validate_power_mean_args(weights, items, atol, weight_is_nonzero)
+    if isclose(order, 0):
+        return prod(
+            power(items, weights, where=weight_is_nonzero),
+            axis=0,
+            where=weight_is_nonzero,
+        )
+    elif order < -100:
+        return amin(items, axis=0, where=weight_is_nonzero, initial=inf)
+    elif order > 100:
+        return amax(items, axis=0, where=weight_is_nonzero, initial=-inf)
+    items_power = power(items, order, where=weight_is_nonzero)
+    items_product = multiply(items_power, weights, where=weight_is_nonzero)
+    items_sum = numpy_sum(items_product, axis=0, where=weight_is_nonzero)
+    return power(items_sum, 1 / order)
+
+
+def unique_correspondence(items, ordered_unique_items=None):
+    """Returns uniqued items and a mapping from items to uniqued items.
+
+    Parameters
+    ----------
+    items: numpy.ndarray
+        Array of items to unique and/or obtain mapping for.
+    ordered_unique_items: Iterable
+        Unique items in desired order. If None, ordering will be
+        established according to numpy.unique.
+
+    Returns
+    -------
+    A tuple with coordinates:
+    0 - numpy.ndarray
+        The ordered unique items.
+    1 - numpy.ndarray
+        The position in the unique items iterable for each item in
+        items.
+    """
+    if ordered_unique_items is None:
+        ordered_unique_items_, item_positions = unique(items, return_inverse=True)
+    else:
+        ordered_unique_items_ = array(ordered_unique_items)
+        item_to_position = {
+            str(item): pos for pos, item in enumerate(ordered_unique_items_)
+        }
+        if len(item_to_position) != len(ordered_unique_items):
+            raise InvalidArgumentError(f"Expected ordered_unique_items to be uniqued.")
+        item_positions = array([item_to_position[str(item)] for item in items])
+    return (ordered_unique_items_, item_positions)
diff --git a/tox.ini b/tox.ini
new file mode 100644
index 00000000..df8b0d2a
--- /dev/null
+++ b/tox.ini
@@ -0,0 +1,17 @@
+[tox]
+envlist = virtual_env
+isolated_build = True
+
+[gh-actions]
+python =
+	3.8: py38
+	3.9: py39
+	3.10: py310
+
+[testenv]
+deps = 
+	pytest
+       	numpy
+       	pandas
+commands = 
+	pytest {posargs:tests}


Test patch:
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
new file mode 100644
index 00000000..742d317c
--- /dev/null
+++ b/.github/workflows/tests.yml
@@ -0,0 +1,28 @@
+name: Tests
+
+on:
+  - push
+    branches: [ main ]
+  - pull_request
+    branches: [ main ]
+
+jobs:
+  test:
+    runs-on: ${{ matrix.os}}
+    strategy:
+      matrix:
+        os: [macos-latest, ubuntu-latest, windows-latest]
+        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
+        pip install tox tox-gh-actions
+    - name: Test with tox
+      run: tox
diff --git a/test_diversity.py b/test_diversity.py
deleted file mode 100644
index b3ad9cbe..00000000
--- a/test_diversity.py
+++ /dev/null
@@ -1,55 +0,0 @@
-import pytest
-import diversity
-from numpy import isclose
-from run import process_input_file
-# FIXME not implemented: used to compare dataframes (for integration test)
-from pandas.testing import assert_frame_equal
-
-
-TEST_QS = [0, 0.5, 1, 1.5, 2, 3, 4, 100, 1000]
-TEST_QDS = 3
-
-
-@pytest.fixture
-def mock_input_file(tmp_path):
-    mock_csv_data = [
-        "name_1,1,0,AAA",
-        "name_2,1,0,BBB",
-        "name_3,1,0,CCC",
-    ]
-    datafile = tmp_path / "input.csv"
-    datafile.write_text("\n".join(mock_csv_data))
-    return str(datafile)
-
-
-@pytest.fixture
-def mock_matrix_file(tmp_path):
-    mock_csv_data = [
-        "1.0,0.0,0.0",
-        "0.0,1.0,0.0",
-        "0.0,0.0,1.0",
-    ]
-    datafile = tmp_path / "matrix.csv"
-    datafile.write_text("\n".join(mock_csv_data))
-    return str(datafile)
-
-
-@pytest.fixture
-def mock_empty_matrix_file(tmp_path):
-    datafile = tmp_path / "empty_matrix.csv"
-    return str(datafile)
-
-
-def test_alpha_from_file(mock_input_file, mock_matrix_file):
-    df = process_input_file(mock_input_file)
-    for q in TEST_QS:
-        qDs = diversity.alpha(df, q, z_filepath=mock_matrix_file)
-        assert isclose(TEST_QDS, qDs)
-
-
-def test_alpha_from_similarity_fn(mock_input_file, mock_empty_matrix_file):
-    df = process_input_file(mock_input_file)
-    for q in TEST_QS:
-        qDs = diversity.alpha(
-            df, q, z_filepath=mock_empty_matrix_file, similarity_fn=diversity.sequence_similarity)
-        assert isclose(TEST_QDS, qDs)
diff --git a/tests/__init__.py b/tests/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/main_test.py b/tests/main_test.py
new file mode 100644
index 00000000..598e917b
--- /dev/null
+++ b/tests/main_test.py
@@ -0,0 +1,165 @@
+"""Tests for diversity.__main__."""
+from argparse import Namespace
+
+from numpy import inf
+from pytest import mark
+
+# from pandas import read_tsv, DataFrame
+# from pandas.testing import assert_frame_equal
+
+# from diversity.metacommunity import Abundance, Metacommunity, Similarity
+from diversity.__main__ import main
+
+MAIN_TEST_CASES = [
+    {
+        "description": "disjoint communities; uniform counts; uniform inter-community similarities; viewpoint 0.",
+        "args": Namespace(
+            input_filepath="counts.tsv",
+            output_filepath="diversities.tsv",
+            similarity_matrix_filepath="similarities.tsv",
+            viewpoint=[0],
+            log_level="WARNING",
+            subcommunity_column="subcommunity",
+            species_column="species",
+            count_column="count",
+        ),
+        "input_filecontents": (
+            "subcommunity\tspecies\tcount\n"
+            "subcommunity_1\tspecies_1\t1\n"
+            "subcommunity_1\tspecies_2\t1\n"
+            "subcommunity_1\tspecies_3\t1\n"
+            "subcommunity_2\tspecies_4\t1\n"
+            "subcommunity_2\tspecies_5\t1\n"
+            "subcommunity_2\tspecies_6\t1\n"
+        ),
+        "similarities_filecontents": (
+            "species_1\tspecies_2\tspecies_3\tspecies_4\tspecies_5\tspecies_6\n"
+            "1.0\t0.5\t0.5\t0.7\t0.7\t0.7\n"
+            "0.5\t1.0\t0.5\t0.7\t0.7\t0.7\n"
+            "0.5\t0.5\t1.0\t0.7\t0.7\t0.7\n"
+            "0.7\t0.7\t0.7\t1.0\t0.5\t0.5\n"
+            "0.7\t0.7\t0.7\t0.5\t1.0\t0.5\n"
+            "0.7\t0.7\t0.7\t0.5\t0.5\t1.0\n"
+        ),
+        "output_filecontents": (
+            "community\tviewpoint\talpha\trho\tbeta\tgamma\tnormalized_alpha\tnormalized_rho\tnormalized_beta\n"
+            "subcommunity_1\t0.00\t3.0000\t2.0500\t0.4878\t1.4634\t1.5000\t1.0250\t0.9756\n"
+            "subcommunity_2\t0.00\t3.0000\t2.0500\t0.4878\t1.4634\t1.5000\t1.0250\t0.9756\n"
+            "metacommunity\t0.00\t3.0000\t2.0500\t0.4878\t1.4634\t1.5000\t1.0250\t0.9756\n"
+        ),
+    },
+    {
+        "description": "overlapping communities; non-uniform counts; non-uniform inter-community similarities; viewpoint 2.",
+        "args": Namespace(
+            input_filepath="foo_counts.tsv",
+            output_filepath="bar_counts.tsv",
+            similarity_matrix_filepath="baz_similarities.tsv",
+            viewpoint=[2, 101, 102, inf],
+            log_level="WARNING",
+            subcommunity_column="subcommunity",
+            species_column="species",
+            count_column="count",
+        ),
+        "input_filecontents": (
+            "species\tsubcommunity\tcount\n"
+            "species_1\tsubcommunity_1\t2\n"
+            "species_2\tsubcommunity_1\t3\n"
+            "species_1\tsubcommunity_2\t5\n"
+            "species_3\tsubcommunity_2\t1\n"
+        ),
+        "similarities_filecontents": (
+            "species_1\tspecies_2\tspecies_3\n"
+            "1.0\t0.5\t0.1\n"
+            "0.5\t1.0\t0.2\n"
+            "0.1\t0.2\t1.0\n"
+        ),
+        "output_filecontents": (
+            "community\tviewpoint\talpha\trho\tbeta\tgamma\tnormalized_alpha\tnormalized_rho\tnormalized_beta\n"
+            "subcommunity_1\t2.00\t2.8947\t1.9194\t0.5210\t1.4745\t1.3158\t0.8724\t1.1462\n"
+            "subcommunity_2\t2.00\t2.4444\t1.6587\t0.6029\t1.4570\t1.3333\t0.9047\t1.1053\n"
+            "metacommunity\t2.00\t2.6304\t1.7678\t0.5627\t1.4649\t1.3253\t0.8898\t1.1235\n"
+            "subcommunity_1\t101.00\t2.7641\t1.6836\t0.5940\t1.2908\t1.2564\t0.7653\t1.3067\n"  # 2.76408362, 1.68357822, 0.593973, 1.29084362, 1.25640164, 0.76526283, 1.30674059
+            "subcommunity_2\t101.00\t2.1608\t1.5610\t0.6406\t1.2814\t1.1786\t0.8515\t1.1744\n"  # 2.16079876, 1.56104879, 0.64059497, 1.28140391, 1.1786175, 0.85148116, 1.17442411
+            "metacommunity\t101.00\t2.1739\t1.5705\t0.5987\t1.2849\t1.1858\t0.7713\t1.1816\n"  # 2.1739359070171167, 1.5705327595028187, 0.5986709767832536, 1.2848640573127161, 1.185766665953564, 0.7713202342679308, 1.1815641039291136
+            "subcommunity_1\t102.00\t2.7500\t1.6750\t0.5970\t1.2791\t1.2500\t0.7614\t1.3134\n"  # 2.75, 1.675, 0.59701493, 1.27906977, 1.25, 0.76136364, 1.31343284
+            "subcommunity_2\t102.00\t2.1569\t1.5333\t0.6522\t1.2791\t1.1765\t0.8364\t1.1957\n"  # 2.15686275, 1.53333333, 0.65217391, 1.27906977, 1.17647059, 0.83636364, 1.19565217
+            "metacommunity\t102.00\t2.1569\t1.5333\t0.5970\t1.2791\t1.1765\t0.7614\t1.1957\n"  # 2.156862745098039, 1.5333333333333337, 0.5970149253731344, 1.2790697674418605, 1.176470588235294, 0.7613636363636362, 1.1956521739130435
+            "subcommunity_1\tinf\t2.7500\t1.6750\t0.5970\t1.2791\t1.2500\t0.7614\t1.3134\n"  # 2.75, 1.675, 0.59701493, 1.27906977, 1.25, 0.76136364, 1.31343284
+            "subcommunity_2\tinf\t2.1569\t1.5333\t0.6522\t1.2791\t1.1765\t0.8364\t1.1957\n"  # 2.15686275, 1.53333333, 0.65217391, 1.27906977, 1.17647059, 0.83636364, 1.19565217
+            "metacommunity\tinf\t2.1569\t1.5333\t0.5970\t1.2791\t1.1765\t0.7614\t1.1957\n"  # 2.156862745098039, 1.5333333333333337, 0.5970149253731344, 1.2790697674418605, 1.176470588235294, 0.7613636363636362, 1.1956521739130435
+        ),
+    },
+    {
+        "description": "overlapping communities; non-uniform counts; non-uniform inter-community similarities; viewpoint 2.",
+        "args": Namespace(
+            input_filepath="foo_counts.tsv",
+            output_filepath="bar_counts.tsv",
+            similarity_matrix_filepath="baz_similarities.tsv",
+            viewpoint=[2, 101, 102, inf],
+            log_level="WARNING",
+            subcommunity_column="subcommunity_",
+            species_column="species_",
+            count_column="count_",
+        ),
+        "input_filecontents": (
+            "species\tsubcommunity_\tcount_\tspecies_\tsubcommunity\tcount\n"
+            "species_5\tsubcommunity_1\t2\tspecies_1\tsubcommunity_x\t100\n"
+            "species_9\tsubcommunity_1\t3\tspecies_2\tsubcommunity_y\t55\n"
+            "species_9\tsubcommunity_2\t5\tspecies_1\tsubcommunity_z\t12\n"
+            "species_0\tsubcommunity_2\t1\tspecies_3\tsubcommunity_z\t19\n"
+        ),
+        "similarities_filecontents": (
+            "species_1\tspecies_2\tspecies_3\n"
+            "1.0\t0.5\t0.1\n"
+            "0.5\t1.0\t0.2\n"
+            "0.1\t0.2\t1.0\n"
+        ),
+        "output_filecontents": (
+            "community\tviewpoint\talpha\trho\tbeta\tgamma\tnormalized_alpha\tnormalized_rho\tnormalized_beta\n"
+            "subcommunity_1\t2.00\t2.8947\t1.9194\t0.5210\t1.4745\t1.3158\t0.8724\t1.1462\n"
+            "subcommunity_2\t2.00\t2.4444\t1.6587\t0.6029\t1.4570\t1.3333\t0.9047\t1.1053\n"
+            "metacommunity\t2.00\t2.6304\t1.7678\t0.5627\t1.4649\t1.3253\t0.8898\t1.1235\n"
+            "subcommunity_1\t101.00\t2.7641\t1.6836\t0.5940\t1.2908\t1.2564\t0.7653\t1.3067\n"  # 2.76408362, 1.68357822, 0.593973, 1.29084362, 1.25640164, 0.76526283, 1.30674059
+            "subcommunity_2\t101.00\t2.1608\t1.5610\t0.6406\t1.2814\t1.1786\t0.8515\t1.1744\n"  # 2.16079876, 1.56104879, 0.64059497, 1.28140391, 1.1786175, 0.85148116, 1.17442411
+            "metacommunity\t101.00\t2.1739\t1.5705\t0.5987\t1.2849\t1.1858\t0.7713\t1.1816\n"  # 2.1739359070171167, 1.5705327595028187, 0.5986709767832536, 1.2848640573127161, 1.185766665953564, 0.7713202342679308, 1.1815641039291136
+            "subcommunity_1\t102.00\t2.7500\t1.6750\t0.5970\t1.2791\t1.2500\t0.7614\t1.3134\n"  # 2.75, 1.675, 0.59701493, 1.27906977, 1.25, 0.76136364, 1.31343284
+            "subcommunity_2\t102.00\t2.1569\t1.5333\t0.6522\t1.2791\t1.1765\t0.8364\t1.1957\n"  # 2.15686275, 1.53333333, 0.65217391, 1.27906977, 1.17647059, 0.83636364, 1.19565217
+            "metacommunity\t102.00\t2.1569\t1.5333\t0.5970\t1.2791\t1.1765\t0.7614\t1.1957\n"  # 2.156862745098039, 1.5333333333333337, 0.5970149253731344, 1.2790697674418605, 1.176470588235294, 0.7613636363636362, 1.1956521739130435
+            "subcommunity_1\tinf\t2.7500\t1.6750\t0.5970\t1.2791\t1.2500\t0.7614\t1.3134\n"  # 2.75, 1.675, 0.59701493, 1.27906977, 1.25, 0.76136364, 1.31343284
+            "subcommunity_2\tinf\t2.1569\t1.5333\t0.6522\t1.2791\t1.1765\t0.8364\t1.1957\n"  # 2.15686275, 1.53333333, 0.65217391, 1.27906977, 1.17647059, 0.83636364, 1.19565217
+            "metacommunity\tinf\t2.1569\t1.5333\t0.5970\t1.2791\t1.1765\t0.7614\t1.1957\n"  # 2.156862745098039, 1.5333333333333337, 0.5970149253731344, 1.2790697674418605, 1.176470588235294, 0.7613636363636362, 1.1956521739130435
+        ),
+    },
+]
+
+
+class TestMain:
+    """Tests __main__.main."""
+
+    def write_file(self, path, contents):
+        """Writes contents into file at path."""
+        with open(path, "w") as file:
+            file.write(contents)
+
+    @mark.parametrize("test_case", MAIN_TEST_CASES)
+    def test_main(self, test_case, tmp_path):
+        """Tests __main__.main."""
+        test_case["args"].input_filepath = tmp_path / test_case["args"].input_filepath
+        test_case["args"].similarity_matrix_filepath = (
+            tmp_path / test_case["args"].similarity_matrix_filepath
+        )
+        test_case["args"].output_filepath = tmp_path / test_case["args"].output_filepath
+
+        self.write_file(
+            test_case["args"].input_filepath,
+            test_case["input_filecontents"],
+        )
+        self.write_file(
+            test_case["args"].similarity_matrix_filepath,
+            test_case["similarities_filecontents"],
+        )
+
+        main(test_case["args"])
+        with open(test_case["args"].output_filepath, "r") as file:
+            output_filecontents = file.read()
+        assert output_filecontents == test_case["output_filecontents"]
diff --git a/tests/metacommunity_test.py b/tests/metacommunity_test.py
new file mode 100644
index 00000000..6f9cbea3
--- /dev/null
+++ b/tests/metacommunity_test.py
@@ -0,0 +1,1491 @@
+"""Tests for diversity.metacommunity."""
+from copy import deepcopy
+from itertools import product
+from warnings import filterwarnings, resetwarnings
+
+from numpy import allclose, array, empty, float64, isclose, ndarray, unique
+from pandas import DataFrame
+from pytest import mark
+
+from diversity.metacommunity import (
+    Abundance,
+    make_similarity,
+    make_metacommunity,
+    Metacommunity,
+    SimilarityFromFile,
+    SimilarityFromFunction,
+    SimilarityFromMemory,
+)
+from diversity.utilities import ArgumentWarning, unique_correspondence
+
+ABUNDANCE_TEST_CASES = [
+    ####################################################################
+    {
+        "description": "default parameters; 2 communities; both contain exclusive species",
+        "counts": array(
+            [
+                ["community_1", "species_1", "2"],
+                ["community_1", "species_2", "3"],
+                ["community_2", "species_1", "4"],
+                ["community_2", "species_3", "1"],
+            ]
+        ),
+        "species_order": None,
+        "subcommunity_order": None,
+        "subcommunity_column": 0,
+        "species_column": 1,
+        "count_column": 2,
+        "expected_species_order": unique_correspondence(
+            array(["species_1", "species_2", "species_1", "species_3"])
+        )[0],
+        "expected_subcommunity_order": unique_correspondence(
+            array(["community_1", "community_1", "community_2", "community_2"])
+        )[0],
+        "subcommunity_species_to_subcommunity_abundance": {
+            ("community_1", "species_1"): 2 / 10,
+            ("community_1", "species_2"): 3 / 10,
+            ("community_1", "species_3"): 0 / 10,
+            ("community_2", "species_1"): 4 / 10,
+            ("community_2", "species_2"): 0 / 10,
+            ("community_2", "species_3"): 1 / 10,
+        },
+        "species_to_metacommunity_abundance": {
+            "species_1": 6 / 10,
+            "species_2": 3 / 10,
+            "species_3": 1 / 10,
+        },
+        "subcommunity_to_normalizing_constants": {
+            "community_1": 5 / 10,
+            "community_2": 5 / 10,
+        },
+        "subcommunity_species_to_normalized_subcommunity_abundance": {
+            ("community_1", "species_1"): 2 / 5,
+            ("community_1", "species_2"): 3 / 5,
+            ("community_1", "species_3"): 0 / 5,
+            ("community_2", "species_1"): 4 / 5,
+            ("community_2", "species_2"): 0 / 5,
+            ("community_2", "species_3"): 1 / 5,
+        },
+    },
+    ####################################################################
+    {
+        "description": "default parameters; 2 communities; one contains exclusive species",
+        "counts": array(
+            [
+                ["community_1", "species_1", "2"],
+                ["community_1", "species_2", "3"],
+                ["community_2", "species_1", "4"],
+                ["community_2", "species_3", "1"],
+                ["community_1", "species_3", "5"],
+            ]
+        ),
+        "species_order": None,
+        "subcommunity_order": None,
+        "subcommunity_column": 0,
+        "species_column": 1,
+        "count_column": 2,
+        "expected_species_order": unique_correspondence(
+            array(["species_1", "species_2", "species_1", "species_3", "species_3"])
+        )[0],
+        "expected_subcommunity_order": unique_correspondence(
+            array(
+                [
+                    "community_1",
+                    "community_1",
+                    "community_2",
+                    "community_2",
+                    "community_1",
+                ]
+            )
+        )[0],
+        "subcommunity_species_to_subcommunity_abundance": {
+            ("community_1", "species_1"): 2 / 15,
+            ("community_1", "species_2"): 3 / 15,
+            ("community_1", "species_3"): 5 / 15,
+            ("community_2", "species_1"): 4 / 15,
+            ("community_2", "species_2"): 0 / 15,
+            ("community_2", "species_3"): 1 / 15,
+        },
+        "species_to_metacommunity_abundance": {
+            "species_1": 6 / 15,
+            "species_2": 3 / 15,
+            "species_3": 6 / 15,
+        },
+        "subcommunity_to_normalizing_constants": {
+            "community_1": 10 / 15,
+            "community_2": 5 / 15,
+        },
+        "subcommunity_species_to_normalized_subcommunity_abundance": {
+            ("community_1", "species_1"): 2 / 10,
+            ("community_1", "species_2"): 3 / 10,
+            ("community_1", "species_3"): 5 / 10,
+            ("community_2", "species_1"): 4 / 5,
+            ("community_2", "species_2"): 0 / 5,
+            ("community_2", "species_3"): 1 / 5,
+        },
+    },
+    ####################################################################
+    {
+        "description": "default parameters; 2 communities; neither contain exclusive species",
+        "counts": array(
+            [
+                ["community_1", "species_1", "2"],
+                ["community_1", "species_2", "3"],
+                ["community_2", "species_1", "4"],
+                ["community_2", "species_3", "1"],
+                ["community_1", "species_3", "5"],
+                ["community_2", "species_2", "1"],
+            ]
+        ),
+        "species_order": None,
+        "subcommunity_order": None,
+        "subcommunity_column": 0,
+        "species_column": 1,
+        "count_column": 2,
+        "expected_species_order": unique_correspondence(
+            array(
+                [
+                    "species_1",
+                    "species_2",
+                    "species_1",
+                    "species_3",
+                    "species_3",
+                    "species_2",
+                ]
+            )
+        )[0],
+        "expected_subcommunity_order": unique_correspondence(
+            array(
+                [
+                    "community_1",
+                    "community_1",
+                    "community_2",
+                    "community_2",
+                    "community_1",
+                    "community_2",
+                ]
+            )
+        )[0],
+        "subcommunity_species_to_subcommunity_abundance": {
+            ("community_1", "species_1"): 2 / 16,
+            ("community_1", "species_2"): 3 / 16,
+            ("community_1", "species_3"): 5 / 16,
+            ("community_2", "species_1"): 4 / 16,
+            ("community_2", "species_2"): 1 / 16,
+            ("community_2", "species_3"): 1 / 16,
+        },
+        "species_to_metacommunity_abundance": {
+            "species_1": 6 / 16,
+            "species_2": 4 / 16,
+            "species_3": 6 / 16,
+        },
+        "subcommunity_to_normalizing_constants": {
+            "community_1": 10 / 16,
+            "community_2": 6 / 16,
+        },
+        "subcommunity_species_to_normalized_subcommunity_abundance": {
+            ("community_1", "species_1"): 2 / 10,
+            ("community_1", "species_2"): 3 / 10,
+            ("community_1", "species_3"): 5 / 10,
+            ("community_2", "species_1"): 4 / 6,
+            ("community_2", "species_2"): 1 / 6,
+            ("community_2", "species_3"): 1 / 6,
+        },
+    },
+    ####################################################################
+    {
+        "description": "default parameters; 2 mutually exclusive communities",
+        "counts": array(
+            [
+                ["community_1", "species_1", "2"],
+                ["community_1", "species_2", "3"],
+                ["community_2", "species_4", "4"],
+                ["community_2", "species_3", "1"],
+            ]
+        ),
+        "species_order": None,
+        "subcommunity_order": None,
+        "subcommunity_column": 0,
+        "species_column": 1,
+        "count_column": 2,
+        "expected_species_order": unique_correspondence(
+            array(["species_1", "species_2", "species_4", "species_3"])
+        )[0],
+        "expected_subcommunity_order": unique_correspondence(
+            array(
+                [
+                    "community_1",
+                    "community_1",
+                    "community_2",
+                    "community_2",
+                ]
+            )
+        )[0],
+        "subcommunity_species_to_subcommunity_abundance": {
+            ("community_1", "species_1"): 2 / 10,
+            ("community_1", "species_2"): 3 / 10,
+            ("community_1", "species_3"): 0 / 10,
+            ("community_1", "species_4"): 0 / 10,
+            ("community_2", "species_1"): 0 / 10,
+            ("community_2", "species_2"): 0 / 10,
+            ("community_2", "species_3"): 1 / 10,
+            ("community_2", "species_4"): 4 / 10,
+        },
+        "species_to_metacommunity_abundance": {
+            "species_1": 2 / 10,
+            "species_2": 3 / 10,
+            "species_3": 1 / 10,
+            "species_4": 4 / 10,
+        },
+        "subcommunity_to_normalizing_constants": {
+            "community_1": 5 / 10,
+            "community_2": 5 / 10,
+        },
+        "subcommunity_species_to_normalized_subcommunity_abundance": {
+            ("community_1", "species_1"): 2 / 5,
+            ("community_1", "species_2"): 3 / 5,
+            ("community_1", "species_3"): 0 / 5,
+            ("community_1", "species_4"): 0 / 5,
+            ("community_2", "species_1"): 0 / 5,
+            ("community_2", "species_2"): 0 / 5,
+            ("community_2", "species_3"): 1 / 5,
+            ("community_2", "species_4"): 4 / 5,
+        },
+    },
+    ####################################################################
+    {
+        "description": "default parameters; 1 community",
+        "counts": array(
+            [
+                ["community_1", "species_1", "2"],
+                ["community_1", "species_3", "3"],
+                ["community_1", "species_2", "5"],
+            ]
+        ),
+        "species_order": None,
+        "subcommunity_order": None,
+        "subcommunity_column": 0,
+        "species_column": 1,
+        "count_column": 2,
+        "expected_species_order": unique_correspondence(
+            array(["species_1", "species_3", "species_2"])
+        )[0],
+        "expected_subcommunity_order": unique_correspondence(
+            array(
+                [
+                    "community_1",
+                    "community_1",
+                    "community_1",
+                ]
+            )
+        )[0],
+        "subcommunity_species_to_subcommunity_abundance": {
+            ("community_1", "species_1"): 2 / 10,
+            ("community_1", "species_2"): 5 / 10,
+            ("community_1", "species_3"): 3 / 10,
+        },
+        "species_to_metacommunity_abundance": {
+            "species_1": 2 / 10,
+            "species_2": 5 / 10,
+            "species_3": 3 / 10,
+        },
+        "subcommunity_to_normalizing_constants": {
+            "community_1": 10 / 10,
+        },
+        "subcommunity_species_to_normalized_subcommunity_abundance": {
+            ("community_1", "species_1"): 2 / 10,
+            ("community_1", "species_2"): 5 / 10,
+            ("community_1", "species_3"): 3 / 10,
+        },
+    },
+    ####################################################################
+    {
+        "description": "nondefault species_order; 2 communities; both contain exclusive species",
+        "counts": array(
+            [
+                ["community_1", "species_1", "2"],
+                ["community_1", "species_2", "3"],
+                ["community_2", "species_1", "4"],
+                ["community_2", "species_3", "1"],
+            ]
+        ),
+        "species_order": array(["species_2", "species_1", "species_3"]),
+        "subcommunity_order": None,
+        "subcommunity_column": 0,
+        "species_column": 1,
+        "count_column": 2,
+        "expected_species_order": array(["species_2", "species_1", "species_3"]),
+        "expected_subcommunity_order": unique_correspondence(
+            array(["community_1", "community_1", "community_2", "community_2"])
+        )[0],
+        "subcommunity_species_to_subcommunity_abundance": {
+            ("community_1", "species_1"): 2 / 10,
+            ("community_1", "species_2"): 3 / 10,
+            ("community_1", "species_3"): 0 / 10,
+            ("community_2", "species_1"): 4 / 10,
+            ("community_2", "species_2"): 0 / 10,
+            ("community_2", "species_3"): 1 / 10,
+        },
+        "species_to_metacommunity_abundance": {
+            "species_1": 6 / 10,
+            "species_2": 3 / 10,
+            "species_3": 1 / 10,
+        },
+        "subcommunity_to_normalizing_constants": {
+            "community_1": 5 / 10,
+            "community_2": 5 / 10,
+        },
+        "subcommunity_species_to_normalized_subcommunity_abundance": {
+            ("community_1", "species_1"): 2 / 5,
+            ("community_1", "species_2"): 3 / 5,
+            ("community_1", "species_3"): 0 / 5,
+            ("community_2", "species_1"): 4 / 5,
+            ("community_2", "species_2"): 0 / 5,
+            ("community_2", "species_3"): 1 / 5,
+        },
+    },
+    ####################################################################
+    {
+        "description": "nondefault subcommunity_order; 3 communities; all contain exclusive species",
+        "counts": array(
+            [
+                ["community_1", "species_2", "3"],
+                ["community_2", "species_1", "4"],
+                ["community_1", "species_1", "2"],
+                ["community_2", "species_3", "1"],
+                ["community_3", "species_1", "1"],
+                ["community_3", "species_4", "1"],
+            ]
+        ),
+        "species_order": None,
+        "subcommunity_order": array(["community_2", "community_1", "community_3"]),
+        "subcommunity_column": 0,
+        "species_column": 1,
+        "count_column": 2,
+        "expected_species_order": unique_correspondence(
+            array(
+                [
+                    "species_1",
+                    "species_2",
+                    "species_1",
+                    "species_3",
+                    "species_1",
+                    "species_4",
+                ]
+            )
+        )[0],
+        "expected_subcommunity_order": array(
+            ["community_2", "community_1", "community_3"]
+        ),
+        "subcommunity_species_to_subcommunity_abundance": {
+            ("community_1", "species_1"): 2 / 12,
+            ("community_1", "species_2"): 3 / 12,
+            ("community_1", "species_3"): 0 / 12,
+            ("community_1", "species_4"): 0 / 12,
+            ("community_2", "species_1"): 4 / 12,
+            ("community_2", "species_2"): 0 / 12,
+            ("community_2", "species_3"): 1 / 12,
+            ("community_2", "species_4"): 0 / 12,
+            ("community_3", "species_1"): 1 / 12,
+            ("community_3", "species_2"): 0 / 12,
+            ("community_3", "species_3"): 0 / 12,
+            ("community_3", "species_4"): 1 / 12,
+        },
+        "species_to_metacommunity_abundance": {
+            "species_1": 7 / 12,
+            "species_2": 3 / 12,
+            "species_3": 1 / 12,
+            "species_4": 1 / 12,
+        },
+        "subcommunity_to_normalizing_constants": {
+            "community_1": 5 / 12,
+            "community_2": 5 / 12,
+            "community_3": 2 / 12,
+        },
+        "subcommunity_species_to_normalized_subcommunity_abundance": {
+            ("community_1", "species_1"): 2 / 5,
+            ("community_1", "species_2"): 3 / 5,
+            ("community_1", "species_3"): 0 / 5,
+            ("community_1", "species_4"): 0 / 5,
+            ("community_2", "species_1"): 4 / 5,
+            ("community_2", "species_2"): 0 / 5,
+            ("community_2", "species_3"): 1 / 5,
+            ("community_2", "species_4"): 0 / 5,
+            ("community_3", "species_1"): 1 / 2,
+            ("community_3", "species_2"): 0 / 2,
+            ("community_3", "species_3"): 0 / 2,
+            ("community_3", "species_4"): 1 / 2,
+        },
+    },
+    ####################################################################
+    {
+        "description": "permuted columns; 2 communities; both contain exclusive species",
+        "counts": array(
+            [
+                ["2", "community_1", "species_1"],
+                ["3", "community_1", "species_2"],
+                ["4", "community_2", "species_1"],
+                ["1", "community_2", "species_3"],
+            ]
+        ),
+        "species_order": None,
+        "subcommunity_order": None,
+        "subcommunity_column": 1,
+        "species_column": 2,
+        "count_column": 0,
+        "expected_species_order": unique_correspondence(
+            array(["species_1", "species_2", "species_1", "species_3"])
+        )[0],
+        "expected_subcommunity_order": unique_correspondence(
+            array(["community_1", "community_1", "community_2", "community_2"])
+        )[0],
+        "subcommunity_species_to_subcommunity_abundance": {
+            ("community_1", "species_1"): 2 / 10,
+            ("community_1", "species_2"): 3 / 10,
+            ("community_1", "species_3"): 0 / 10,
+            ("community_2", "species_1"): 4 / 10,
+            ("community_2", "species_2"): 0 / 10,
+            ("community_2", "species_3"): 1 / 10,
+        },
+        "species_to_metacommunity_abundance": {
+            "species_1": 6 / 10,
+            "species_2": 3 / 10,
+            "species_3": 1 / 10,
+        },
+        "subcommunity_to_normalizing_constants": {
+            "community_1": 5 / 10,
+            "community_2": 5 / 10,
+        },
+        "subcommunity_species_to_normalized_subcommunity_abundance": {
+            ("community_1", "species_1"): 2 / 5,
+            ("community_1", "species_2"): 3 / 5,
+            ("community_1", "species_3"): 0 / 5,
+            ("community_2", "species_1"): 4 / 5,
+            ("community_2", "species_2"): 0 / 5,
+            ("community_2", "species_3"): 1 / 5,
+        },
+    },
+]
+
+
+class TestAbundance:
+    """Tests diversity.metacommunity.Abundance."""
+
+    def make_abundance(self, test_case):
+        """Simple initializer without modifying test_case."""
+        test_case_ = deepcopy(test_case)
+        return Abundance(
+            counts=test_case_["counts"],
+            species_order=test_case_["species_order"],
+            subcommunity_order=test_case_["subcommunity_order"],
+            subcommunity_column=test_case_["subcommunity_column"],
+            species_column=test_case_["species_column"],
+            count_column=test_case_["count_column"],
+        )
+
+    @mark.parametrize("test_case", ABUNDANCE_TEST_CASES)
+    def test_init(self, test_case):
+        """Tests initializer."""
+        abundance = self.make_abundance(test_case)
+        assert (abundance.species_order == test_case["expected_species_order"]).all()
+        assert (
+            abundance.subcommunity_order == test_case["expected_subcommunity_order"]
+        ).all()
+
+    @mark.parametrize("test_case", ABUNDANCE_TEST_CASES)
+    def test_subcommunity_abundance(self, test_case):
+        """Tests .subcommunity_abundance."""
+        abundance = self.make_abundance(test_case)
+        subcommunity_abundance = abundance.subcommunity_abundance
+        for (i, species), (j, subcommunity) in product(
+            enumerate(test_case["expected_species_order"]),
+            enumerate(test_case["expected_subcommunity_order"]),
+        ):
+            assert isclose(
+                subcommunity_abundance[i, j],
+                test_case["subcommunity_species_to_subcommunity_abundance"][
+                    (subcommunity, species)
+                ],
+            ), (
+                f"\n(i,j): {(i,j)};"
+                f"\n(species,subcommunity): {(species,subcommunity)}."
+            )
+
+    @mark.parametrize("test_case", ABUNDANCE_TEST_CASES)
+    def test_metacommunity_abundance(self, test_case):
+        """Tests .metacommunity_abundance."""
+        abundance = self.make_abundance(test_case)
+        metacommunity_abundance = abundance.metacommunity_abundance
+        for i, species in enumerate(test_case["expected_species_order"]):
+            assert isclose(
+                metacommunity_abundance[i, 0],
+                test_case["species_to_metacommunity_abundance"][species],
+            ), f"\ni:{i};\nspecies: {species}."
+
+    @mark.parametrize("test_case", ABUNDANCE_TEST_CASES)
+    def test_subcommunity_normalizing_constants(self, test_case):
+        """Tests .subcommunity_normalizing_constants."""
+        abundance = self.make_abundance(test_case)
+        subcommunity_normalizing_constants = (
+            abundance.subcommunity_normalizing_constants
+        )
+        for j, subcommunity in enumerate(test_case["expected_subcommunity_order"]):
+            assert isclose(
+                subcommunity_normalizing_constants[j],
+                test_case["subcommunity_to_normalizing_constants"][subcommunity],
+            )
+
+    @mark.parametrize("test_case", ABUNDANCE_TEST_CASES)
+    def test_normalized_subcommunity_abundance(self, test_case):
+        """Tests .normalized_subcommunity_abundance."""
+        abundance = self.make_abundance(test_case)
+        normalized_subcommunity_abundance = abundance.normalized_subcommunity_abundance
+        for (i, species), (j, subcommunity) in product(
+            enumerate(test_case["expected_species_order"]),
+            enumerate(test_case["expected_subcommunity_order"]),
+        ):
+            assert isclose(
+                normalized_subcommunity_abundance[i, j],
+                test_case["subcommunity_species_to_normalized_subcommunity_abundance"][
+                    (subcommunity, species)
+                ],
+            ), (
+                f"\n(i,j): {(i,j)};"
+                f"\n(species,subcommunity): {(species,subcommunity)};"
+            )
+
+
+def arrange_values(ordered_names, name_to_row):
+    """Arranges matrix rows according to ordered_names ordering."""
+    matrix = empty(
+        shape=(len(ordered_names), name_to_row[ordered_names[0]].shape[0]),
+        dtype=float64,
+    )
+    for i, name in enumerate(ordered_names):
+        matrix[i][:] = name_to_row[name]
+    return matrix
+
+
+SIMILARITY_FROM_FILE_TEST_CASES = [
+    {
+        "description": "tsv file; 2 communities",
+        "similarity_matrix_filepath": "similarities_file.tsv",
+        "expected_species_order": array(["species_3", "species_1", "species_2"]),
+        "species_to_relative_abundances": {
+            "species_1": array([1 / 1000, 1 / 100]),
+            "species_2": array([1 / 10, 1 / 1]),
+            "species_3": array([10, 100]),
+        },
+        "species_to_weighted_similarities": {
+            "species_1": array([1.051, 10.51]),
+            "species_2": array([2.1005, 21.005]),
+            "species_3": array([10.0201, 100.201]),
+        },
+        "similarities_filecontents": (
+            "species_3\tspecies_1\tspecies_2\n"
+            "1\t0.1\t0.2\n"
+            "0.1\t1\t0.5\n"
+            "0.2\t0.5\t1\n"
+        ),
+    },
+    {
+        "description": "csv file; 2 communities",
+        "similarity_matrix_filepath": "similarities_file.csv",
+        "expected_species_order": array(["species_3", "species_1", "species_2"]),
+        "species_to_relative_abundances": {
+            "species_1": array([1 / 1000, 1 / 100]),
+            "species_2": array([1 / 10, 1 / 1]),
+            "species_3": array([10, 100]),
+        },
+        "species_to_weighted_similarities": {
+            "species_1": array([1.051, 10.51]),
+            "species_2": array([2.1005, 21.005]),
+            "species_3": array([10.0201, 100.201]),
+        },
+        "similarities_filecontents": (
+            "species_3,species_1,species_2\n" "1,0.1,0.2\n" "0.1,1,0.5\n" "0.2,0.5,1\n"
+        ),
+    },
+    {
+        "description": "no file extension; 1 community",
+        "similarity_matrix_filepath": "similarities_file",
+        "expected_species_order": array(["species_3", "species_1", "species_2"]),
+        "species_to_relative_abundances": {
+            "species_1": array([1 / 1000]),
+            "species_2": array([1 / 10]),
+            "species_3": array([10]),
+        },
+        "species_to_weighted_similarities": {
+            "species_1": array([1.051]),
+            "species_2": array([2.1005]),
+            "species_3": array([10.0201]),
+        },
+        "similarities_filecontents": (
+            "species_3\tspecies_1\tspecies_2\n"
+            "1\t0.1\t0.2\n"
+            "0.1\t1\t0.5\n"
+            "0.2\t0.5\t1\n"
+        ),
+    },
+]
+
+
+class TestSimilarityFromFile:
+    """Tests diversity.metacommunity.Similarity."""
+
+    @mark.parametrize("test_case", SIMILARITY_FROM_FILE_TEST_CASES)
+    def test_init(self, test_case, tmp_path):
+        """Tests initializer."""
+        test_case["similarity_matrix_filepath"] = (
+            tmp_path / test_case["similarity_matrix_filepath"]
+        )
+        with open(test_case["similarity_matrix_filepath"], "w") as file:
+            file.write(test_case["similarities_filecontents"])
+        if test_case["similarity_matrix_filepath"].suffix == "":
+            filterwarnings("ignore", category=ArgumentWarning)
+        similarity = SimilarityFromFile(
+            similarity_matrix_filepath=test_case["similarity_matrix_filepath"],
+        )
+        resetwarnings()
+        assert (similarity.species_order == test_case["expected_species_order"]).all()
+
+    @mark.parametrize("test_case", SIMILARITY_FROM_FILE_TEST_CASES)
+    def test_calculate_weighted_similarities(self, test_case, tmp_path):
+        """Tests .calculate_weighted_similarities."""
+        test_case["similarity_matrix_filepath"] = (
+            tmp_path / test_case["similarity_matrix_filepath"]
+        )
+        with open(test_case["similarity_matrix_filepath"], "w") as file:
+            file.write(test_case["similarities_filecontents"])
+        if test_case["similarity_matrix_filepath"].suffix == "":
+            filterwarnings("ignore", category=ArgumentWarning)
+        similarity = SimilarityFromFile(
+            similarity_matrix_filepath=test_case["similarity_matrix_filepath"],
+        )
+        resetwarnings()
+        relative_abundances = arrange_values(
+            test_case["expected_species_order"],
+            test_case["species_to_relative_abundances"],
+        )
+        weighted_similarities = similarity.calculate_weighted_similarities(
+            relative_abundances
+        )
+        expected_weighted_similarities = arrange_values(
+            test_case["expected_species_order"],
+            test_case["species_to_weighted_similarities"],
+        )
+        assert weighted_similarities.shape == relative_abundances.shape
+        assert allclose(weighted_similarities, expected_weighted_similarities)
+        with open(test_case["similarity_matrix_filepath"], "r") as file:
+            similarities_filecontents = file.read()
+        assert similarities_filecontents == test_case["similarities_filecontents"]
+
+
+SIMILARITY_FROM_FUNCTION_TEST_CASES = [
+    {
+        "description": "similarity function; 2 communities; tsv similarities file",
+        "similarity_function": lambda a, b: 1 / sum(a * b),
+        "features": array([[1, 2], [3, 5], [7, 11]]),
+        "species_order": array(["species_3", "species_1", "species_2"]),
+        "expected_species_order": array(["species_3", "species_1", "species_2"]),
+        "species_to_relative_abundances": {
+            "species_3": array([1 / 1000, 1 / 100]),
+            "species_1": array([1 / 10, 1 / 1]),
+            "species_2": array([10, 100]),
+        },
+        "species_to_weighted_similarities": {
+            "species_3": array([0.35271989, 3.52719894]),
+            "species_1": array([0.13459705, 1.34597047]),
+            "species_2": array([0.0601738, 0.60173802]),
+        },
+        # similarity matrix
+        # [0.2          , 0.07692307692, 0.03448275862]
+        # [0.07692307692, 0.02941176471, 0.01315789474]
+        # [0.03448275862, 0.01315789474, 0.005882352941]
+    },
+    {
+        "description": "similarity function; 2 communities; csv similarities file",
+        "similarity_function": lambda a, b: 1 / sum(a * b),
+        "features": array([[1, 2], [3, 5], [7, 11]]),
+        "species_order": array(["species_3", "species_1", "species_2"]),
+        "expected_species_order": array(["species_3", "species_1", "species_2"]),
+        "species_to_relative_abundances": {
+            "species_3": array([1 / 1000, 1 / 100]),
+            "species_1": array([1 / 10, 1 / 1]),
+            "species_2": array([10, 100]),
+        },
+        "species_to_weighted_similarities": {
+            "species_3": array([0.35271989, 3.52719894]),
+            "species_1": array([0.13459705, 1.34597047]),
+            "species_2": array([0.0601738, 0.60173802]),
+        },
+        # similarity matrix
+        # [0.2          , 0.07692307692, 0.03448275862]
+        # [0.07692307692, 0.02941176471, 0.01315789474]
+        # [0.03448275862, 0.01315789474, 0.005882352941]
+    },
+    {
+        "description": "similarity function; 1 community; similarities file without extension",
+        "similarity_function": lambda a, b: 1 / sum(a * b),
+        "features": array([[1, 2], [3, 5], [7, 11]]),
+        "species_order": array(["species_3", "species_1", "species_2"]),
+        "expected_species_order": array(["species_3", "species_1", "species_2"]),
+        "species_to_relative_abundances": {
+            "species_3": array([1 / 1000]),
+            "species_1": array([1 / 10]),
+            "species_2": array([10]),
+        },
+        "species_to_weighted_similarities": {
+            "species_3": array([0.35271989]),
+            "species_1": array([0.13459705]),
+            "species_2": array([0.0601738]),
+        },
+        # similarity matrix
+        # [0.2          , 0.07692307692, 0.03448275862]
+        # [0.07692307692, 0.02941176471, 0.01315789474]
+        # [0.03448275862, 0.01315789474, 0.005882352941]
+    },
+]
+
+
+class TestSimilarityFromFunction:
+    """Tests diversity.metacommunity.Similarity."""
+
+    @mark.parametrize("test_case", SIMILARITY_FROM_FUNCTION_TEST_CASES)
+    def test_init(self, test_case, tmp_path):
+        """Tests initializer."""
+        similarity = SimilarityFromFunction(
+            similarity_function=test_case["similarity_function"],
+            features=test_case["features"],
+            species_order=test_case["species_order"],
+        )
+        assert (similarity.species_order == test_case["expected_species_order"]).all()
+
+    @mark.parametrize("test_case", SIMILARITY_FROM_FUNCTION_TEST_CASES)
+    def test_calculate_weighted_similarities(self, test_case, tmp_path):
+        """Tests .calculate_weighted_similarities."""
+        similarity = SimilarityFromFunction(
+            similarity_function=test_case["similarity_function"],
+            features=test_case["features"],
+            species_order=test_case["species_order"],
+        )
+        relative_abundances = arrange_values(
+            test_case["expected_species_order"],
+            test_case["species_to_relative_abundances"],
+        )
+        weighted_similarities = similarity.calculate_weighted_similarities(
+            relative_abundances
+        )
+        expected_weighted_similarities = arrange_values(
+            test_case["expected_species_order"],
+            test_case["species_to_weighted_similarities"],
+        )
+        assert weighted_similarities.shape == relative_abundances.shape
+        assert allclose(weighted_similarities, expected_weighted_similarities)
+
+
+SIMILARITY_FROM_MEMORY_TEST_CASES = [
+    {
+        "description": "similarities in memory; 2 communities",
+        "similarity_matrix": array(
+            [
+                [1, 0.5, 0.1],
+                [0.5, 1, 0.2],
+                [0.1, 0.2, 1],
+            ]
+        ),
+        "species_order": array(["species_1", "species_2", "species_3"]),
+        "expected_species_order": array(["species_1", "species_2", "species_3"]),
+        "species_to_relative_abundances": {
+            "species_1": array([1 / 1000, 1 / 100]),
+            "species_2": array([1 / 10, 1 / 1]),
+            "species_3": array([10, 100]),
+        },
+        "species_to_weighted_similarities": {
+            "species_1": array([1.051, 10.51]),
+            "species_2": array([2.1005, 21.005]),
+            "species_3": array([10.0201, 100.201]),
+        },
+    },
+    {
+        "description": "similarities in memory; 1 community",
+        "similarity_matrix": array(
+            [
+                [1, 0.5, 0.1],
+                [0.5, 1, 0.2],
+                [0.1, 0.2, 1],
+            ]
+        ),
+        "similarities_filepath": None,
+        "similarity_function": None,
+        "features": None,
+        "species_order": array(["species_1", "species_2", "species_3"]),
+        "expected_species_order": array(["species_1", "species_2", "species_3"]),
+        "species_to_relative_abundances": {
+            "species_1": array([1 / 1000]),
+            "species_2": array([1 / 10]),
+            "species_3": array([10]),
+        },
+        "species_to_weighted_similarities": {
+            "species_1": array([1.051]),
+            "species_2": array([2.1005]),
+            "species_3": array([10.0201]),
+        },
+        "similarities_filecontents": None,
+    },
+]
+
+
+class TestSimilarityFromMemory:
+    """Tests diversity.metacommunity.Similarity."""
+
+    def arrange_values(self, ordered_names, name_to_row):
+        """Arranges rows according to ordered_names ordering."""
+        matrix = empty(
+            shape=(len(ordered_names), len(name_to_row[ordered_names[0]])),
+            dtype=float64,
+        )
+        for i, name in enumerate(ordered_names):
+            matrix[i] = name_to_row[name]
+        return matrix
+
+    @mark.parametrize("test_case", SIMILARITY_FROM_MEMORY_TEST_CASES)
+    def test_init(self, test_case):
+        """Tests initializer."""
+        similarity = SimilarityFromMemory(
+            similarity_matrix=test_case["similarity_matrix"],
+            species_order=test_case["species_order"],
+        )
+        assert (similarity.species_order == test_case["expected_species_order"]).all()
+
+    @mark.parametrize("test_case", SIMILARITY_FROM_MEMORY_TEST_CASES)
+    def test_calculate_weighted_similarities(self, test_case, tmp_path):
+        """Tests .calculate_weighted_similarities."""
+        similarity = SimilarityFromMemory(
+            similarity_matrix=test_case["similarity_matrix"],
+            species_order=test_case["species_order"],
+        )
+        relative_abundances = self.arrange_values(
+            test_case["expected_species_order"],
+            test_case["species_to_relative_abundances"],
+        )
+        weighted_similarities = similarity.calculate_weighted_similarities(
+            relative_abundances
+        )
+        expected_weighted_similarities = self.arrange_values(
+            test_case["expected_species_order"],
+            test_case["species_to_weighted_similarities"],
+        )
+        assert weighted_similarities.shape == relative_abundances.shape
+        assert allclose(weighted_similarities, expected_weighted_similarities)
+
+
+CREATE_SIMILARITY_TEST_CASES = [
+    {
+        "description": "from memory; no additional parameters",
+        "similarity_matrix": array([[1.0, 0.5, 0.1], [0.5, 1.0, 0.2], [0.1, 0.2, 1.0]]),
+        "species_order": array(["species_1", "species_2", "species_3"]),
+        "similarity_matrix_filepath": None,
+        "similarity_function": None,
+        "features": None,
+        "similarity_matrix_file_content": None,
+        "expected_type": SimilarityFromMemory,
+    },
+    {
+        "description": "from memory; with additional ignored parameters",
+        "similarity_matrix": array([[1.0, 0.5, 0.1], [0.5, 1.0, 0.2], [0.1, 0.2, 1.0]]),
+        "species_order": array(["species_1", "species_2", "species_3"]),
+        "similarity_matrix_filepath": None,
+        "similarity_function": None,
+        "features": array(
+            [[3.2, 5.3, 1.2, 9.4], [4.2, 7.4, 9.5, 7.3], [4.2, 6.2, 6.4, 7.3]]
+        ),
+        "similarity_matrix_file_content": None,
+        "expected_type": SimilarityFromMemory,
+    },
+    {
+        "description": "from function",
+        "similarity_matrix": None,
+        "species_order": array(["species_1", "species_2", "species_3"]),
+        "similarity_matrix_filepath": "foo_similarities.tsv",
+        "similarity_function": lambda x, y: 1 / sum(x * y),
+        "features": array(
+            [[3.2, 5.3, 1.2, 9.4], [4.2, 7.4, 9.5, 7.3], [4.2, 6.2, 6.4, 7.3]]
+        ),
+        "similarity_matrix_file_content": None,
+        "expected_type": SimilarityFromFunction,
+    },
+    {
+        "description": "from file; no additional parameters",
+        "similarity_matrix": None,
+        "species_order": None,
+        "similarity_matrix_filepath": "foo_similarities.tsv",
+        "similarity_function": None,
+        "features": None,
+        "similarity_matrix_file_content": "s3\ts1\ts2\n",
+        "expected_type": SimilarityFromFile,
+    },
+    {
+        "description": "from file; with additional ignored parameters",
+        "similarity_matrix": None,
+        "species_order": array(["species_1", "species_2", "species_3"]),
+        "similarity_matrix_filepath": "foo_similarities.tsv",
+        "similarity_function": None,
+        "features": None,
+        "similarity_matrix_file_content": "s3\ts1\ts2\n",
+        "expected_type": SimilarityFromFile,
+    },
+]
+
+
+class TestCreateSimilarity:
+    """Tests metacommunity.create_similarity."""
+
+    @mark.parametrize("test_case", CREATE_SIMILARITY_TEST_CASES)
+    def test_create_similarity(self, test_case, tmp_path):
+        """Tests create_similarity test cases."""
+        if test_case["expected_type"] == SimilarityFromFile:
+            absolute_path = tmp_path / test_case["similarity_matrix_filepath"]
+            with open(absolute_path, "w") as file:
+                file.write(test_case["similarity_matrix_file_content"])
+            test_case["similarity_matrix_filepath"] = absolute_path
+        similarity_object = make_similarity(
+            similarity_matrix=test_case["similarity_matrix"],
+            species_order=test_case["species_order"],
+            similarity_matrix_filepath=test_case["similarity_matrix_filepath"],
+            similarity_function=test_case["similarity_function"],
+            features=test_case["features"],
+        )
+        assert isinstance(similarity_object, test_case["expected_type"])
+
+
+METACOMMUNITY_TEST_CASES = [
+    {
+        "description": "disjoint communities; uniform counts; uniform inter-community similarities; viewpoint 0.",
+        "similarity": SimilarityFromMemory(
+            similarity_matrix=array(
+                [
+                    [1.0, 0.5, 0.5, 0.7, 0.7, 0.7],
+                    [0.5, 1.0, 0.5, 0.7, 0.7, 0.7],
+                    [0.5, 0.5, 1.0, 0.7, 0.7, 0.7],
+                    [0.7, 0.7, 0.7, 1.0, 0.5, 0.5],
+                    [0.7, 0.7, 0.7, 0.5, 1.0, 0.5],
+                    [0.7, 0.7, 0.7, 0.5, 0.5, 1.0],
+                ]
+            ),
+            species_order=array(
+                [
+                    "species_1",
+                    "species_2",
+                    "species_3",
+                    "species_4",
+                    "species_5",
+                    "species_6",
+                ]
+            ),
+        ),
+        "abundance": Abundance(
+            counts=array(
+                [
+                    ["subcommunity_1", "species_1", "1"],
+                    ["subcommunity_1", "species_2", "1"],
+                    ["subcommunity_1", "species_3", "1"],
+                    ["subcommunity_2", "species_4", "1"],
+                    ["subcommunity_2", "species_5", "1"],
+                    ["subcommunity_2", "species_6", "1"],
+                ]
+            ),
+            species_order=array(
+                [
+                    "species_1",
+                    "species_2",
+                    "species_3",
+                    "species_4",
+                    "species_5",
+                    "species_6",
+                ]
+            ),
+        ),
+        "viewpoint": 0,
+        "subcommunity_to_alpha": {
+            "subcommunity_1": 3.0,
+            "subcommunity_2": 3.0,
+        },
+        "subcommunity_to_rho": {
+            "subcommunity_1": 2.05,
+            "subcommunity_2": 2.05,
+        },
+        "subcommunity_to_beta": {
+            "subcommunity_1": 0.487805,
+            "subcommunity_2": 0.487805,
+        },
+        "subcommunity_to_gamma": {
+            "subcommunity_1": 1.463415,
+            "subcommunity_2": 1.463415,
+        },
+        "subcommunity_to_normalized_alpha": {
+            "subcommunity_1": 1.5,
+            "subcommunity_2": 1.5,
+        },
+        "subcommunity_to_normalized_rho": {
+            "subcommunity_1": 1.025,
+            "subcommunity_2": 1.025,
+        },
+        "subcommunity_to_normalized_beta": {
+            "subcommunity_1": 0.97561,
+            "subcommunity_2": 0.97561,
+        },
+        "metacommunity_alpha": 3.0,
+        "metacommunity_rho": 2.05,
+        "metacommunity_beta": 0.487805,
+        "metacommunity_gamma": 1.463415,
+        "metacommunity_normalized_alpha": 1.5,
+        "metacommunity_normalized_rho": 1.025,
+        "metacommunity_normalized_beta": 0.97561,
+    },
+    {
+        "description": "overlapping communities; non-uniform counts; non-uniform inter-community similarities; viewpoint 2.",
+        "similarity": SimilarityFromMemory(
+            similarity_matrix=array(
+                [[1.0, 0.5, 0.1], [0.5, 1.0, 0.2], [0.1, 0.2, 1.0]]
+            ),
+            species_order=array(
+                [
+                    "species_1",
+                    "species_2",
+                    "species_3",
+                ]
+            ),
+        ),
+        "abundance": Abundance(
+            counts=array(
+                [
+                    ["subcommunity_1", "species_1", "2"],
+                    ["subcommunity_1", "species_2", "3"],
+                    ["subcommunity_2", "species_1", "5"],
+                    ["subcommunity_2", "species_3", "1"],
+                ]
+            ),
+            species_order=array(
+                [
+                    "species_1",
+                    "species_2",
+                    "species_3",
+                ]
+            ),
+        ),
+        "viewpoint": 2,
+        "subcommunity_to_alpha": {
+            "subcommunity_1": 2.89473684,
+            "subcommunity_2": 2.44444444,
+        },
+        "subcommunity_to_rho": {
+            "subcommunity_1": 1.91938708,
+            "subcommunity_2": 1.65870021,
+        },
+        "subcommunity_to_beta": {
+            "subcommunity_1": 0.52099965,
+            "subcommunity_2": 0.6028817,
+        },
+        "subcommunity_to_gamma": {
+            "subcommunity_1": 1.47453083,
+            "subcommunity_2": 1.45695364,
+        },
+        "subcommunity_to_normalized_alpha": {
+            "subcommunity_1": 1.31578947,
+            "subcommunity_2": 1.33333333,
+        },
+        "subcommunity_to_normalized_rho": {
+            "subcommunity_1": 0.87244867,
+            "subcommunity_2": 0.90474557,
+        },
+        "subcommunity_to_normalized_beta": {
+            "subcommunity_1": 1.14619924,
+            "subcommunity_2": 1.10528311,
+        },
+        "metacommunity_alpha": 2.630434782608696,
+        "metacommunity_rho": 1.7678383245514573,
+        "metacommunity_beta": 0.5626846957892072,
+        "metacommunity_gamma": 1.4648910411622276,
+        "metacommunity_normalized_alpha": 1.3253012048192772,
+        "metacommunity_normalized_rho": 0.8897736389973379,
+        "metacommunity_normalized_beta": 1.1235132485772303,
+    },
+]
+
+
+class TestMetacommunity:
+    """Tests metacommunity.Metacommunity."""
+
+    def make_metacommunity(self, test_case):
+        """Simple initializer that doesn't modify test_case."""
+        test_case_ = deepcopy(test_case)
+        return Metacommunity(
+            similarity=test_case_["similarity"], abundance=test_case_["abundance"]
+        )
+
+    @mark.parametrize("test_case", METACOMMUNITY_TEST_CASES)
+    def test_subcommunity_alpha(self, test_case):
+        metacommunity = self.make_metacommunity(test_case)
+        subcommunity_order = test_case["abundance"].subcommunity_order
+        subcommunity_alpha = metacommunity.subcommunity_alpha(test_case["viewpoint"])
+        for i, subcommunity in enumerate(subcommunity_order):
+            assert isclose(
+                subcommunity_alpha[i], test_case["subcommunity_to_alpha"][subcommunity]
+            ), f"\n(i, subcommunity): {(i,subcommunity)}"
+
+    @mark.parametrize("test_case", METACOMMUNITY_TEST_CASES)
+    def test_subcommunity_rho(self, test_case):
+        metacommunity = self.make_metacommunity(test_case)
+        subcommunity_order = test_case["abundance"].subcommunity_order
+        subcommunity_rho = metacommunity.subcommunity_rho(test_case["viewpoint"])
+        for i, subcommunity in enumerate(subcommunity_order):
+            assert isclose(
+                subcommunity_rho[i], test_case["subcommunity_to_rho"][subcommunity]
+            ), f"\n(i, subcommunity): {(i,subcommunity)}"
+
+    @mark.parametrize("test_case", METACOMMUNITY_TEST_CASES)
+    def test_subcommunity_beta(self, test_case):
+        metacommunity = self.make_metacommunity(test_case)
+        subcommunity_order = test_case["abundance"].subcommunity_order
+        subcommunity_beta = metacommunity.subcommunity_beta(test_case["viewpoint"])
+        for i, subcommunity in enumerate(subcommunity_order):
+            assert isclose(
+                subcommunity_beta[i], test_case["subcommunity_to_beta"][subcommunity]
+            ), f"\n(i, subcommunity): {(i,subcommunity)}"
+
+    @mark.parametrize("test_case", METACOMMUNITY_TEST_CASES)
+    def test_subcommunity_gamma(self, test_case):
+        metacommunity = self.make_metacommunity(test_case)
+        subcommunity_order = test_case["abundance"].subcommunity_order
+        subcommunity_gamma = metacommunity.subcommunity_gamma(test_case["viewpoint"])
+        for i, subcommunity in enumerate(subcommunity_order):
+            assert isclose(
+                subcommunity_gamma[i], test_case["subcommunity_to_gamma"][subcommunity]
+            ), f"\n(i, subcommunity): {(i,subcommunity)}"
+
+    @mark.parametrize("test_case", METACOMMUNITY_TEST_CASES)
+    def test_normalized_subcommunity_alpha(self, test_case):
+        metacommunity = self.make_metacommunity(test_case)
+        subcommunity_order = test_case["abundance"].subcommunity_order
+        normalized_subcommunity_alpha = metacommunity.normalized_subcommunity_alpha(
+            test_case["viewpoint"]
+        )
+        for i, subcommunity in enumerate(subcommunity_order):
+            assert isclose(
+                normalized_subcommunity_alpha[i],
+                test_case["subcommunity_to_normalized_alpha"][subcommunity],
+            ), f"\n(i, subcommunity): {(i,subcommunity)}"
+
+    @mark.parametrize("test_case", METACOMMUNITY_TEST_CASES)
+    def test_normalized_subcommunity_rho(self, test_case):
+        metacommunity = self.make_metacommunity(test_case)
+        subcommunity_order = test_case["abundance"].subcommunity_order
+        normalized_subcommunity_rho = metacommunity.normalized_subcommunity_rho(
+            test_case["viewpoint"]
+        )
+        for i, subcommunity in enumerate(subcommunity_order):
+            assert isclose(
+                normalized_subcommunity_rho[i],
+                test_case["subcommunity_to_normalized_rho"][subcommunity],
+            ), f"\n(i, subcommunity): {(i,subcommunity)}"
+
+    @mark.parametrize("test_case", METACOMMUNITY_TEST_CASES)
+    def test_normalized_subcommunity_beta(self, test_case):
+        metacommunity = self.make_metacommunity(test_case)
+        subcommunity_order = test_case["abundance"].subcommunity_order
+        normalized_subcommunity_beta = metacommunity.normalized_subcommunity_beta(
+            test_case["viewpoint"]
+        )
+        for i, subcommunity in enumerate(subcommunity_order):
+            assert isclose(
+                normalized_subcommunity_beta[i],
+                test_case["subcommunity_to_normalized_beta"][subcommunity],
+            ), f"\n(i, subcommunity): {(i,subcommunity)}"
+
+    @mark.parametrize("test_case", METACOMMUNITY_TEST_CASES)
+    def test_metacommunity_alpha(self, test_case):
+        metacommunity = self.make_metacommunity(test_case)
+        metacommunity_alpha = metacommunity.metacommunity_alpha(test_case["viewpoint"])
+        assert isclose(metacommunity_alpha, test_case["metacommunity_alpha"])
+
+    @mark.parametrize("test_case", METACOMMUNITY_TEST_CASES)
+    def test_metacommunity_rho(self, test_case):
+        metacommunity = self.make_metacommunity(test_case)
+        metacommunity_rho = metacommunity.metacommunity_rho(test_case["viewpoint"])
+        assert isclose(metacommunity_rho, test_case["metacommunity_rho"])
+
+    @mark.parametrize("test_case", METACOMMUNITY_TEST_CASES)
+    def test_metacommunity_beta(self, test_case):
+        metacommunity = self.make_metacommunity(test_case)
+        metacommunity_beta = metacommunity.metacommunity_beta(test_case["viewpoint"])
+        assert isclose(metacommunity_beta, test_case["metacommunity_beta"])
+
+    @mark.parametrize("test_case", METACOMMUNITY_TEST_CASES)
+    def test_metacommunity_gamma(self, test_case):
+        metacommunity = self.make_metacommunity(test_case)
+        metacommunity_gamma = metacommunity.metacommunity_gamma(test_case["viewpoint"])
+        assert isclose(metacommunity_gamma, test_case["metacommunity_gamma"])
+
+    @mark.parametrize("test_case", METACOMMUNITY_TEST_CASES)
+    def test_metacommunity_normalized_alpha(self, test_case):
+        metacommunity = self.make_metacommunity(test_case)
+        metacommunity_normalized_alpha = metacommunity.normalized_metacommunity_alpha(
+            test_case["viewpoint"]
+        )
+        assert isclose(
+            metacommunity_normalized_alpha, test_case["metacommunity_normalized_alpha"]
+        )
+
+    @mark.parametrize("test_case", METACOMMUNITY_TEST_CASES)
+    def test_metacommunity_normalized_rho(self, test_case):
+        metacommunity = self.make_metacommunity(test_case)
+        metacommunity_normalized_rho = metacommunity.normalized_metacommunity_rho(
+            test_case["viewpoint"]
+        )
+        assert isclose(
+            metacommunity_normalized_rho, test_case["metacommunity_normalized_rho"]
+        )
+
+    @mark.parametrize("test_case", METACOMMUNITY_TEST_CASES)
+    def test_metacommunity_normalized_beta(self, test_case):
+        metacommunity = self.make_metacommunity(test_case)
+        metacommunity_normalized_beta = metacommunity.normalized_metacommunity_beta(
+            test_case["viewpoint"]
+        )
+        assert isclose(
+            metacommunity_normalized_beta, test_case["metacommunity_normalized_beta"]
+        )
+
+
+MAKE_METACOMMUNITY_TEST_CASES = [
+    {
+        "description": "numpy.ndarray counts; similarities from memory; no additional parameters",
+        "counts": array(
+            [
+                ["community_1", "species_1", "2"],
+                ["community_1", "species_2", "3"],
+                ["community_2", "species_1", "4"],
+                ["community_2", "species_3", "1"],
+            ]
+        ),
+        "similarity_matrix": array([[1.0, 0.5, 0.1], [0.5, 1.0, 0.2], [0.1, 0.2, 1.0]]),
+        "species_order": array(["species_1", "species_2", "species_3"]),
+        "similarity_matrix_filepath": None,
+        "similarity_function": None,
+        "features": None,
+        "expected_similarity_type": SimilarityFromMemory,
+    },
+    {
+        "description": "numpy.ndarray counts; similarities from memory; with additional ignored parameters",
+        "counts": array(
+            [
+                ["community_1", "species_1", "2"],
+                ["community_1", "species_2", "3"],
+                ["community_2", "species_1", "4"],
+                ["community_2", "species_3", "1"],
+            ]
+        ),
+        "similarity_matrix": array([[1.0, 0.5, 0.1], [0.5, 1.0, 0.2], [0.1, 0.2, 1.0]]),
+        "species_order": array(["species_1", "species_2", "species_3"]),
+        "similarity_matrix_filepath": None,
+        "similarity_function": None,
+        "features": array(
+            [[3.2, 5.3, 1.2, 9.4], [4.2, 7.4, 9.5, 7.3], [4.2, 6.2, 6.4, 7.3]]
+        ),
+        "expected_similarity_type": SimilarityFromMemory,
+    },
+    {
+        "description": "numpy.ndarray counts; similarities from function",
+        "counts": array(
+            [
+                ["community_1", "species_1", "2"],
+                ["community_1", "species_2", "3"],
+                ["community_2", "species_1", "4"],
+                ["community_2", "species_3", "1"],
+            ]
+        ),
+        "similarity_matrix": None,
+        "species_order": array(["species_1", "species_2", "species_3"]),
+        "similarity_matrix_filepath": "foo_similarities.tsv",
+        "similarity_function": lambda x, y: 1 / sum(x * y),
+        "features": array(
+            [[3.2, 5.3, 1.2, 9.4], [4.2, 7.4, 9.5, 7.3], [4.2, 6.2, 6.4, 7.3]]
+        ),
+        "expected_similarity_type": SimilarityFromFunction,
+    },
+    {
+        "description": "numpy.ndarray counts; similarities from file; no additional parameters",
+        "counts": array(
+            [
+                ["community_1", "species_1", "2"],
+                ["community_1", "species_2", "3"],
+                ["community_2", "species_1", "4"],
+                ["community_2", "species_3", "1"],
+            ]
+        ),
+        "similarity_matrix": None,
+        "species_order": None,
+        "similarity_matrix_filepath": "foo_similarities.tsv",
+        "similarity_function": None,
+        "features": None,
+        "expected_similarity_type": SimilarityFromFile,
+    },
+    {
+        "description": "numpy.ndarray counts; similarities from memory; with additional ignored parameters",
+        "counts": array(
+            [
+                ["community_1", "species_1", "2"],
+                ["community_1", "species_2", "3"],
+                ["community_2", "species_1", "4"],
+                ["community_2", "species_3", "1"],
+            ]
+        ),
+        "similarity_matrix": None,
+        "species_order": array(["species_1", "species_2", "species_3"]),
+        "similarity_matrix_filepath": "foo_similarities.tsv",
+        "similarity_function": None,
+        "features": None,
+        "expected_similarity_type": SimilarityFromFile,
+    },
+    {
+        "description": "pandas.DataFrame counts; similarities from memory; no additional parameters",
+        "counts": DataFrame(
+            {
+                "community": [
+                    "community_1",
+                    "community_1",
+                    "community_2",
+                    "community_2",
+                ],
+                "species": ["species_1", "species_2", "species_1", "species_3"],
+                "count": [2, 3, 4, 1],
+            }
+        ),
+        "similarity_matrix": array([[1.0, 0.5, 0.1], [0.5, 1.0, 0.2], [0.1, 0.2, 1.0]]),
+        "species_order": array(["species_1", "species_2", "species_3"]),
+        "similarity_matrix_filepath": None,
+        "similarity_function": None,
+        "features": None,
+        "expected_similarity_type": SimilarityFromMemory,
+    },
+    {
+        "description": "pandas.DataFrame counts; similarities from memory; with additional ignored parameters",
+        "counts": DataFrame(
+            {
+                "community": [
+                    "community_1",
+                    "community_1",
+                    "community_2",
+                    "community_2",
+                ],
+                "species": ["species_1", "species_2", "species_1", "species_3"],
+                "count": [2, 3, 4, 1],
+            }
+        ),
+        "similarity_matrix": array([[1.0, 0.5, 0.1], [0.5, 1.0, 0.2], [0.1, 0.2, 1.0]]),
+        "species_order": array(["species_1", "species_2", "species_3"]),
+        "similarity_matrix_filepath": None,
+        "similarity_function": None,
+        "features": array(
+            [[3.2, 5.3, 1.2, 9.4], [4.2, 7.4, 9.5, 7.3], [4.2, 6.2, 6.4, 7.3]]
+        ),
+        "expected_similarity_type": SimilarityFromMemory,
+    },
+    {
+        "description": "pandas.DataFrame counts; similarities from function",
+        "counts": DataFrame(
+            {
+                "community": [
+                    "community_1",
+                    "community_1",
+                    "community_2",
+                    "community_2",
+                ],
+                "species": ["species_1", "species_2", "species_1", "species_3"],
+                "count": [2, 3, 4, 1],
+            }
+        ),
+        "similarity_matrix": None,
+        "species_order": array(["species_1", "species_2", "species_3"]),
+        "similarity_matrix_filepath": "foo_similarities.tsv",
+        "similarity_function": lambda x, y: 1 / sum(x * y),
+        "features": array(
+            [[3.2, 5.3, 1.2, 9.4], [4.2, 7.4, 9.5, 7.3], [4.2, 6.2, 6.4, 7.3]]
+        ),
+        "expected_similarity_type": SimilarityFromFunction,
+    },
+    {
+        "description": "pandas.DataFrame counts; similarities from file; no additional parameters",
+        "counts": DataFrame(
+            {
+                "community": [
+                    "community_1",
+                    "community_1",
+                    "community_2",
+                    "community_2",
+                ],
+                "species": ["species_1", "species_2", "species_1", "species_3"],
+                "count": [2, 3, 4, 1],
+            }
+        ),
+        "similarity_matrix": None,
+        "species_order": None,
+        "similarity_matrix_filepath": "foo_similarities.tsv",
+        "similarity_function": None,
+        "features": None,
+        "expected_similarity_type": SimilarityFromFile,
+    },
+    {
+        "description": "pandas.DataFrame counts; similarities from memory; with additional ignored parameters",
+        "counts": DataFrame(
+            {
+                "community": [
+                    "community_1",
+                    "community_1",
+                    "community_2",
+                    "community_2",
+                ],
+                "species": ["species_1", "species_2", "species_1", "species_3"],
+                "count": [2, 3, 4, 1],
+            }
+        ),
+        "similarity_matrix": None,
+        "species_order": array(["species_1", "species_2", "species_3"]),
+        "similarity_matrix_filepath": "foo_similarities.tsv",
+        "similarity_function": None,
+        "features": None,
+        "expected_similarity_type": SimilarityFromFile,
+    },
+]
+
+
+class TestMakeMetacommunity:
+    """Tests metacommunity.make_metacommunity."""
+
+    @mark.parametrize("test_case", MAKE_METACOMMUNITY_TEST_CASES)
+    def test_create_similarity(self, test_case, tmp_path):
+        """Tests make_metacommunity test cases."""
+        if test_case["similarity_matrix_filepath"] is not None:
+            test_case["similarity_matrix_filepath"] = (
+                tmp_path / test_case["similarity_matrix_filepath"]
+            )
+        if test_case["expected_similarity_type"] == SimilarityFromFile:
+            with open(test_case["similarity_matrix_filepath"], "w") as file:
+                if isinstance(test_case["counts"], ndarray):
+                    unique_species = unique(test_case["counts"][:, 1])
+                else:
+                    unique_species = unique(test_case["counts"].iloc[:, 1].to_numpy())
+
+                file.write("\t".join(species for species in unique_species) + "\n")
+        metacommunity = make_metacommunity(
+            counts=test_case["counts"],
+            similarity_matrix=test_case["similarity_matrix"],
+            similarity_matrix_filepath=test_case["similarity_matrix_filepath"],
+            similarity_function=test_case["similarity_function"],
+            features=test_case["features"],
+            species_order=test_case["species_order"],
+        )
+        assert isinstance(
+            metacommunity.similarity, test_case["expected_similarity_type"]
+        )
diff --git a/tests/utilities_test.py b/tests/utilities_test.py
new file mode 100644
index 00000000..dfd48f86
--- /dev/null
+++ b/tests/utilities_test.py
@@ -0,0 +1,834 @@
+"""Tests for diversity.utilities."""
+from copy import deepcopy
+
+from numpy import allclose, array, inf
+from pytest import mark, raises
+
+from diversity.utilities import InvalidArgumentError, power_mean, unique_correspondence
+
+
+POWER_MEAN_TEST_CASES = [
+    {
+        "description": "No zero weights; nonzero order; -100 <= order <= 100; 2-d data.",
+        "order": 3,
+        "weights": array(
+            [
+                [1.0e-07, 1.0e-07, 1.0e-07],
+                [2.0e-08, 2.0e-08, 2.0e-08],
+                [1.0e-01, 1.0e-01, 1.0e-01],
+            ]
+        ),
+        "items": array(
+            [
+                [0.30, 0.69, 0.53],
+                [0.74, 0.14, 0.53],
+                [0.71, 0.13, 0.86],
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": array([0.32955284, 0.06034367, 0.39917668]),
+        "expect_raise": False,
+    },
+    {
+        "description": "No zero weights; nonzero order; -100 <= order <= 100; 2-d data.",
+        "order": 1,
+        "weights": array(
+            [
+                [1.0e-07, 1.0e-07, 1.0e-07],
+                [2.0e-08, 2.0e-08, 2.0e-08],
+                [1.0e-01, 1.0e-01, 1.0e-01],
+            ]
+        ),
+        "items": array(
+            [
+                [0.30, 0.69, 0.53],
+                [0.74, 0.14, 0.53],
+                [0.71, 0.13, 0.86],
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": array([0.07100004, 0.01300007, 0.08600006]),
+        "expect_raise": False,
+    },
+    {
+        "description": "No zero weights; nonzero order; -100 <= order <= 100; 2-d data.",
+        "order": -100,
+        "weights": array(
+            [
+                [1.0e-07, 1.0e-07, 1.0e-07],
+                [2.0e-08, 2.0e-08, 2.0e-08],
+                [1.0e-01, 1.0e-01, 1.0e-01],
+            ]
+        ),
+        "items": array(
+            [
+                [0.30, 0.69, 0.53],
+                [0.74, 0.14, 0.53],
+                [0.71, 0.13, 0.86],
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": array([0.35246927, 0.13302809, 0.62156143]),
+        "expect_raise": False,
+    },
+    {
+        "description": "Some zero weights; nonzero order; -100 <= order <= 100; 2-d data.",
+        "order": 2,
+        "weights": array(
+            [
+                [1.0e-09, 1.0e-09, 1.0e-09],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [1.0e-07, 1.0e-07, 1.0e-04],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [2.0e-08, 2.0e-08, 4.0e-08],
+                [1.0e-01, 1.0e-01, 2.0e-01],
+            ]
+        ),
+        "items": array(
+            [
+                [0.52, 0.73, 0.85],
+                [0.56, 0.28, 0.99],
+                [0.30, 0.69, 0.53],
+                [0.65, 0.65, 0.65],
+                [0.74, 0.14, 0.53],
+                [0.71, 0.13, 0.86],
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": array([0.22452176, 0.04111019, 0.3846402231]),
+        "expect_raise": False,
+    },
+    {
+        "description": "Zero weight column; nonzero order; -100 <= order <= 100; 2-d data.",
+        "order": 2,
+        "weights": array(
+            [
+                [1.0e-09, 1.0e-09, 1.0e-02],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [1.0e-10, 1.0e-10, 1.0e-02],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [2.0e-11, 2.0e-11, 2.0e-02],
+                [0.5e-08, 0.5e-08, 0.5e-01],
+            ]
+        ),
+        "items": array(
+            [
+                [0.52, 0.73, 0.85],
+                [0.56, 0.28, 0.99],
+                [0.30, 0.69, 0.53],
+                [0.65, 0.65, 0.65],
+                [0.74, 0.14, 0.53],
+                [0.71, 0.13, 0.86],
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": None,
+        "expect_raise": True,
+    },
+    {
+        "description": "No zero weights; zero order; 2-d data.",
+        "order": 0,
+        "weights": array(
+            [
+                [1.0e-07, 1.0e-07, 1.0e-07],
+                [2.0e-08, 2.0e-08, 2.0e-08],
+                [1.0e-01, 1.0e-01, 1.0e-01],
+            ]
+        ),
+        "items": array(
+            [
+                [0.30, 0.69, 0.53],
+                [0.74, 0.14, 0.53],
+                [0.71, 0.13, 0.86],
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": array([0.96633071, 0.8154443, 0.9850308]),
+        "expect_raise": False,
+    },
+    {
+        "description": "Some zero weights; zero order; 2-d data.",
+        "order": 0,
+        "weights": array(
+            [
+                [1.0e-09, 1.0e-09, 1.0e-09],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [1.0e-07, 1.0e-07, 1.0e-04],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [2.0e-08, 2.0e-08, 4.0e-08],
+                [1.0e-01, 1.0e-01, 2.0e-01],
+            ]
+        ),
+        "items": array(
+            [
+                [0.52, 0.73, 0.85],
+                [0.56, 0.28, 0.99],
+                [0.30, 0.69, 0.53],
+                [0.65, 0.65, 0.65],
+                [0.74, 0.14, 0.53],
+                [0.71, 0.13, 0.86],
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": array([0.96633071, 0.8154443, 0.9702242087]),
+        "expect_raise": False,
+    },
+    {
+        "description": "All zero weights; zero order; 2-d data.",
+        "order": 0,
+        "weights": array(
+            [
+                [1.0e-09, 1.0e-09, 1.0e-09],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [1.0e-10, 1.0e-10, 1.0e-10],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [2.0e-11, 2.0e-11, 2.0e-11],
+                [0.5e-08, 0.5e-08, 0.5e-08],
+            ]
+        ),
+        "items": array(
+            [
+                [0.52, 0.73, 0.85],
+                [0.56, 0.28, 0.99],
+                [0.30, 0.69, 0.53],
+                [0.65, 0.65, 0.65],
+                [0.74, 0.14, 0.53],
+                [0.71, 0.13, 0.86],
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": None,
+        "expect_raise": True,
+    },
+    {
+        "description": "No zero weights; order < -100; 2-d data.",
+        "order": -101,
+        "weights": array(
+            [
+                [1.0e-07, 1.0e-07, 1.0e-07],
+                [2.0e-08, 2.0e-08, 2.0e-08],
+                [1.0e-01, 1.0e-01, 1.0e-01],
+            ]
+        ),
+        "items": array(
+            [
+                [0.30, 0.69, 0.53],
+                [0.74, 0.14, 0.53],
+                [0.71, 0.13, 0.86],
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": array([0.30, 0.13, 0.53]),
+        "expect_raise": False,
+    },
+    {
+        "description": "Some zero weights; order == -inf; 2-d data.",
+        "order": -inf,
+        "weights": array(
+            [
+                [1.0e-09, 1.0e-09, 1.0e-09],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [1.0e-07, 1.0e-07, 1.0e-04],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [2.0e-08, 2.0e-08, 4.0e-08],
+                [1.0e-01, 1.0e-01, 2.0e-01],
+            ]
+        ),
+        "items": array(
+            [
+                [0.28, 0.10, 0.51],
+                [0.56, 0.28, 0.99],
+                [0.30, 0.69, 0.53],
+                [0.28, 0.09, 0.14],
+                [0.74, 0.14, 0.53],
+                [0.71, 0.13, 0.86],
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": array([0.30, 0.13, 0.53]),
+        "expect_raise": False,
+    },
+    {
+        "description": "Zero weight columns; order < -100; 2-d data.",
+        "order": -382,
+        "weights": array(
+            [
+                [1.0e-03, 1.0e-01, 1.0e-09],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [1.0e-00, 1.0e-00, 1.0e-10],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [2.0e-01, 2.0e-01, 2.0e-11],
+                [0.5e-00, 0.5e-02, 0.5e-08],
+            ]
+        ),
+        "items": array(
+            [
+                [0.52, 0.73, 0.85],
+                [0.56, 0.28, 0.99],
+                [0.30, 0.69, 0.53],
+                [0.65, 0.65, 0.65],
+                [0.74, 0.14, 0.53],
+                [0.71, 0.13, 0.86],
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": None,
+        "expect_raise": True,
+    },
+    {
+        "description": "No zero weights; order > 100; 2-d data.",
+        "order": 101,
+        "weights": array(
+            [
+                [1.0e-07, 1.0e-07, 1.0e-07],
+                [2.0e-08, 2.0e-08, 2.0e-08],
+                [1.0e-01, 1.0e-01, 1.0e-01],
+            ]
+        ),
+        "items": array(
+            [
+                [0.30, 0.69, 0.53],
+                [0.74, 0.14, 0.53],
+                [0.71, 0.13, 0.86],
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": array([0.74, 0.69, 0.86]),
+        "expect_raise": False,
+    },
+    {
+        "description": "Some zero weights; order == inf; 2-d data.",
+        "order": inf,
+        "weights": array(
+            [
+                [1.0e-09, 1.0e-09, 1.0e-09],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [1.0e-07, 1.0e-07, 1.0e-04],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [2.0e-08, 2.0e-08, 4.0e-08],
+                [1.0e-01, 1.0e-01, 2.0e-01],
+            ]
+        ),
+        "items": array(
+            [
+                [0.99, 0.79, 0.90],
+                [0.56, 0.28, 0.99],
+                [0.30, 0.69, 0.53],
+                [0.89, 0.09, 0.14],
+                [0.74, 0.14, 0.53],
+                [0.71, 0.13, 0.86],
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": array([0.74, 0.69, 0.86]),
+        "expect_raise": False,
+    },
+    {
+        "description": "All zero weights; order > 100; 2-d data.",
+        "order": 364,
+        "weights": array(
+            [
+                [1.0e-09, 1.0e-09, 1.0e-09],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [1.0e-10, 1.0e-10, 1.0e-10],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [2.0e-11, 2.0e-11, 2.0e-11],
+                [0.5e-08, 0.5e-08, 0.5e-08],
+            ]
+        ),
+        "items": array(
+            [
+                [0.52, 0.73, 0.85],
+                [0.56, 0.28, 0.99],
+                [0.30, 0.69, 0.53],
+                [0.65, 0.65, 0.65],
+                [0.74, 0.14, 0.53],
+                [0.71, 0.13, 0.86],
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": None,
+        "expect_raise": True,
+    },
+    {
+        "description": "Some zero weights; nonzero order; -100 <= order <= 100; 1-d data",
+        "order": 2,
+        "weights": array(
+            [
+                1.0e-09,
+                0.0e-00,
+                1.0e-07,
+                0.0e-00,
+                2.0e-08,
+                1.0e-01,
+            ]
+        ),
+        "items": array(
+            [
+                0.52,
+                0.56,
+                0.30,
+                0.65,
+                0.74,
+                0.71,
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": array([0.22452176]),
+        "expect_raise": False,
+    },
+    {
+        "description": "Some zero weights; zero order; 1-d data.",
+        "order": 0,
+        "weights": array(
+            [
+                1.0e-09,
+                0.0e-00,
+                1.0e-07,
+                0.0e-00,
+                2.0e-08,
+                1.0e-01,
+            ]
+        ),
+        "items": array(
+            [
+                0.52,
+                0.56,
+                0.30,
+                0.65,
+                0.74,
+                0.71,
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": array([0.96633071]),
+        "expect_raise": False,
+    },
+    {
+        "description": "Some zero weights; order == -inf; 1-d data.",
+        "order": -inf,
+        "weights": array(
+            [
+                1.0e-09,
+                0.0e-00,
+                1.0e-07,
+                0.0e-00,
+                2.0e-08,
+                1.0e-01,
+            ]
+        ),
+        "items": array(
+            [
+                0.52,
+                0.56,
+                0.30,
+                0.28,
+                0.74,
+                0.71,
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": array([0.30]),
+        "expect_raise": False,
+    },
+    {
+        "description": "Some zero weights; order == inf; 1-d data.",
+        "order": inf,
+        "weights": array(
+            [
+                1.0e-09,
+                0.0e-00,
+                1.0e-07,
+                0.0e-00,
+                2.0e-08,
+                1.0e-01,
+            ]
+        ),
+        "items": array(
+            [
+                0.52,
+                0.56,
+                0.30,
+                0.89,
+                0.74,
+                0.71,
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": array([0.74]),
+        "expect_raise": False,
+    },
+    {
+        "description": "Zero weights; nonzero order; -100 <= order <= 100; 1-d data",
+        "order": 2,
+        "weights": array(
+            [
+                1.0e-09,
+                0.0e-00,
+                1.0e-10,
+                0.0e-00,
+                0.5e-08,
+                1.0e-09,
+            ]
+        ),
+        "items": array(
+            [
+                0.52,
+                0.56,
+                0.30,
+                0.65,
+                0.74,
+                0.71,
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": None,
+        "expect_raise": True,
+    },
+    {
+        "description": "Some zero weights; nonzero order; -100 <= order <= 100; 2-d data; atol == 1e-9.",
+        "order": 2,
+        "weights": array(
+            [
+                [1.1e-09, 1.1e-09, 1.1e-09],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [1.0e-07, 1.0e-07, 1.0e-07],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [2.0e-08, 2.0e-08, 4.0e-08],
+                [1.0e-06, 1.0e-06, 2.0e-06],
+            ]
+        ),
+        "items": array(
+            [
+                [0.52, 0.73, 0.85],
+                [0.56, 0.28, 0.99],
+                [0.30, 0.69, 0.53],
+                [0.65, 0.65, 0.65],
+                [0.74, 0.14, 0.53],
+                [0.71, 0.13, 0.86],
+            ]
+        ),
+        "atol": 1e-9,
+        "expected_result": array([0.0007241198, 0.0002559066, 0.001232607298]),
+        "expect_raise": False,
+    },
+    {
+        "description": "Some zero weights; zero order; 2-d data; atol == 1e-2.",
+        "order": 0,
+        "weights": array(
+            [
+                [1.1e-02, 1.1e-02, 1.1e-02],
+                [0.0e-00, 0.0e-00, 1.0e-3],
+                [1.0e-01, 1.0e-01, 1.0e-01],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [4.0e-02, 4.0e-01, 4.0e-01],
+                [2.0e-01, 2.0e-02, 2.0e-02],
+            ]
+        ),
+        "items": array(
+            [
+                [0.52, 0.73, 0.85],
+                [0.56, 0.28, 0.99],
+                [0.30, 0.69, 0.53],
+                [0.65, 0.65, 0.65],
+                [0.74, 0.14, 0.53],
+                [0.71, 0.13, 0.86],
+            ]
+        ),
+        "atol": 1e-2,
+        "expected_result": array([0.8120992339, 0.4198668065, 0.7245218911]),
+        "expect_raise": False,
+    },
+    {
+        "description": "Some zero weights; order == -inf; 2-d data; atol == 1e-9.",
+        "order": -inf,
+        "weights": array(
+            [
+                [1.1e-09, 1.1e-09, 1.1e-09],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [1.0e-07, 1.0e-07, 1.0e-07],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [2.0e-08, 2.0e-08, 4.0e-08],
+                [1.0e-01, 1.0e-01, 2.0e-06],
+            ]
+        ),
+        "items": array(
+            [
+                [0.28, 0.10, 0.51],
+                [0.56, 0.28, 0.99],
+                [0.30, 0.69, 0.53],
+                [0.28, 0.09, 0.14],
+                [0.74, 0.14, 0.53],
+                [0.71, 0.13, 0.86],
+            ]
+        ),
+        "atol": 1e-9,
+        "expected_result": array([0.28, 0.10, 0.51]),
+        "expect_raise": False,
+    },
+    {
+        "description": "Some zero weights; order == inf; 2-d data; atol == 1e-9.",
+        "order": inf,
+        "weights": array(
+            [
+                [1.1e-09, 1.1e-09, 1.1e-09],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [1.0e-07, 1.0e-07, 1.0e-07],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [2.0e-08, 2.0e-08, 4.0e-08],
+                [1.0e-01, 1.0e-01, 2.0e-06],
+            ]
+        ),
+        "items": array(
+            [
+                [0.99, 0.79, 0.90],
+                [0.56, 0.28, 0.99],
+                [0.30, 0.69, 0.53],
+                [0.89, 0.09, 0.14],
+                [0.74, 0.14, 0.53],
+                [0.71, 0.13, 0.86],
+            ]
+        ),
+        "atol": 1e-9,
+        "expected_result": array([0.99, 0.79, 0.90]),
+        "expect_raise": False,
+    },
+    {
+        "description": "Some zero weights; nonzero order; -100 <= order <= 100; 2-d data; mismatching number of rows.",
+        "order": 2,
+        "weights": array(
+            [
+                [1.1e-09, 1.1e-09, 1.1e-09],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [1.0e-07, 1.0e-07, 1.0e-07],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [2.0e-08, 2.0e-08, 4.0e-08],
+                [1.0e-06, 1.0e-06, 2.0e-06],
+                [2.0e-05, 1.2e-02, 4.5e-07],
+            ]
+        ),
+        "items": array(
+            [
+                [0.52, 0.73, 0.85],
+                [0.56, 0.28, 0.99],
+                [0.30, 0.69, 0.53],
+                [0.65, 0.65, 0.65],
+                [0.74, 0.14, 0.53],
+                [0.71, 0.13, 0.86],
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": None,
+        "expect_raise": True,
+    },
+    {
+        "description": "Some zero weights; nonzero order; -100 <= order <= 100; 2-d data; mismatching number of columns.",
+        "order": 2,
+        "weights": array(
+            [
+                [1.1e-09, 1.1e-09, 1.1e-09],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [1.0e-07, 1.0e-07, 1.0e-07],
+                [0.0e-00, 0.0e-00, 0.0e-00],
+                [2.0e-08, 2.0e-08, 4.0e-08],
+                [1.0e-06, 1.0e-06, 2.0e-06],
+            ]
+        ),
+        "items": array(
+            [
+                [0.52, 0.73, 0.85, 1.23],
+                [0.56, 0.28, 0.99, 4.56],
+                [0.30, 0.69, 0.53, 7.89],
+                [0.65, 0.65, 0.65, 0.12],
+                [0.74, 0.14, 0.53, 3.45],
+                [0.71, 0.13, 0.86, 6.78],
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": None,
+        "expect_raise": True,
+    },
+    {
+        "description": "Some zero weights; nonzero order; -100 <= order <= 100; 1/2-d data; mismatching dimensions.",
+        "order": 2,
+        "weights": array(
+            [
+                1.1e-09,
+                0.0e-00,
+                1.0e-07,
+                0.0e-00,
+                2.0e-08,
+                1.0e-06,
+            ]
+        ),
+        "items": array(
+            [
+                [0.52, 0.73, 0.85],
+                [0.56, 0.28, 0.99],
+                [0.30, 0.69, 0.53],
+                [0.65, 0.65, 0.65],
+                [0.74, 0.14, 0.53],
+                [0.71, 0.13, 0.86],
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": None,
+        "expect_raise": True,
+    },
+    {
+        "description": "Some zero weights; nonzero order; -100 <= order <= 100; 1-d data; mismatching number of rows.",
+        "order": 2,
+        "weights": array(
+            [
+                1.1e-09,
+                0.0e-00,
+                1.0e-07,
+                0.0e-00,
+                2.0e-08,
+                1.0e-06,
+                2.0e-05,
+            ]
+        ),
+        "items": array(
+            [
+                0.52,
+                0.56,
+                0.30,
+                0.65,
+                0.74,
+                0.71,
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": None,
+        "expect_raise": True,
+    },
+    {
+        "description": "Some zero weights; nonzero order; -100 <= order <= 100; >2-d.",
+        "order": 2,
+        "weights": array(
+            [
+                [[1.1e-09, 2.2e-08], [3.3e-07, 4.4e-06]],
+                [[5.5e-06, 6.6e-05], [7.7e-04, 8.8e-03]],
+            ]
+        ),
+        "items": array(
+            [
+                [[1.1e-09, 2.2e-08], [3.3e-07, 4.4e-06]],
+                [[5.5e-06, 6.6e-05], [7.7e-04, 8.8e-03]],
+            ]
+        ),
+        "atol": 1e-8,
+        "expected_result": None,
+        "expect_raise": True,
+    },
+]
+
+
+class TestPowerMean:
+    """Tests metacommunity.utilities.power_mean."""
+
+    @mark.parametrize("test_case", POWER_MEAN_TEST_CASES)
+    def test_power_mean(self, test_case):
+        """Tests power_mean test cases."""
+        if test_case["expect_raise"]:
+            with raises(InvalidArgumentError):
+                power_mean(
+                    order=test_case["order"],
+                    weights=test_case["weights"],
+                    items=test_case["items"],
+                    atol=test_case["atol"],
+                )
+        else:
+            actual_result = power_mean(
+                order=test_case["order"],
+                weights=test_case["weights"],
+                items=test_case["items"],
+                atol=test_case["atol"],
+            )
+            assert allclose(actual_result, test_case["expected_result"])
+            assert allclose(test_case["expected_result"], actual_result)
+
+
+UNIQUE_CORRESPONDENCE_TEST_CASES = [
+    {
+        "description": "No ordering imposed; non-empty non-unique items",
+        "items": array([1, 2, 1, 5, 3, 2, 4, 100, 3, 4]),
+        "ordered_unique_items": None,
+        "n_unique": 6,
+    },
+    {
+        "description": "No ordering imposed; non-empty unique items",
+        "items": array(["foo", "bar", "zip", "zap", "zongo", "wakka"]),
+        "ordered_unique_items": None,
+        "n_unique": 6,
+    },
+    {
+        "description": "No ordering imposed; empty items",
+        "items": array([]),
+        "ordered_unique_items": None,
+        "n_unique": 0,
+    },
+    {
+        "description": "Ordering consists of uniques only; non-empty non-unique items",
+        "items": array(["foo", "bar", "fooo", "zip", "bar", "foo", "foo"]),
+        "ordered_unique_items": array(["zip", "fooo", "bar", "foo"]),
+        "n_unique": 4,
+    },
+    {
+        "description": "Ordering consists of uniques only; non-empty unique items",
+        "items": array([124, 64, 12, 32]),
+        "ordered_unique_items": array([64, 12, 32, 124]),
+        "n_unique": 4,
+    },
+    {
+        "description": "Ordering consists of uniques only; empty items",
+        "items": array([]),
+        "ordered_unique_items": array([]),
+        "n_unique": 0,
+    },
+    {
+        "description": "Ordering properly conains uniques; non-empty non-unique items",
+        "items": array(["foo", "bar", "fooo", "zip", "bar", "foo", "foo"]),
+        "ordered_unique_items": array(
+            ["zip", "shazam", "foo", "bar", "bazinga", "fooo"]
+        ),
+        "n_unique": 4,
+    },
+    {
+        "description": "Ordering properly conains uniques; non-empty unique items",
+        "items": array(["foo", "bar", "fooo", "zip"]),
+        "ordered_unique_items": array(
+            ["zip", "shazam", "foo", "bar", "bazinga", "fooo"]
+        ),
+        "n_unique": 4,
+    },
+    {
+        "description": "Ordering properly conains uniques; empty items",
+        "items": array([]),
+        "ordered_unique_items": array(
+            ["zip", "shazam", "foo", "bar", "bazinga", "fooo"]
+        ),
+        "n_unique": 0,
+    },
+]
+
+
+class TestUniqueCorrespondence:
+    """Tests diversity.utilities.unique_correspondence."""
+
+    @mark.parametrize("test_case", UNIQUE_CORRESPONDENCE_TEST_CASES)
+    def test_unique_correspondence(self, test_case):
+        """Tests unique_correspondence test cases."""
+        original_items = deepcopy(test_case["items"])
+        original_ordered_unique_items = deepcopy(test_case["ordered_unique_items"])
+        result_unique_items, result_item_positions = unique_correspondence(
+            items=test_case["items"],
+            ordered_unique_items=test_case["ordered_unique_items"],
+        )
+        assert len(set(result_unique_items)) == len(result_unique_items)
+        assert set(original_items).issubset(set(result_unique_items))
+        for unique_pos, item in zip(result_item_positions, original_items):
+            assert item == result_unique_items[unique_pos]
+        if original_ordered_unique_items is not None:
+            assert all(result_unique_items == original_ordered_unique_items)