Repository: motiwari/BanditPAM (language: python)

Issue:
New error in GHA
`--no-use-pep517 error: It is not possible to use --no-use-pep517
without setuptools and wheel installed.`
Change `exact` in `buildTarget` and `swapTarget` to bool

New error in GHA
`--no-use-pep517 error: It is not possible to use --no-use-pep517
without setuptools and wheel installed.`


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/build_linux_wheels.yml b/.github/workflows/build_linux_wheels.yml
index 439238f2..bedd0efb 100644
--- a/.github/workflows/build_linux_wheels.yml
+++ b/.github/workflows/build_linux_wheels.yml
@@ -27,8 +27,8 @@ jobs:
       - name: Install Python dependencies 
         run: |
           python -m pip install --upgrade pip
-          pip install pytest
-          pip install -r requirements.txt
+          python -m pip install pytest
+          python -m pip install -r requirements.txt
 
       - name: Build sdist
         run: python setup.py sdist --formats=gztar,zip
diff --git a/.github/workflows/build_mac_wheels.yml b/.github/workflows/build_mac_wheels.yml
index 7be1909b..f060448c 100644
--- a/.github/workflows/build_mac_wheels.yml
+++ b/.github/workflows/build_mac_wheels.yml
@@ -7,8 +7,8 @@ on:
       - published
 jobs:
   build_wheels_macos:
-    name: Build wheels on macos-10.15
-    runs-on: macos-10.15
+    name: Build wheels on macos-latest
+    runs-on: macos-latest
     strategy:
       fail-fast: false
     steps:
diff --git a/.github/workflows/run_cmake_build.yml b/.github/workflows/run_cmake_build.yml
index 09d2831b..db6bf557 100644
--- a/.github/workflows/run_cmake_build.yml
+++ b/.github/workflows/run_cmake_build.yml
@@ -38,6 +38,14 @@ jobs:
           mkdir -p data
           curl -XGET https://motiwari.com/banditpam_data/MNIST_1k.csv > data/MNIST_1k.csv
 
+      - name: Clone for PR
+        if: ${{ github.event_name == 'pull_request' }}
+        run: git clone -b $GITHUB_HEAD_REF https://github.com/motiwari/BanditPAM
+
+      - name: Clone for Push
+        if: ${{ github.event_name != 'pull_request' }}
+        run: git clone -b ${GITHUB_REF#refs/heads/} https://github.com/motiwari/BanditPAM
+
       # On the MacOS machine, when installing from source,
       # we use LLVM Clang
       # TODO(@motiwari): Clean this environment variables up
@@ -50,7 +58,6 @@ jobs:
           export DYLD_LIBRARY_PATH=/usr/local/opt/libomp/lib:/usr/local/opt/libomp/include
           export CC=/usr/local/opt/llvm/bin/clang
           export CXX=/usr/local/opt/llvm/bin/clang++
-          git clone -b ${GITHUB_REF#refs/heads/} https://github.com/motiwari/BanditPAM
           cd BanditPAM
           mkdir build
           cd build
diff --git a/.github/workflows/run_style_checks.yml b/.github/workflows/run_style_checks.yml
index 7737feeb..f2411375 100644
--- a/.github/workflows/run_style_checks.yml
+++ b/.github/workflows/run_style_checks.yml
@@ -22,8 +22,8 @@ jobs:
       - name: Install Python dependencies
         run: |
           python -m pip install --upgrade pip
-          pip install flake8 pytest cpplint
-          pip install -r requirements.txt
+          python -m pip install flake8 pytest cpplint
+          python -m pip install -r requirements.txt
 
       - name: Lint with flake8
         run: |
diff --git a/README.md b/README.md
index 40ee0926..4486c02b 100644
--- a/README.md
+++ b/README.md
@@ -24,7 +24,7 @@ Mo Tiwari, Martin Jinye Zhang, James Mayclin, Sebastian Thrun, Chris Piech, Ilan
 ```
 
 # Requirements
-# TL;DR run `pip install banditpam` or `install.packages(banditpam)` and jump to the [examples](https://github.com/motiwari/BanditPAM#example-1-synthetic-data-from-a-gaussian-mixture-model). 
+# TL;DR run `python -m pip install banditpam` or `install.packages(banditpam)` and jump to the [examples](https://github.com/motiwari/BanditPAM#example-1-synthetic-data-from-a-gaussian-mixture-model). 
 
 If you have any difficulties, please see the [platform-specific guides](https://github.com/motiwari/BanditPAM#platform-specific-installation-guides) and file a Github issue if you have additional trouble.
 
@@ -41,8 +41,8 @@ If you have any difficulties, please see the [platform-specific guides](https://
 ## Install the repo and its dependencies:
 This can be done either through PyPI (recommended)
 ```python
-/BanditPAM/: pip install -r requirements.txt
-/BanditPAM/: pip install banditpam
+/BanditPAM/: python -m pip install -r requirements.txt
+/BanditPAM/: python -m pip install banditpam
 ```
 OR through the source code via
 ```python
@@ -50,8 +50,8 @@ OR through the source code via
 /BanditPAM/: cd headers/carma
 /BanditPAM/: mkdir build && cd build && cmake -DCARMA_INSTALL_LIB=ON .. && sudo cmake --build . --config Release --target install
 /BanditPAM/: cd ../../..
-/BanditPAM/: pip install -r requirements.txt
-/BanditPAM/: sudo pip install .
+/BanditPAM/: python -m pip install -r requirements.txt
+/BanditPAM/: sudo python -m pip install .
 ```
 
 ### Example 1: Synthetic data from a Gaussian Mixture Model
@@ -201,7 +201,7 @@ The package currently supports a number of distance metrics, including all $L_p$
 
 If you're willing to write a little C++, you only need to add a few lines to [kmedoids_algorithm.cpp](https://github.com/motiwari/BanditPAM/blob/main/src/kmedoids_algorithm.cpp#L560-L615) and [kmedoids_algorithm.hpp](https://github.com/motiwari/BanditPAM/blob/main/headers/kmedoids_algorithm.hpp#L136-L142) to implement your distance metric / pairwise dissimilarity!
 
-Then, be sure to re-install the repository with a `pip install .` (note the trailing `.`).
+Then, be sure to re-install the repository with a `python -m pip install .` (note the trailing `.`).
 
 The maintainers of this repository are working on permitting arbitrary dissimilarity metrics that users write in Python, as well; see [#4](https://github.com/motiwari/BanditPAM/issues/4).
 
diff --git a/docs/install_linux.md b/docs/install_linux.md
index e7e656cb..62af6fd9 100644
--- a/docs/install_linux.md
+++ b/docs/install_linux.md
@@ -12,7 +12,7 @@ Please ensure the following dependencies are installed:
  - CARMA: via the instructions in [its guide](https://github.com/RUrlus/carma#installation)
  - Python3: if not installed, we recommend installing Python3 via [Anaconda](https://www.anaconda.com/products/individual), which is CPython compiled with `clang`
  - `pip` for your Python3 installation: this should be completed if installing via Anaconda above
- - The necessary python packages: via `pip install -r requirements.txt`
+ - The necessary python packages: via `python -m pip install -r requirements.txt`
  
  (NOT RECOMMENDED): Instead of LLVM's `clang`, you can also use another C++ compiler and point your `CC` environment variable to it. Please ensure this is the same compiler used to compile your Python installation if using CPython. If you open a `python` REPL it will show the compiler used during the language installation:
 
@@ -30,8 +30,8 @@ Do not attempt to compile the BanditPAM extension with `gcc` if your Python used
 ## BanditPAM Installation
 
 BanditPAM can then be installed via one of the following ways:
-1) Running `pip install banditpam`, OR
-2) Running `pip install .` in the home directory (`/BanditPAM`)
+1) Running `python -m pip install banditpam`, OR
+2) Running `python -m pip install .` in the home directory (`/BanditPAM`)
 
 ## Known Issues 
 The following is a list of issues seen when installing BanditPAM on Linux. To report a bug, please file an issue at https://github.com/motiwari/BanditPAM/
diff --git a/docs/install_mac.md b/docs/install_mac.md
index d3da97bb..04c56e60 100644
--- a/docs/install_mac.md
+++ b/docs/install_mac.md
@@ -11,7 +11,7 @@ Please ensure the following dependencies are installed:
  - CARMA: via the instructions in [its guide](https://github.com/RUrlus/carma#installation)
  - Python3: if not installed, we recommend installing Python3 via [Anaconda](https://www.anaconda.com/products/individual), which is CPython compiled with `clang`
  - `pip` for your Python3 installation; this should be completed if installing via Anaconda above
- - The necessary python packages: via `pip install -r requirements.txt`
+ - The necessary python packages: via `python -m pip install -r requirements.txt`
 
 (NOT RECOMMENDED): Instead of LLVM's `clang`, you can also use another C++ compiler and point your `CC` environment variable to it. Please ensure this is the same compiler used to compile your Python installation if using CPython. If you open a `python` REPL it will show the compiler used during the language installation:
 
@@ -29,8 +29,8 @@ For the default Python2 installed on Mac, the `Apple clang` compiler is used; fo
 ## BanditPAM Installation
 
 BanditPAM can then be installed via one of the following ways:
-1) Running `pip install banditpam`, OR
-2) Running `pip install .` in the home directory (`/BanditPAM`)
+1) Running `python -m pip install banditpam`, OR
+2) Running `python -m pip install .` in the home directory (`/BanditPAM`)
 
 ## Known Issues 
 The following is a list of issues seen when installing BanditPAM on MacOS. To report a bug, please file an issue at https://github.com/motiwari/BanditPAM/
diff --git a/docs/install_windows.md b/docs/install_windows.md
index c446935b..6822cf1c 100644
--- a/docs/install_windows.md
+++ b/docs/install_windows.md
@@ -1,38 +1,30 @@
 # Installation Tutorial for Windows
 
-The following is a description of the installation process of BanditPAM for Windows. This assumes that:
+The following is a more detailed description of the installation process of BanditPAM for Windows.
  
 ## Prerequisites
 Please ensure the following dependencies are installed:
- - A C++ compiler; we recommend LLVM's `clang`: via the [LLVM installation instructions](https://clang.llvm.org/get_started.html)
- - OpenMP: if using LLVM's `clang`, then OpenMP is already enabled
- - `CMake`: via the [CMake installation instructions](https://cmake.org/install/)
- - Armadillo: via the [Armadillo installation instructions](http://arma.sourceforge.net/download.html)
+ - Visual Studio (2022 is recommended) with "Desktop development with C++" workload: via the [Visual Studio download page](https://visualstudio.microsoft.com/vs/)
+ - The Armadillo library: via `git clone https://gitlab.com/conradsnicta/armadillo-code.git armadillo` in the `headers` directory
  - CARMA: via the instructions in [its guide](https://github.com/RUrlus/carma#installation)
- - Python3: if not installed, we recommend installing Python3 via [Anaconda](https://www.anaconda.com/products/individual), which is CPython compiled with `clang`
+ - Python3: if not installed, we recommend installing Python3 via [Anaconda](https://www.anaconda.com/products/individual)
  - `pip` for your Python3 installation: this should be completed if installing via Anaconda above
- - The necessary python packages: via `pip install -r requirements.txt`
- 
- (NOT RECOMMENDED): Instead of LLVM's `clang`, you can also use another C++ compiler and point your `CC` environment variable to it. Please ensure this is the same compiler used to compile your Python installation if using CPython. If you open a `python` REPL it will show the compiler used during the language installation:
-
- ```
- >> python
-Python 3.7.9 (default, Mar  1 2021, 13:32:26)
-[Clang 11.0.0 (clang-1100.0.33.17)] :: Intel Corporation on darwin
-Type "help", "copyright", "credits" or "license" for more information.
-Intel(R) Distribution for Python is brought to you by Intel Corporation.
-Please check out: https://software.intel.com/en-us/python-distribution
-```
-
-Do not attempt to compile the BanditPAM extension with `gcc` if your Python used `clang`, and vice versa.
+ - The necessary python packages: via `python -m pip install -r requirements.txt`
 
 ## BanditPAM Installation
 
-BanditPAM can then be installed via one of the following ways:
-1) Running `pip install banditpam`, OR
-2) Running `pip install .` in the home directory (`/BanditPAM`)
+CMake Build:
 
-## Known Issues 
-The following is a list of issues seen when installing BanditPAM on Windows. To report a bug, please file an issue at https://github.com/motiwari/BanditPAM/
+1) Run `scripts/retrieve_windows_cmake_files.sh` to retrieve the files necessary for the Windows CMake build
+2) Add the Visual Studio IDE location to PATH in Environment Variables (e.g. `C:\Program Files\Microsoft Visual Studio\2022\Community\Common7\IDE`).
+3) Add `${project_dir}\headers\armadillo\examples\lib_win64` to PATH 
+4) Run `devenv BanditPAM.sln /Build "Release|x64"` after `cmake ..` in a prompt other than Git Bash
+5) The `.exe` will be located in `build/src/Release`
 
-- You may run into an error like `error: ‘class arma::Mat<float>’ has no member named ‘n_alloc’`; this is fundamentally due to an [error in CARMA](https://github.com/RUrlus/carma/pull/98). If this isn't fixed by the time you're reading this, follow the fix described in [this issue](https://github.com/motiwari/BanditPAM/issues/169).
+Python Build:
+1) Run `python -m pip install banditpam`, OR
+2) Follow these steps:
+   1) Add the location of `cl.exe` to PATH in Environment Variables (e.g. `C:\Program Files\Microsoft Visual Studio\2022\Community\VC\Tools\MSVC\14.34.31933\bin\Hostx64\x64`).
+   2) Run `python -m pip install .` in the home directory (`/BanditPAM`)
+   2) Add the file `C:\Program Files\Microsoft Visual Studio\2022\Community\VC\Tools\MSVC\14.34.31933\bin\Hostx86\x64\clang_rt.asan_dynamic-x86_64.dll` to `build\lib.win-amd64-cpython-310`
+   3) Run `python -m pip install .` in the home directory (`/BanditPAM`)
\ No newline at end of file
diff --git a/headers/algorithms/banditpam.hpp b/headers/algorithms/banditpam.hpp
index 6489de89..a9915703 100644
--- a/headers/algorithms/banditpam.hpp
+++ b/headers/algorithms/banditpam.hpp
@@ -22,8 +22,8 @@ class BanditPAM : public km::KMedoids {
    * @param inputData Input data to cluster
    */
   void fitBanditPAM(
-    const arma::fmat& inputData,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat);
+          const arma::fmat &inputData,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat);
 
   /**
    * @brief Empirical estimation of standard deviation of arm returns
@@ -37,10 +37,10 @@ class BanditPAM : public km::KMedoids {
    * @returns Estimate of each arm's standard deviation
    */
   arma::frowvec buildSigma(
-    const arma::fmat& data,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-    const arma::frowvec& bestDistances,
-    const bool useAbsolute);
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const arma::frowvec &bestDistances,
+          const bool useAbsolute);
 
   /**
    * @brief Estimates the mean reward for each arm in the BUILD steps.
@@ -50,17 +50,17 @@ class BanditPAM : public km::KMedoids {
    * @param bestDistances Contains best distances from each point to medoids
    * @param useAbsolute Flag to use the absolute distance to each arm instead
    * of improvement over prior loss; necessary for the first BUILD step
-   * @param exact 0 if using standard batch size; size of dataset otherwise
+   * @param exact false if using standard batch size; true otherwise
    *
    * @returns Estimate of each arm's change in loss
    */
   arma::frowvec buildTarget(
-    const arma::fmat& data,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-    const arma::uvec* target,
-    const arma::frowvec* bestDistances,
-    const bool useAbsolute,
-    const size_t exact);
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const arma::uvec *target,
+          const arma::frowvec *bestDistances,
+          const bool useAbsolute,
+          const bool exact);
 
   /**
    * @brief Performs the BUILD step of BanditPAM.
@@ -77,10 +77,10 @@ class BanditPAM : public km::KMedoids {
    * @param medoids Matrix that contains the coordinates of each medoid
    */
   void build(
-    const arma::fmat& data,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-    arma::urowvec* medoidIndices,
-    arma::fmat* medoids);
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          arma::urowvec *medoidIndices,
+          arma::fmat *medoids);
 
   /**
    * @brief Empirical estimation of standard deviation of arm returns
@@ -94,11 +94,11 @@ class BanditPAM : public km::KMedoids {
    * @returns Estimate of each arm's standard deviation
    */
   arma::fmat swapSigma(
-    const arma::fmat& data,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-    const arma::frowvec* bestDistances,
-    const arma::frowvec* secondBestDistances,
-    const arma::urowvec* assignments);
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const arma::frowvec *bestDistances,
+          const arma::frowvec *secondBestDistances,
+          const arma::urowvec *assignments);
 
   /**
    * @brief Estimates the mean reward for each arm in SWAP step.
@@ -110,19 +110,19 @@ class BanditPAM : public km::KMedoids {
    * @param bestDistances Contains best distances from each point to medoids
    * @param secondBestDistances Contains second best distances from each point to medoids
    * @param assignments Assignments of datapoints to their closest medoid
-   * @param exact 0 if using standard batch size; size of dataset otherwise
+   * @param exact false if using standard batch size; true otherwise
    *
    * @returns Estimate of each arm's change in loss
    */
   arma::fmat swapTarget(
-    const arma::fmat& data,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-    const arma::urowvec* medoidIndices,
-    const arma::uvec* targets,
-    const arma::frowvec* bestDistances,
-    const arma::frowvec* secondBestDistances,
-    const arma::urowvec* assignments,
-    const size_t exact);
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const arma::urowvec *medoidIndices,
+          const arma::uvec *targets,
+          const arma::frowvec *bestDistances,
+          const arma::frowvec *secondBestDistances,
+          const arma::urowvec *assignments,
+          const bool exact);
 
   /**
   * @brief Performs the SWAP step of BanditPAM.
@@ -141,11 +141,11 @@ class BanditPAM : public km::KMedoids {
   * @param assignments Array of containing the medoid each point is closest to
   */
   void swap(
-    const arma::fmat& data,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-    arma::urowvec* medoidIndices,
-    arma::fmat* medoids,
-    arma::urowvec* assignments);
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          arma::urowvec *medoidIndices,
+          arma::fmat *medoids,
+          arma::urowvec *assignments);
 };
 }  // namespace km
 #endif  // HEADERS_ALGORITHMS_BANDITPAM_HPP_
diff --git a/headers/algorithms/banditpam_orig.hpp b/headers/algorithms/banditpam_orig.hpp
index 5165cdde..265e1be2 100644
--- a/headers/algorithms/banditpam_orig.hpp
+++ b/headers/algorithms/banditpam_orig.hpp
@@ -22,8 +22,8 @@ class BanditPAM_orig : public km::KMedoids {
    * @param inputData Input data to cluster
    */
   void fitBanditPAM_orig(
-      const arma::fmat& inputData,
-      std::optional<std::reference_wrapper<const arma::fmat>> distMat);
+          const arma::fmat &inputData,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat);
 
   /**
    * @brief Empirical estimation of standard deviation of arm returns
@@ -37,10 +37,10 @@ class BanditPAM_orig : public km::KMedoids {
    * @returns Estimate of each arm's standard deviation
    */
   arma::frowvec buildSigma(
-    const arma::fmat& data,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-    const arma::frowvec& bestDistances,
-    const bool useAbsolute);
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const arma::frowvec &bestDistances,
+          const bool useAbsolute);
 
   /**
    * @brief Estimates the mean reward for each arm in the BUILD steps.
@@ -50,17 +50,17 @@ class BanditPAM_orig : public km::KMedoids {
    * @param bestDistances Contains best distances from each point to medoids
    * @param useAbsolute Flag to use the absolute distance to each arm instead
    * of improvement over prior loss; necessary for the first BUILD step
-   * @param exact 0 if using standard batch size; size of dataset otherwise
+   * @param exact false if using standard batch size; true otherwise
    *
    * @returns Estimate of each arm's change in loss
    */
   arma::frowvec buildTarget(
-    const arma::fmat& data,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-    const arma::uvec* target,
-    const arma::frowvec* bestDistances,
-    const bool useAbsolute,
-    const size_t exact);
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const arma::uvec *target,
+          const arma::frowvec *bestDistances,
+          const bool useAbsolute,
+          const bool exact);
 
   /**
    * @brief Performs the BUILD step of BanditPAM_orig.
@@ -77,10 +77,10 @@ class BanditPAM_orig : public km::KMedoids {
    * @param medoids Matrix that contains the coordinates of each medoid
    */
   void build(
-    const arma::fmat& data,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-    arma::urowvec* medoidIndices,
-    arma::fmat* medoids);
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          arma::urowvec *medoidIndices,
+          arma::fmat *medoids);
 
   /**
    * @brief Empirical estimation of standard deviation of arm returns
@@ -94,11 +94,11 @@ class BanditPAM_orig : public km::KMedoids {
    * @returns Estimate of each arm's standard deviation
    */
   arma::fmat swapSigma(
-    const arma::fmat& data,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-    const arma::frowvec* bestDistances,
-    const arma::frowvec* secondBestDistances,
-    const arma::urowvec* assignments);
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const arma::frowvec *bestDistances,
+          const arma::frowvec *secondBestDistances,
+          const arma::urowvec *assignments);
 
   /**
    * @brief Estimates the mean reward for each arm in SWAP step.
@@ -111,19 +111,19 @@ class BanditPAM_orig : public km::KMedoids {
    * @param secondBestDistances Contains second best distances
    * from each point to medoids
    * @param assignments Assignments of datapoints to their closest medoid
-   * @param exact 0 if using standard batch size; size of dataset otherwise
+   * @param exact false if using standard batch size; true otherwise
    *
    * @returns Estimate of each arm's change in loss
    */
   arma::fvec swapTarget(
-    const arma::fmat& data,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-    const arma::urowvec* medoidIndices,
-    const arma::uvec* targets,
-    const arma::frowvec* bestDistances,
-    const arma::frowvec* secondBestDistances,
-    const arma::urowvec* assignments,
-    const size_t exact);
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const arma::urowvec *medoidIndices,
+          const arma::uvec *targets,
+          const arma::frowvec *bestDistances,
+          const arma::frowvec *secondBestDistances,
+          const arma::urowvec *assignments,
+          const bool exact);
 
   /**
   * @brief Performs the SWAP step of BanditPAM_orig.
@@ -142,11 +142,11 @@ class BanditPAM_orig : public km::KMedoids {
   * @param assignments Array of containing the medoid each point is closest to
   */
   void swap(
-    const arma::fmat& data,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-    arma::urowvec* medoidIndices,
-    arma::fmat* medoids,
-    arma::urowvec* assignments);
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          arma::urowvec *medoidIndices,
+          arma::fmat *medoids,
+          arma::urowvec *assignments);
 };
 }  // namespace km
 #endif  // HEADERS_ALGORITHMS_BANDITPAM_ORIG_HPP_
diff --git a/headers/algorithms/fastpam1.hpp b/headers/algorithms/fastpam1.hpp
index a633196e..d2469e1d 100644
--- a/headers/algorithms/fastpam1.hpp
+++ b/headers/algorithms/fastpam1.hpp
@@ -22,12 +22,12 @@ class FastPAM1 : public km::KMedoids {
    * @param inputData Input data to cluster
    */
   void fitFastPAM1(
-    const arma::fmat& inputData,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat);
+          const arma::fmat &inputData,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat);
 
   /**
    * @brief Performs the BUILD step of FastPAM1.
-   * 
+   *
    * Loops over all datapoint and checks each's distance to every other
    * datapoint in the dataset, then adds the point with the lowest overall
    * loss to the set of medoids.
@@ -37,16 +37,16 @@ class FastPAM1 : public km::KMedoids {
    * as medoids are identified
    */
   void buildFastPAM1(
-    const arma::fmat& data,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-    arma::urowvec* medoidIndices);
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          arma::urowvec *medoidIndices);
 
   /**
    * @brief Performs the SWAP step of FastPAM1.
    *
    * Loops over all (medoid, non-medoid) pairs and computes the change in loss
    * when the points are swapped in and out of the medoid set. Then updates
-   * the list of medoids by performing the swap that would lower the overall 
+   * the list of medoids by performing the swap that would lower the overall
    * loss the most, provided at least one such swap would reduce the loss.
    *
    * @param data Transposed input data to cluster
@@ -55,10 +55,10 @@ class FastPAM1 : public km::KMedoids {
    * @param assignments Array of containing the medoid each point is closest to
    */
   void swapFastPAM1(
-    const arma::fmat& data,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-    arma::urowvec* medoidIndices,
-    arma::urowvec* assignments);
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          arma::urowvec *medoidIndices,
+          arma::urowvec *assignments);
 };
 }  // namespace km
 #endif  // HEADERS_ALGORITHMS_FASTPAM1_HPP_
diff --git a/headers/algorithms/kmedoids_algorithm.hpp b/headers/algorithms/kmedoids_algorithm.hpp
index 9a412bf5..687b52d1 100644
--- a/headers/algorithms/kmedoids_algorithm.hpp
+++ b/headers/algorithms/kmedoids_algorithm.hpp
@@ -29,16 +29,16 @@ class KMedoids {
   //  arguments in kmedoids_pywrapper.cpp, otherwise undefined behavior can
   //  result (variables being initialized with others' values)
   KMedoids(
-    size_t nMedoids = 5,
-    const std::string& algorithm = "BanditPAM",
-    size_t maxIter = 100,
-    size_t buildConfidence = 3,
-    size_t swapConfidence = 5,
-    bool useCache = true,
-    bool usePerm = true,
-    size_t cacheWidth = 1000,
-    bool parallelize = true,
-    size_t seed = 0);
+          size_t nMedoids = 5,
+          const std::string &algorithm = "BanditPAM",
+          size_t maxIter = 100,
+          size_t buildConfidence = 3,
+          size_t swapConfidence = 5,
+          bool useCache = true,
+          bool usePerm = true,
+          size_t cacheWidth = 1000,
+          bool parallelize = true,
+          size_t seed = 0);
 
   ~KMedoids();
 
@@ -47,59 +47,59 @@ class KMedoids {
    *
    * @param inputData Input data to cluster
    * @param loss The loss function used during medoid computation
-   * 
+   *
    * @throws if the input data is empty.
    */
   void fit(
-    const arma::fmat& inputData,
-    const std::string& loss,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat);
+          const arma::fmat &inputData,
+          const std::string &loss,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat);
 
   /**
    * @brief Returns the medoids at the end of the BUILD step.
-   * 
+   *
    * @returns Medoids at the end of the BUILD step
    */
   arma::urowvec getMedoidsBuild() const;
 
   /**
    * @brief Returns the medoids after all SWAP steps have been run.
-   * 
+   *
    * @returns Medoids at the end of all SWAP steps
    */
   arma::urowvec getMedoidsFinal() const;
 
   /**
    * @brief Returns the medoid assignments for each datapoint.
-   * 
+   *
    * @returns Cluster assignment for each point
    */
   arma::urowvec getLabels() const;
 
   /**
    * @brief Returns the number of SWAP steps performed.
-   * 
+   *
    * @returns Number of swap steps performed
    */
   size_t getSteps() const;
 
   /**
    * @brief Returns the number of medoids, k.
-   * 
+   *
    * @returns Current value of k, the number of medoids/clusters
    */
   size_t getNMedoids() const;
 
   /**
    * @brief Sets the number of medoids, k.
-   * 
+   *
    * @param newNMedoids The new number, k, of medoids to use
    */
   void setNMedoids(size_t newNMedoids);
 
   /**
    * @brief Returns the algorithm being used for k-medoids clustering.
-   * 
+   *
    * @returns The name of the algorithm being used
    */
   std::string getAlgorithm() const;
@@ -109,11 +109,11 @@ class KMedoids {
    *
    * @param newAlgorithm The new algorithm to use: "BanditPAM", "PAM", or "FastPAM1"
    */
-  void setAlgorithm(const std::string& newAlgorithm);
+  void setAlgorithm(const std::string &newAlgorithm);
 
   /**
    * @brief Returns the maximum number of SWAP steps during clustering.
-   * 
+   *
    * @returns The maximum number of SWAP steps that can be
    * performed during clustering
    */
@@ -121,7 +121,7 @@ class KMedoids {
 
   /**
    * @brief Sets the maximum number of SWAP steps during clustering.
-   * 
+   *
    * @param newMaxIter New maximum number of iterations to use
    */
   void setMaxIter(size_t newMaxIter);
@@ -129,7 +129,7 @@ class KMedoids {
   /**
    * @brief Returns the buildConfidence, a parameter that affects the width
    * of the confidence intervals during the BUILD step.
-   * 
+   *
    * @returns The current value of buildConfidence
    */
   size_t getBuildConfidence() const;
@@ -139,7 +139,7 @@ class KMedoids {
    * of the confidence intervals during the BUILD step.
    *
    * @param newBuildConfidence The new buildConfidence to use
-   * 
+   *
    * @throws If attempting to set buildConfidence when not using BanditPAM
    */
   void setBuildConfidence(size_t newBuildConfidence);
@@ -147,7 +147,7 @@ class KMedoids {
   /**
    * @brief Returns the swapConfidence, a parameter that affects the width
    * of the confidence intervals during the BUILD step.
-   * 
+   *
    * @returns The current value of swapConfidence
    */
   size_t getSwapConfidence() const;
@@ -157,7 +157,7 @@ class KMedoids {
    * of the confidence intervals during the BUILD step.
    *
    * @param newSwapConfidence The new swapConfidence to use
-   * 
+   *
    * @throws If attempting to set buildConfidence when not using BanditPAM
    */
   void setSwapConfidence(size_t newSwapConfidence);
@@ -180,7 +180,7 @@ class KMedoids {
    * @brief Sets the loss function to use during clustering.
    *
    * @param loss Loss function to be used e.g. "L2"
-   * 
+   *
    * @throws If the loss function is unrecognized
    */
   void setLossFn(std::string loss);
@@ -201,6 +201,15 @@ class KMedoids {
    */
   float getAverageLoss() const;
 
+  /**
+   * @brief Get the loss at the end of the BUILD step of the algorithm
+   *
+   * @returns The loss at the end of the BUILD step of the algorithm
+   *
+   * @throws If the BUILD step has not been run yet
+   */
+  float getBuildLoss() const;
+
   /**
    * @brief Returns whether a distance cache is being used
    *
@@ -321,7 +330,7 @@ class KMedoids {
   float getTimePerSwap() const;
 
   /// The cache which stores pairwise distance computations
-  float* cache;
+  float *cache;
 
   /// The permutation in which to sample the reference points
   arma::uvec permutation;
@@ -350,118 +359,118 @@ class KMedoids {
    * @param assignments Assignments of datapoints to their closest medoid
    */
   void calcBestDistancesSwap(
-    const arma::fmat& data,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-    const arma::urowvec* medoidIndices,
-    arma::frowvec* bestDistances,
-    arma::frowvec* secondBestDistances,
-    arma::urowvec* assignments,
-    const bool swapPerformed = true);
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const arma::urowvec *medoidIndices,
+          arma::frowvec *bestDistances,
+          arma::frowvec *secondBestDistances,
+          arma::urowvec *assignments,
+          const bool swapPerformed = true);
 
   /**
-   * @brief Calculate the average loss for the given choice of medoids. 
-   * 
+   * @brief Calculate the average loss for the given choice of medoids.
+   *
    * @param data Transposed data to cluster
    * @param medoidIndices Indices of the medoids in the dataset
-   * 
+   *
    * @returns The average loss, i.e., the average distance from each point to its
    * nearest medoid
    */
   float calcLoss(
-    const arma::fmat& data,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-    const arma::urowvec* medoidIndices);
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const arma::urowvec *medoidIndices);
 
   /**
    * @brief A wrapper around the given loss function that caches distances
    * between the given points.
-   * 
+   *
    * NOTE: if you change useCacheFunctionOverride, also change useCache
-   * 
+   *
    * @param data Transposed data to cluster
    * @param i Index of first datapoint
    * @param j Index of second datapoint
    * @param useCacheFunctionOverride Whether to use the cache in this function (by default, uses value of useCache)
-   * 
+   *
    * @returns The distance between points i and j
    */
   float cachedLoss(
-    const arma::fmat& data,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-    const size_t i,
-    const size_t j,
-    const size_t category,
-    const bool useCacheFunctionOverride = true);
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const size_t i,
+          const size_t j,
+          const size_t category,
+          const bool useCacheFunctionOverride = true);
 
   /// If using an L_p loss, the value of p
   size_t lp;
 
   /**
-   * @brief Computes the Lp distance between the 
+   * @brief Computes the Lp distance between the
    * datapoints of indices i and j in the dataset
    *
    * @param data Transposed data to cluster
    * @param i Index of first datapoint
    * @param j Index of second datapoint
-   * 
+   *
    * @returns The Lp distance between points i and j
    */
   float LP(
-    const arma::fmat& data,
-    const size_t i,
-    const size_t j) const;
+          const arma::fmat &data,
+          const size_t i,
+          const size_t j) const;
 
   /**
-   * @brief Computes the L-infinity distance between the 
+   * @brief Computes the L-infinity distance between the
    * datapoints of indices i and j in the dataset
    *
    * @param data Transposed data to cluster
    * @param i Index of first datapoint
    * @param j Index of second datapoint
-   * 
+   *
    * @returns The L-infinity distance between points i and j
    */
-  float LINF(const arma::fmat& data,
-    const size_t i,
-    const size_t j) const;
+  float LINF(const arma::fmat &data,
+             const size_t i,
+             const size_t j) const;
 
   /**
-   * @brief Computes the cosine distance between the 
+   * @brief Computes the cosine distance between the
    * datapoints of indices i and j in the dataset
    *
    * @param data Transposed data to cluster
    * @param i Index of first datapoint
    * @param j Index of second datapoint
-   * 
+   *
    * @returns The cosine distance between points i and j
    */
-  float cos(const arma::fmat& data,
-    const size_t i,
-    const size_t j) const;
+  float cos(const arma::fmat &data,
+            const size_t i,
+            const size_t j) const;
 
   /**
-   * @brief Computes the Manhattan distance between the 
+   * @brief Computes the Manhattan distance between the
    * datapoints of indices i and j in the dataset
    *
    * @param data Transposed data to cluster
    * @param i Index of first datapoint
    * @param j Index of second datapoint
-   * 
+   *
    * @returns The Manhattan distance between points i and j
    */
-  float manhattan(const arma::fmat& data,
-    const size_t i,
-    const size_t j) const;
+  float manhattan(const arma::fmat &data,
+                  const size_t i,
+                  const size_t j) const;
 
   /**
-   * @brief Checks whether algorithm choice is valid. The given 
-   * algorithm must be either "BanditPAM", "PAM", or "FastPAM1". 
+   * @brief Checks whether algorithm choice is valid. The given
+   * algorithm must be either "BanditPAM", "PAM", or "FastPAM1".
    *
    * @param algorithm Name of the k-medoids algorithm to use
-   * 
+   *
    * @throws If the algorithm is invalid.
    */
-  void checkAlgorithm(const std::string& algorithm) const;
+  void checkAlgorithm(const std::string &algorithm) const;
 
   /// Number of medoids to use -- the "k" in k-medoids
   size_t nMedoids;
@@ -486,10 +495,10 @@ class KMedoids {
 
   /// Function pointer to the loss function to use
   float (KMedoids::*lossFn)(
-    const arma::fmat& data,
-    const size_t i,
-    const size_t j)
-    const;
+          const arma::fmat &data,
+          const size_t i,
+          const size_t j)
+  const;
 
   /// Number of SWAP steps performed
   size_t steps = 0;
@@ -521,6 +530,9 @@ class KMedoids {
   /// Contains the average loss at the last step of the algorithm
   float averageLoss = 0.0;
 
+  /// Contains the loss at the end of the BUILD step of the algorithm
+  float buildLoss = 0.0;
+
   /// Number of points to sample per reference batch
   size_t batchSize = 100;
 
diff --git a/headers/algorithms/pam.hpp b/headers/algorithms/pam.hpp
index 1b5340e4..61d417d1 100644
--- a/headers/algorithms/pam.hpp
+++ b/headers/algorithms/pam.hpp
@@ -22,8 +22,8 @@ class PAM : public km::KMedoids {
   * @param inputData Input data to cluster
   */
   void fitPAM(
-  const arma::fmat& inputData,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat);
+          const arma::fmat &inputData,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat);
 
   /**
   * @brief Performs the BUILD step of PAM.
@@ -37,9 +37,9 @@ class PAM : public km::KMedoids {
   * as medoids are identified
   */
   void buildPAM(
-    const arma::fmat& data,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-    arma::urowvec* medoidIndices);
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          arma::urowvec *medoidIndices);
 
   /**
   * @brief Performs the SWAP steps of BanditPAM.
@@ -55,10 +55,10 @@ class PAM : public km::KMedoids {
   * @param assignments Array of containing the medoid each point is closest to
   */
   void swapPAM(
-    const arma::fmat& data,
-    std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-    arma::urowvec* medoidIndices,
-    arma::urowvec* assignments);
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          arma::urowvec *medoidIndices,
+          arma::urowvec *assignments);
 };
 }  // namespace km
 #endif  // HEADERS_ALGORITHMS_PAM_HPP_
diff --git a/headers/carma b/headers/carma
index 75057ddd..2fbc2e6f 160000
--- a/headers/carma
+++ b/headers/carma
@@ -1,1 +1,1 @@
-Subproject commit 75057dddb7d6f7422db70a5d937709e10fb25dc5
+Subproject commit 2fbc2e6faf2e40e41003c06cbb13744405732b5f
diff --git a/headers/python_bindings/kmedoids_pywrapper.hpp b/headers/python_bindings/kmedoids_pywrapper.hpp
index ede77534..1102a825 100644
--- a/headers/python_bindings/kmedoids_pywrapper.hpp
+++ b/headers/python_bindings/kmedoids_pywrapper.hpp
@@ -29,9 +29,9 @@ class KMedoidsWrapper : public km::KMedoids {
    * @param k The number of medoids to compute
    */
   void fitPython(
-    const pybind11::array_t<float>& inputData,
-    const std::string& loss,
-    pybind11::kwargs kw);
+          const pybind11::array_t<float> &inputData,
+          const std::string &loss,
+          pybind11::kwargs kw);
 
   /**
    * @brief Returns the build medoids
@@ -39,7 +39,7 @@ class KMedoidsWrapper : public km::KMedoids {
    * Returns as a numpy array the build medoids at the end of the BUILD step
    * after KMedoids::fit has been called.
    */
-  pybind11::array_t<arma::uword> getMedoidsBuildPython();
+  pybind11::array_t <arma::uword> getMedoidsBuildPython();
 
   /**
    * @brief Returns the final medoids
@@ -47,7 +47,7 @@ class KMedoidsWrapper : public km::KMedoids {
    * Returns as a numpy array the final medoids at the end of the SWAP step
    * after KMedoids::fit has been called.
    */
-  pybind11::array_t<arma::uword> getMedoidsFinalPython();
+  pybind11::array_t <arma::uword> getMedoidsFinalPython();
 
   /**
    * @brief Returns the medoid assignments for each datapoint
@@ -55,7 +55,7 @@ class KMedoidsWrapper : public km::KMedoids {
    * Returns as a numpy array the medoid each input datapoint is assigned to
    * after KMedoids::fit is called and the final medoids have been identified
    */
-  pybind11::array_t<arma::uword> getLabelsPython();
+  pybind11::array_t <arma::uword> getLabelsPython();
 
   /**
    * @brief Returns the number of swap steps
@@ -73,6 +73,14 @@ class KMedoidsWrapper : public km::KMedoids {
    */
   float getLossPython();
 
+  /**
+   * @brief Returns the loss after the BUILD step
+   *
+   * The loss after the BUILD step, i.e., the average distance from each point to its
+   * nearest medoid after the BUILD step
+   */
+  float getBuildLossPython();
+
   /**
    * @brief Returns the number of distance computations (sample complexity) used by .fit()
    *
@@ -127,74 +135,82 @@ class KMedoidsWrapper : public km::KMedoids {
 
 // TODO(@motiwari): Encapsulate these
 
-/**
- * @brief Binding for the C++ function KMedoids::fit
- */
-void fit_python(pybind11::class_<km::KMedoidsWrapper> *);
+  /**
+  * @brief Binding for the C++ function KMedoids::fit
+  */
+  void fit_python(pybind11::class_ <km::KMedoidsWrapper> *cls);
 
-/**
- * @brief Binding for the C++ function KMedoids::getMedoidsBuild()
- */
-void build_medoids_python(pybind11::class_<km::KMedoidsWrapper> *);
+  /**
+  * @brief Binding for the C++ function KMedoids::getMedoidsBuild()
+  */
+  void build_medoids_python(pybind11::class_ <km::KMedoidsWrapper> *cls);
 
-/**
- * @brief Binding for the C++ function KMedoids::getMedoidsFinal()
- */
-void medoids_python(pybind11::class_<km::KMedoidsWrapper> *);
+  /**
+  * @brief Binding for the C++ function KMedoids::getMedoidsFinal()
+  */
+  void medoids_python(pybind11::class_ <km::KMedoidsWrapper> *cls);
 
-/**
- * @brief Binding for the C++ function KMedoids::getLabels()
- */
-void labels_python(pybind11::class_<km::KMedoidsWrapper> *);
+  /**
+  * @brief Binding for the C++ function KMedoids::getLabels()
+  */
+  void labels_python(pybind11::class_ <km::KMedoidsWrapper> *cls);
 
-/**
- * @brief Binding for the C++ function KMedoids::getSteps()
- */
-void steps_python(pybind11::class_<km::KMedoidsWrapper> *);
+  /**
+  * @brief Binding for the C++ function KMedoids::getSteps()
+  */
+  void steps_python(pybind11::class_ <km::KMedoidsWrapper> *cls);
 
-/**
- * @brief Binding for the C++ function KMedoids::calcLoss()
- */
-void loss_python(pybind11::class_<km::KMedoidsWrapper> *);
+  /**
+  * @brief Binding for the C++ function KMedoids::calcLoss()
+  */
+  void loss_python(pybind11::class_ <km::KMedoidsWrapper> *cls);
 
-/**
- * @brief Binding for the C++ function KMedoids::getDistanceComputations()
- */
-void distance_computations_python(pybind11::class_<km::KMedoidsWrapper> *);
+  /**
+  * @brief Binding for the C++ function KMedoids::calcLoss() (after BUILD step)
+  */
+  void build_loss_python(pybind11::class_ <km::KMedoidsWrapper> *cls);
 
-// TODO(@motiwari): Add docstring
-void misc_distance_computations_python(pybind11::class_<km::KMedoidsWrapper> *);
+  /**
+  * @brief Binding for the C++ function KMedoids::getDistanceComputations()
+  */
+  void
+  distance_computations_python(pybind11::class_ <km::KMedoidsWrapper> *cls);
 
-// TODO(@motiwari): Add docstring
-void build_distance_computations_python(
-    pybind11::class_<km::KMedoidsWrapper> *);
+  // TODO(@motiwari): Add docstring
+  void misc_distance_computations_python(
+        pybind11::class_ <km::KMedoidsWrapper> *cls);
 
-// TODO(@motiwari): Add docstring
-void swap_distance_computations_python(pybind11::class_<km::KMedoidsWrapper> *);
+  // TODO(@motiwari): Add docstring
+  void build_distance_computations_python(
+        pybind11::class_ <km::KMedoidsWrapper> *cls);
 
-/**
- * @brief Binding for the C++ function KMedoids::getCacheWrites()
- */
-void cache_writes_python(pybind11::class_<km::KMedoidsWrapper> *);
+  // TODO(@motiwari): Add docstring
+  void swap_distance_computations_python(
+        pybind11::class_ <km::KMedoidsWrapper> *cls);
 
-/**
- * @brief Binding for the C++ function KMedoids::getCacheHits()
- */
-void cache_hits_python(pybind11::class_<km::KMedoidsWrapper> *);
+  /**
+  * @brief Binding for the C++ function KMedoids::getCacheWrites()
+  */
+  void cache_writes_python(pybind11::class_ <km::KMedoidsWrapper> *cls);
 
-/**
- * @brief Binding for the C++ function KMedoids::getCacheMisses()
- */
-void cache_misses_python(pybind11::class_<km::KMedoidsWrapper> *);
+  /**
+  * @brief Binding for the C++ function KMedoids::getCacheHits()
+  */
+  void cache_hits_python(pybind11::class_ <km::KMedoidsWrapper> *cls);
 
-/**
- * @brief Binding for the C++ function KMedoids::getTotalSwapTime()
- */
-void total_swap_time_python(pybind11::class_<km::KMedoidsWrapper> *);
+  /**
+  * @brief Binding for the C++ function KMedoids::getCacheMisses()
+  */
+  void cache_misses_python(pybind11::class_ <km::KMedoidsWrapper> *cls);
 
-/**
- * @brief Binding for the C++ function KMedoids::getTimePerSwap()
- */
-void time_per_swap_python(pybind11::class_<km::KMedoidsWrapper> *);
+  /**
+  * @brief Binding for the C++ function KMedoids::getTotalSwapTime()
+  */
+  void total_swap_time_python(pybind11::class_ <km::KMedoidsWrapper> *cls);
+
+  /**
+  * @brief Binding for the C++ function KMedoids::getTimePerSwap()
+  */
+  void time_per_swap_python(pybind11::class_ <km::KMedoidsWrapper> *cls);
 }  // namespace km
 #endif  // HEADERS_PYTHON_BINDINGS_KMEDOIDS_PYWRAPPER_HPP_
diff --git a/scripts/cache_measurements.py b/scripts/cache_measurements.py
index b3d29804..1fa6f9c9 100644
--- a/scripts/cache_measurements.py
+++ b/scripts/cache_measurements.py
@@ -56,13 +56,21 @@ def test_cache_stats():
             misses_1000,
         ) = get_cache_statistics(kmed=kmed, X=X, loss="L2", cache_width=1000)
         time_750, width_750, writes_750, hits_750, misses_750 = \
-            get_cache_statistics(kmed=kmed, X=X, loss="L2", cache_width=750)
+            get_cache_statistics(
+                kmed=kmed, X=X, loss="L2", cache_width=750
+            )
         time_500, width_500, writes_500, hits_500, misses_500 = \
-            get_cache_statistics(kmed=kmed, X=X, loss="L2", cache_width=500)
+            get_cache_statistics(
+                kmed=kmed, X=X, loss="L2", cache_width=500
+            )
         time_250, width_250, writes_250, hits_250, misses_250 = \
-            get_cache_statistics(kmed=kmed, X=X, loss="L2", cache_width=250)
+            get_cache_statistics(
+                kmed=kmed, X=X, loss="L2", cache_width=250
+            )
         time_0, width_0, writes_0, hits_0, misses_0 = \
-            get_cache_statistics(kmed=kmed, X=X, loss="L2", cache_width=0)
+            get_cache_statistics(
+                kmed=kmed, X=X, loss="L2", cache_width=0
+            )
 
         assert (
             hits_1000 > hits_750 > hits_500 > hits_250 > hits_0
@@ -70,10 +78,10 @@ def test_cache_stats():
         assert (
             misses_1000 < misses_750 < misses_500 < misses_250 < misses_0
         ), "Cache misses should decrease as cache size increases"
-        assert (
-            misses_1000 == 0
-        ), "There should be no cache misses at width=1000 as the" \
-           " whole dataset fits in memory"
+        assert misses_1000 == 0, (
+            "There should be no cache misses at width=1000 as the"
+            " whole dataset fits in memory"
+        )
         assert (
             hits_0 == 0
         ), "There should be no cache hits at width=0 as there is no cache"
@@ -84,8 +92,8 @@ def test_cache_stats():
         assert width_250 == 250, "Cache width should be 250 when set to 250"
         assert width_500 == 500, "Cache width should be 500 when set to 500"
         assert width_750 == 750, "Cache width should be 750 when set to 750"
-        assert width_1000 == 1000, "Cache width should be 1000 when set " \
-                                   "to 1000"
+        assert width_1000 == 1000, "Cache width should be 1000 when set to " \
+                                   "1000"
 
     def test_parallelization():
         X = np.loadtxt(os.path.join("data", "MNIST_10k.csv"))
diff --git a/scripts/comparison_with_fasterpam.py b/scripts/comparison_with_fasterpam.py
index b324e155..9d6043be 100644
--- a/scripts/comparison_with_fasterpam.py
+++ b/scripts/comparison_with_fasterpam.py
@@ -23,8 +23,11 @@ def benchmark(data, f, n=1):
         v = np.array(v)
         min, avg = v.min(), v.mean()
         ste = v.std(ddof=1) / np.sqrt(len(v)) if len(v) > 1 else 0.0
-        print("{:16s} min={:-10.2f} mean={:-10.2f} ±{:-.2f}"
-              .format(k, min, avg, ste))
+        print(
+            "{:16s} min={:-10.2f} mean={:-10.2f} ±{:-.2f}".format(
+                k, min, avg, ste
+            )
+        )
 
 
 def run_fasterpam(data, seed):
@@ -35,6 +38,7 @@ def run_fasterpam(data, seed):
     print("FasterPAM took ", r.n_iter, " iterations")
     meds, lbl = data[r.medoids], r.labels
     verified_loss = np.sqrt(((data - meds[lbl]) ** 2).sum(axis=1)).sum()
+
     return {
         "time (ms)": (end - start),
         "verified loss": verified_loss,
@@ -65,7 +69,8 @@ def run_old_bandit(data, seed):
         n_medoids=5,
         parallelize=True,
         algorithm="BanditPAM_orig",
-        dist_mat=diss)
+        dist_mat=diss
+    )
     print(km.algorithm)
     km.seed = seed
     start = time.time()
diff --git a/scripts/retrieve_windows_cmake_files.sh b/scripts/retrieve_windows_cmake_files.sh
new file mode 100644
index 00000000..c3afa432
--- /dev/null
+++ b/scripts/retrieve_windows_cmake_files.sh
@@ -0,0 +1,25 @@
+#!/bin/bash
+
+# This script is used to retrieve the files necessary for the Windows CMake build
+
+echo "Retrieving files for Windows CMake build..."
+
+# Move up one directory
+cd ..
+
+# Clone the GitHub repository
+git clone https://github.com/ThrunGroup/BanditPAM_Windows
+
+# Move unistd.h to headers/unistd.h
+mv BanditPAM_Windows/unistd.h headers/unistd.h
+#
+# Move getopt.h to headers/getopt.h
+mv BanditPAM_Windows/getopt.h headers/getopt.h
+#
+# Move getopt.cpp to src/getopt.cpp
+mv BanditPAM_Windows/getopt.cpp src/getopt.cpp
+
+# Delete BanditPAM_Windows directory
+rm -rf BanditPAM_Windows
+
+echo "Done!"
\ No newline at end of file
diff --git a/setup.py b/setup.py
index 58509bbc..9ac7d8e4 100644
--- a/setup.py
+++ b/setup.py
@@ -1,3 +1,4 @@
+import platform
 import sys
 import os
 import tempfile
@@ -8,7 +9,7 @@
 import distutils.sysconfig
 import distutils.spawn
 
-__version__ = "4.0.1"
+__version__ = "4.0.2"
 
 # TODO(@motiwari): Move this to a separate file
 GHA = "GITHUB_ACTIONS"
@@ -56,14 +57,19 @@ def compiler_check():
         )
     except KeyError:
         # The 'CC' environment variable hasn't been set.
-        # In this case, search for clang and gcc
+        # In this case, search for compilers that we can use.
         # Borrowed from https://github.com/clab/dynet/blob/master/setup.py
-        if distutils.spawn.find_executable("clang") is not None:
+        if distutils.spawn.find_executable("cl") is not None:
+            return "msvc"
+        elif distutils.spawn.find_executable("clang") is not None:
             return "clang"
         elif distutils.spawn.find_executable("gcc") is not None:
             return "gcc"
 
-    raise Exception("No C++ compiler was found. Please install LLVM clang.")
+    raise Exception(
+        "No C++ compiler was found. Please ensure you have "
+        "MSVC, LLVM clang, or GCC."
+    )
 
 
 def has_flag(compiler: str, flagname: str):
@@ -92,11 +98,11 @@ def cpp_flag(compiler: str):
     The newer version is prefered over c++11 (when it is available).
     """
     compiler_name = compiler_check()
-    if compiler_name == "clang":
-        flags = ["-std=c++17", "-std=c++14", "-std=c++11"]
-    else:
-        # Assume gcc
+    if compiler_name == "clang" or compiler_name == "gcc":
         flags = ["-std=c++17"]
+    else:
+        # Assume msvc
+        flags = ["/std:c++17"]  # required for std::optional
 
     for flag in flags:
         if has_flag(compiler, flag):
@@ -118,6 +124,7 @@ def check_brew_package(pkg_name: str):
             Please run `brew install %s`"
             % (pkg_name, pkg_name)
         )
+
     return output.decode().strip()
 
 
@@ -199,6 +206,7 @@ def check_armadillo_install_linux():
                 "armadillo-code",
             ),
         )
+
     return output.decode().strip()
 
 
@@ -216,6 +224,7 @@ def check_linux_package_installation(pkg_name: str):
             libgdbm-dev libc6-dev libbz2-dev libffi-dev zlib1g-dev"
             % (pkg_name)
         )
+
     return output.decode().strip()
 
 
@@ -263,8 +272,8 @@ def setup_colab(delete_source=False):
             + repo_location
         )
         os.system(
-            repo_location +
-            os.path.join(
+            repo_location
+            + os.path.join(
                 "scripts",
                 "colab_files",
                 "colab_install_armadillo.sh",
@@ -353,7 +362,11 @@ class BuildExt(build_ext):
         assert compiler_check() == "clang" or os.environ.get(
             GHA, False
         ), "Need to install LLVM clang!"
-        darwin_opts = ["-stdlib=libc++", "-mmacosx-version-min=10.14", "-O3"]
+        darwin_opts = [
+            "-stdlib=libc++",
+            "-mmacosx-version-min=10.14",
+            "-O3",
+        ]
         c_opts["unix"] += darwin_opts
         l_opts["unix"] += darwin_opts
     elif sys.platform == "linux" or sys.platform == "linux2":
@@ -363,6 +376,9 @@ class BuildExt(build_ext):
         linux_opts = ["-O3"]
         c_opts["unix"] += linux_opts
         l_opts["unix"] += linux_opts
+    # Currently necessary (unsure why)
+    elif sys.platform == "win32":
+        c_opts["msvc"] += ["/fsanitize=address"]
 
     def build_extensions(self):
         ct = self.compiler.compiler_type
@@ -374,21 +390,21 @@ def build_extensions(self):
         opts.append(cpp_flag(self.compiler))
         opts.append("-O3")
         if sys.platform == "darwin" and os.environ.get(GHA, False):
-            opts.append('-Xpreprocessor')  # NEEDS TO BE WITH NEXT LINE
-            opts.append('-fopenmp')  # NEEDS TO BE WITH PREVIOUS LINE
+            opts.append("-Xpreprocessor")  # NEEDS TO BE WITH NEXT LINE
+            opts.append("-fopenmp")  # NEEDS TO BE WITH PREVIOUS LINE
 
-            opts.append('-lomp')  # Potentially unused?
-            opts.append('-I/usr/local/opt/libomp/include')
-            opts.append('-L/usr/local/opt/libomp/lib')  # Potentially unused?
-        else:
+            opts.append("-lomp")  # Potentially unused?
+            opts.append("-I/usr/local/opt/libomp/include")
+            opts.append("-L/usr/local/opt/libomp/lib")  # Unused?
+        elif sys.platform != "win32":
             opts.append("-fopenmp")
 
         compiler_name = compiler_check()
-        if (sys.platform == "darwin" and os.environ.get(GHA, False)):
-            link_opts.append('-lomp')  # Potentially unused?
-            link_opts.append('-I/usr/local/opt/libomp/include')
-            link_opts.append('-L/usr/local/opt/libomp/lib')  # Unused?
-        else:
+        if sys.platform == "darwin" and os.environ.get(GHA, False):
+            link_opts.append("-lomp")  # Potentially unused?
+            link_opts.append("-I/usr/local/opt/libomp/include")
+            link_opts.append("-L/usr/local/opt/libomp/lib")  # Unused?
+        elif sys.platform != "win32":
             if compiler_name == "clang":
                 link_opts.append("-lomp")
             else:  # gcc
@@ -400,12 +416,16 @@ def build_extensions(self):
 
         for ext in self.extensions:
             ext.define_macros = [
-                ("VERSION_INFO", '"{}"'.format(
-                    self.distribution.get_version()
-                ))
+                ("VERSION_INFO",
+                 '"{}"'.format(self.distribution.get_version()))
             ]
             ext.extra_compile_args = opts
+            ext.extra_compile_args += []  # []["-arch", "x86_64"]
+
             ext.extra_link_args = link_opts
+            ext.extra_link_args += [
+                "-v",
+            ]  # "-arch", "x86_64"]
 
         build_ext.build_extensions(self)
 
@@ -422,8 +442,7 @@ def main():
             os.path.join("headers", "carma", "include", "carma_bits"),
             os.path.join("/", "usr", "local", "include"),
         ]
-
-    else:  # OSX
+    elif sys.platform == "darwin":  # OSX
         include_dirs = [
             get_pybind_include(),
             get_numpy_include(),
@@ -436,8 +455,9 @@ def main():
             # To include carma when the BanditPAM repo hasnt been initialized
             os.path.join("/", "usr", "local", "include"),
             os.path.join("/", "usr", "local", "include", "carma"),
-            os.path.join("/", "usr", "local", "include", "carma",
-                         "carma_bits"),
+            os.path.join(
+                "/", "usr", "local", "include", "carma", "carma_bits"
+            ),
             # When building from source on M1 Macs, may need these dirs
             # Currently, we should never be building from source on an M1 Mac,
             # Only cross-compiling from an Intel Mac
@@ -448,32 +468,79 @@ def main():
             os.path.join("/", "opt", "homebrew", "lib"),
             os.path.join("/", "opt", "homebrew", "opt"),
             os.path.join("/", "opt", "homebrew", "opt", "armadillo"),
-            os.path.join("/", "opt", "homebrew", "opt", "armadillo",
-                         "include"),
+            os.path.join(
+                "/", "opt", "homebrew", "opt", "armadillo", "include"
+            ),
             os.path.join(
                 "/", "opt", "homebrew", "opt", "armadillo", "include",
                 "armadillo_bits"
             ),
             # Needed for Mac Github Runners
             # for macos-10.15
-            os.path.join("/", "usr", "local", "Cellar", "libomp",
-                         "15.0.2", "include"),
+            os.path.join(
+                "/", "usr", "local", "Cellar", "libomp", "15.0.2", "include"
+            ),
             # for macos-latest
-            os.path.join("/", "usr", "local", "Cellar", "libomp",
-                         "15.0.7", "include"),
+            os.path.join(
+                "/", "usr", "local", "Cellar", "libomp", "15.0.7", "include"
+            ),
         ]
+    elif sys.platform == "win32":  # WIN32
+        include_dirs = [
+            get_pybind_include(),
+            get_numpy_include(),
+            "headers",
+            os.path.join("headers", "algorithms"),
+            os.path.join("headers", "python_bindings"),
+            os.path.join("headers", "carma", "include"),
+            os.path.join("headers", "carma", "include", "carma_bits"),
+            os.path.join("headers", "armadillo", "include"),
+            os.path.join("headers", "armadillo", "include", "armadillo_bits"),
+        ]
+    else:
+        raise Exception("Unrecognized platform")
 
     compiler_name = compiler_check()
-    if (sys.platform == "darwin" and os.environ.get(GHA, False)):
+    if sys.platform == "darwin" and os.environ.get(GHA, False):
         # On Mac Github Runners, we should NOT include gomp or omp here
         # due to build errors.
         libraries = ["armadillo", "omp"]
+    elif sys.platform == "win32":
+        libraries = ["libopenblas"]
     else:
         if compiler_name == "clang":
             libraries = ["armadillo", "omp"]
         else:  # gcc
             libraries = ["armadillo", "gomp"]
 
+    cpp_args = None
+    if sys.platform == "win32":
+        cpp_args = ["/std:c++17"]
+        library_dirs = [
+            # for windows
+            os.path.join(os.getcwd(), r"headers\armadillo\examples\lib_win64"),
+        ]
+    else:
+        cpp_args = [
+            "-static-libstdc++"
+        ]  # TODO(@motiwari): Modify this based on gcc or clang
+        library_dirs = [
+            os.path.join("/", "usr", "local", "lib"),
+            os.path.join(
+                "/", "usr", "local", "Cellar", "libomp", "15.0.2", "lib"
+            ),
+            os.path.join(
+                "/", "usr", "local", "Cellar", "libomp", "15.0.7", "lib"
+            ),
+        ]
+        if sys.platform == "darwin" and platform.processor() == "arm":  # M1
+            library_dirs.append(
+                os.path.join("/", "opt", "homebrew", "opt", "armadillo", "lib")
+            )
+            library_dirs.append(
+                os.path.join("/", "opt", "homebrew", "opt", "libomp", "lib")
+            )
+
     ext_modules = [
         Extension(
             "banditpam",
@@ -483,40 +550,49 @@ def main():
                 os.path.join("src", "algorithms", "banditpam.cpp"),
                 os.path.join("src", "algorithms", "banditpam_orig.cpp"),
                 os.path.join("src", "algorithms", "fastpam1.cpp"),
-                os.path.join("src", "python_bindings",
-                             "kmedoids_pywrapper.cpp"),
+                os.path.join(
+                    "src", "python_bindings", "kmedoids_pywrapper.cpp"
+                ),
                 os.path.join("src", "python_bindings", "medoids_python.cpp"),
-                os.path.join("src", "python_bindings",
-                             "build_medoids_python.cpp"),
+                os.path.join(
+                    "src", "python_bindings", "build_medoids_python.cpp"
+                ),
                 os.path.join("src", "python_bindings", "fit_python.cpp"),
                 os.path.join("src", "python_bindings", "labels_python.cpp"),
                 os.path.join("src", "python_bindings", "steps_python.cpp"),
                 os.path.join("src", "python_bindings", "loss_python.cpp"),
                 os.path.join("src", "python_bindings", "cache_python.cpp"),
-                os.path.join("src", "python_bindings",
-                             "swap_times_python.cpp"),
+                os.path.join(
+                    "src", "python_bindings", "swap_times_python.cpp"
+                ),
             ],
             include_dirs=include_dirs,
-            library_dirs=[
-                os.path.join("/", "usr", "local", "lib"),
-                # For Mac Github runners that install locally
-                # (not build wheels) - testing for macos-10.15
-                os.path.join("/", "usr", "local", "Cellar", "libomp",
-                             "15.0.2", "lib"),
-
-                # for macos-latest
-                os.path.join("/", "usr", "local", "Cellar", "libomp",
-                             "15.0.7", "lib"),
-            ],
+            library_dirs=library_dirs,
             libraries=libraries,
             language="c++1z",  # TODO: modify this based on cpp_flag(compiler)
-            extra_compile_args=["-static-libstdc++"],
+            extra_compile_args=cpp_args,
+            extra_link_args=[],  # Wrong pass (BuildExt sets extra_link_args)
         )
     ]
 
     with open(os.path.join("docs", "long_desc.rst"), encoding="utf-8") as f:
         long_description = f.read()
 
+    my_data_files = [("docs", [os.path.join("docs", "long_desc.rst")])]
+    if sys.platform == "win32":
+        my_data_files.append(
+            (
+                "",
+                [
+                    os.path.join(
+                        os.getcwd(),
+                        r"headers\armadillo\examples\lib_win64"
+                        + r"\libopenblas.dll",
+                        )
+                ],
+            )
+        )
+
     setup(
         name="banditpam",
         version=__version__,
@@ -529,7 +605,7 @@ def main():
         long_description=long_description,
         ext_modules=ext_modules,
         setup_requires=["pybind11>=2.5.0", "numpy>=1.18"],
-        data_files=[("docs", [os.path.join("docs", "long_desc.rst")])],
+        data_files=my_data_files,
         include_package_data=True,
         cmdclass={"build_ext": BuildExt},
         zip_safe=False,
@@ -543,8 +619,9 @@ def main():
             os.path.join("headers", "algorithms", "banditpam.hpp"),
             os.path.join("headers", "algorithms", "fastpam1.hpp"),
             os.path.join("headers", "algorithms", "pam.hpp"),
-            os.path.join("headers", "python_bindings",
-                         "kmedoids_pywrapper.hpp"),
+            os.path.join(
+                "headers", "python_bindings", "kmedoids_pywrapper.hpp"
+            ),
         ],
     )
 
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index ff19474e..7b3dc6a2 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -3,24 +3,55 @@ include_directories(${PROJECT_SOURCE_DIR}/headers/algorithms)
 include_directories(${PROJECT_SOURCE_DIR}/headers/python_bindings)
 include_directories(${PROJECT_SOURCE_DIR}/headers/carma/include)
 include_directories(${PROJECT_SOURCE_DIR}/headers/carma/include/carma_bits)
+# As per README instructions, we expected the user to have cloned armadillo locally
+if(WIN32)
+    include_directories(${PROJECT_SOURCE_DIR}/headers/armadillo/include)
+    include_directories(
+            ${PROJECT_SOURCE_DIR}/headers/armadillo/include/armadillo_bits)
+    set(OPENBLAS_dir ${PROJECT_SOURCE_DIR}/headers/armadillo/examples/lib_win64)
+endif()
 
-# For Mac Github Runner to find OpenMP -- potentially unnecessary after fixing -Xpreprocessor -fopenmp issue
-include_directories(/usr/local/Cellar/libomp/15.0.2/include)
-include_directories(/usr/local/Cellar/libomp/15.0.7/include)
-include_directories(/usr/local/Cellar/libomp)
-include_directories(/usr/local/opt/libomp/lib)
-include_directories(/usr/local/opt/libomp)
-include_directories(/usr/local/opt)
+# For Mac Github Runner to find OpenMP -- potentially unnecessary after fixing
+# -Xpreprocessor -fopenmp issue
+if(APPLE OR UNIX)
+    include_directories(/usr/local/Cellar/libomp/15.0.2/include)
+    include_directories(/usr/local/Cellar/libomp/15.0.7/include)
+    include_directories(/usr/local/Cellar/libomp)
+    include_directories(/usr/local/opt/libomp/include)
+    include_directories(/usr/local/opt/libomp/lib)
+    include_directories(/usr/local/opt/libomp)
+    include_directories(/usr/local/opt)
+endif()
 
-add_executable(BanditPAM main.cpp)
+if(WIN32)
+    add_executable(BanditPAM getopt.cpp main.cpp )
+elseif(APPLE OR UNIX)
+    add_executable(BanditPAM main.cpp )
+endif()
+add_library(
+        BanditPAM_LIB algorithms/kmedoids_algorithm.cpp
+        algorithms/pam.cpp
+        algorithms/banditpam.cpp
+        algorithms/banditpam_orig.cpp
+        algorithms/fastpam1.cpp)
 
-add_library(BanditPAM_LIB algorithms/kmedoids_algorithm.cpp algorithms/pam.cpp algorithms/banditpam.cpp algorithms/banditpam_orig.cpp algorithms/fastpam1.cpp)
 target_link_libraries(BanditPAM PUBLIC BanditPAM_LIB)
 
-find_package(OpenMP REQUIRED)
-target_link_libraries(BanditPAM_LIB PUBLIC OpenMP::OpenMP_CXX)
-target_link_libraries(BanditPAM PUBLIC OpenMP::OpenMP_CXX)
+if(WIN32)
+    target_link_directories(BanditPAM PUBLIC ${OPENBLAS_dir})
+    target_link_libraries(BanditPAM_LIB libopenblas)
+endif()
 
-find_package(Armadillo REQUIRED)
-include_directories(${ARMADILLO_INCLUDE_DIRS})
-target_link_libraries(BanditPAM PUBLIC ${ARMADILLO_LIBRARIES})
+if(APPLE OR UNIX)
+    find_package(Armadillo REQUIRED)
+    include_directories(${ARMADILLO_INCLUDE_DIRS})
+    target_link_libraries(BanditPAM PUBLIC ${ARMADILLO_LIBRARIES})
+endif()
+
+if(WIN32)
+    # Windows uses backslashes for paths while Unix systems use forward slashes
+    string(REPLACE "/" "\\"  OPENBLAS_dir_win ${OPENBLAS_dir})
+    # Copy all DLLs from OPENBLAS_dir_win to the output directory
+    add_custom_command(TARGET BanditPAM_LIB POST_BUILD
+                COMMAND copy ${OPENBLAS_dir_win}\\*.dll $(OutDir))
+endif()
diff --git a/src/algorithms/banditpam.cpp b/src/algorithms/banditpam.cpp
index 0e9fcd5f..aefbc271 100644
--- a/src/algorithms/banditpam.cpp
+++ b/src/algorithms/banditpam.cpp
@@ -13,543 +13,559 @@
 #include <vector>
 
 namespace km {
-void BanditPAM::fitBanditPAM(
-  const arma::fmat& inputData,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat) {
-  data = arma::trans(inputData);
-
-  // Note: even if we are using a distance matrix, we compute the permutation
-  // in the block below because it is used elsewhere in the call stack
-  // TODO(@motiwari): Remove need for data or permutation through when using
-  //  a distance matrix
-  if (this->useCache) {
-    size_t n = data.n_cols;
-    size_t m = fmin(n, cacheWidth);
-    cache = new float[n * m];
+  void BanditPAM::fitBanditPAM(
+          const arma::fmat &inputData,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat) {
+    data = arma::trans(inputData);
+
+    // Note: even if we are using a distance matrix, we compute the permutation
+    // in the block below because it is used elsewhere in the call stack
+    // TODO(@motiwari): Remove need for data or permutation through when using
+    //  a distance matrix
+    if (this->useCache) {
+      size_t n = data.n_cols;
+      size_t m = fmin(n, cacheWidth);
+      cache = new float[n * m];
+
+      #pragma omp parallel for if (this->parallelize)
+      for (size_t idx = 0; idx < m * n; idx++) {
+        cache[idx] = -1;  // TODO(@motiwari): need better value here
+      }
 
-    #pragma omp parallel for if (this->parallelize)
-    for (size_t idx = 0; idx < m*n; idx++) {
-      cache[idx] = -1;  // TODO(@motiwari): need better value here
+      permutation = arma::randperm(n);
+      permutationIdx = 0;
+      reindex = {};  // TODO(@motiwari): Can this intialization be removed?
+      // TODO(@motiwari): Can we parallelize this?
+      for (size_t counter = 0; counter < m; counter++) {
+        reindex[permutation[counter]] = counter;
+      }
     }
 
-    permutation = arma::randperm(n);
-    permutationIdx = 0;
-    reindex = {};  // TODO(@motiwari): Can this intialization be removed?
-    // TODO(@motiwari): Can we parallelize this?
-    for (size_t counter = 0; counter < m; counter++) {
-      reindex[permutation[counter]] = counter;
+    arma::fmat medoidMatrix(data.n_rows, nMedoids);
+    arma::urowvec medoidIndices(nMedoids);
+    steps = 0;
+    BanditPAM::build(data, distMat, &medoidIndices, &medoidMatrix);
+
+    buildLoss = KMedoids::calcLoss(data, distMat, &medoidIndices);
+
+    medoidIndicesBuild = medoidIndices;
+    arma::urowvec assignments(data.n_cols);
+    if (nMedoids > 1) {
+      BanditPAM::swap(
+              data,
+              distMat,
+              &medoidIndices,
+              &medoidMatrix,
+              &assignments);
     }
+
+    medoidIndicesFinal = medoidIndices;
+    labels = assignments;
   }
 
-  arma::fmat medoidMatrix(data.n_rows, nMedoids);
-  arma::urowvec medoidIndices(nMedoids);
-  steps = 0;
-  BanditPAM::build(data, distMat, &medoidIndices, &medoidMatrix);
-
-  medoidIndicesBuild = medoidIndices;
-  arma::urowvec assignments(data.n_cols);
-
-  BanditPAM::swap(data, distMat, &medoidIndices, &medoidMatrix, &assignments);
-
-  medoidIndicesFinal = medoidIndices;
-  labels = assignments;
-}
-
-arma::frowvec BanditPAM::buildSigma(
-  const arma::fmat& data,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-  const arma::frowvec& bestDistances,
-  const bool useAbsolute) {
-  size_t N = data.n_cols;
-  arma::uvec referencePoints;
-  // TODO(@motiwari): Make this wraparound properly as
-  //  last batch_size elements are dropped
-  if (usePerm) {
-    if ((permutationIdx + batchSize - 1) >= N) {
-      permutationIdx = 0;
+  arma::frowvec BanditPAM::buildSigma(
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const arma::frowvec &bestDistances,
+          const bool useAbsolute) {
+    size_t N = data.n_cols;
+    arma::uvec referencePoints;
+    // TODO(@motiwari): Make this wraparound properly as
+    //  last batch_size elements are dropped
+    if (usePerm) {
+      if ((permutationIdx + batchSize - 1) >= N) {
+        permutationIdx = 0;
+      }
+      // inclusive of both indices
+      referencePoints = permutation.subvec(
+              permutationIdx,
+              permutationIdx + batchSize - 1);
+      permutationIdx += batchSize;
+    } else {
+      referencePoints = arma::randperm(N, batchSize);
     }
-    // inclusive of both indices
-    referencePoints = permutation.subvec(
-      permutationIdx,
-      permutationIdx + batchSize - 1);
-    permutationIdx += batchSize;
-  } else {
-    referencePoints = arma::randperm(N, batchSize);
-  }
 
-  arma::fvec sample(batchSize);
-  arma::frowvec updated_sigma(N);
-  #pragma omp parallel for if (this->parallelize)
-  for (size_t i = 0; i < N; i++) {
-    for (size_t j = 0; j < batchSize; j++) {
-      // 0 for MISC
-      float cost =
-          KMedoids::cachedLoss(data, distMat, i, referencePoints(j), 0);
-      if (useAbsolute) {
-        sample(j) = cost;
-      } else {
-        sample(j) = cost < bestDistances(referencePoints(j))
-                          ? cost : bestDistances(referencePoints(j));
-        sample(j) -= bestDistances(referencePoints(j));
+    arma::fvec sample(batchSize);
+    arma::frowvec updated_sigma(N);
+    #pragma omp parallel for if (this->parallelize)
+    for (size_t i = 0; i < N; i++) {
+      for (size_t j = 0; j < batchSize; j++) {
+        // 0 for MISC
+        float cost =
+                KMedoids::cachedLoss(data, distMat, i,
+                                     referencePoints(j), 0);
+        if (useAbsolute) {
+          sample(j) = cost;
+        } else {
+          sample(j) = cost < bestDistances(referencePoints(j))
+                        ? cost : bestDistances(referencePoints(j));
+          sample(j) -= bestDistances(referencePoints(j));
+        }
       }
+      updated_sigma(i) = arma::stddev(sample);
     }
-    updated_sigma(i) = arma::stddev(sample);
-  }
-  return updated_sigma;
-}
-
-arma::frowvec BanditPAM::buildTarget(
-  const arma::fmat& data,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-  const arma::uvec* target,
-  const arma::frowvec* bestDistances,
-  const bool useAbsolute,
-  const size_t exact = 0) {
-  size_t N = data.n_cols;
-  size_t tmpBatchSize = batchSize;
-  if (exact > 0) {
-    tmpBatchSize = N;
+    return updated_sigma;
   }
-  arma::frowvec results(target->n_rows, arma::fill::zeros);
-  arma::uvec referencePoints;
-  // TODO(@motiwari): Make this wraparound properly
-  //  as last batch_size elements are dropped
-  if (usePerm) {
-    if ((permutationIdx + tmpBatchSize - 1) >= N) {
-      permutationIdx = 0;
+
+  arma::frowvec BanditPAM::buildTarget(
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const arma::uvec *target,
+          const arma::frowvec *bestDistances,
+          const bool useAbsolute,
+          const bool exact = false) {
+    size_t N = data.n_cols;
+    size_t tmpBatchSize = batchSize;
+    if (exact) {
+      tmpBatchSize = N;
+    }
+    arma::frowvec results(target->n_rows, arma::fill::zeros);
+    arma::uvec referencePoints;
+    // TODO(@motiwari): Make this wraparound properly
+    //  as last batch_size elements are dropped
+    if (usePerm) {
+      if ((permutationIdx + tmpBatchSize - 1) >= N) {
+        permutationIdx = 0;
+      }
+      // inclusive of both indices
+      referencePoints = permutation.subvec(
+              permutationIdx,
+              permutationIdx + tmpBatchSize - 1);
+      permutationIdx += tmpBatchSize;
+    } else {
+      referencePoints = arma::randperm(N, tmpBatchSize);
     }
-    // inclusive of both indices
-    referencePoints = permutation.subvec(
-      permutationIdx,
-      permutationIdx + tmpBatchSize - 1);
-    permutationIdx += tmpBatchSize;
-  } else {
-    referencePoints = arma::randperm(N, tmpBatchSize);
-  }
 
-  #pragma omp parallel for if (this->parallelize)
-  for (size_t i = 0; i < target->n_rows; i++) {
-    float total = 0;
-    for (size_t j = 0; j < referencePoints.n_rows; j++) {
-      float cost =
-        KMedoids::cachedLoss(
-            data,
-            distMat,
-            (*target)(i),
-            referencePoints(j),
-            1);  // 1 for BUILD
-      if (useAbsolute) {
-        total += cost;
-      } else {
-        total += cost < (*bestDistances)(referencePoints(j))
-                      ? cost : (*bestDistances)(referencePoints(j));
-        total -= (*bestDistances)(referencePoints(j));
+    #pragma omp parallel for if (this->parallelize)
+    for (size_t i = 0; i < target->n_rows; i++) {
+        float total = 0;
+      for (size_t j = 0; j < referencePoints.n_rows; j++) {
+          float cost =
+            KMedoids::cachedLoss(
+                    data,
+                    distMat,
+                    (*target)(i),
+                    referencePoints(j),
+                    1);  // 1 for BUILD
+          if (useAbsolute) {
+            total += cost;
+          } else {
+            total += cost < (*bestDistances)(referencePoints(j))
+                       ? cost : (*bestDistances)(referencePoints(j));
+            total -= (*bestDistances)(referencePoints(j));
+          }
       }
+      results(i) = total / tmpBatchSize;
     }
-    results(i) = total / tmpBatchSize;
+    return results;
   }
-  return results;
-}
-
-void BanditPAM::build(
-  const arma::fmat& data,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-  arma::urowvec* medoidIndices,
-  arma::fmat* medoids) {
-  size_t N = data.n_cols;
-  arma::frowvec N_mat(N);
-  N_mat.fill(N);
-  size_t p = N;
-  bool useAbsolute = true;
-  arma::frowvec estimates(N, arma::fill::zeros);
-  arma::frowvec bestDistances(N);
-  bestDistances.fill(std::numeric_limits<float>::infinity());
-  arma::frowvec sigma(N);
-  arma::urowvec candidates(N, arma::fill::ones);
-  arma::frowvec lcbs(N);
-  arma::frowvec ucbs(N);
-  arma::frowvec numSamples(N, arma::fill::zeros);
-  arma::frowvec exactMask(N, arma::fill::zeros);
-
-  // TODO(@motiwari): #pragma omp parallel for if (this->parallelize)?
-  for (size_t k = 0; k < nMedoids; k++) {
-    // instantiate medoids one-by-one
-    permutationIdx = 0;
-    candidates.fill(1);
-    numSamples.fill(0);
-    exactMask.fill(0);
-    estimates.fill(0);
-    // compute std dev amongst batch of reference points
-    sigma = buildSigma(data, distMat, bestDistances, useAbsolute);
-
-    while (arma::sum(candidates) > precision) {
-      // TODO(@motiwari): Do not need a matrix for this comparison,
-      //  use broadcasting
-      arma::umat compute_exactly =
-        ((numSamples + batchSize) >= N_mat) != exactMask;
-      if (arma::accu(compute_exactly) > 0) {
-        arma::uvec targets = find(compute_exactly);
+
+  void BanditPAM::build(
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          arma::urowvec *medoidIndices,
+          arma::fmat *medoids) {
+    size_t N = data.n_cols;
+    arma::frowvec N_mat(N);
+    N_mat.fill(N);
+    size_t p = N;
+    bool useAbsolute = true;
+    arma::frowvec estimates(N, arma::fill::zeros);
+    arma::frowvec bestDistances(N);
+    bestDistances.fill(std::numeric_limits<float>::infinity());
+    arma::frowvec sigma(N);
+    arma::urowvec candidates(N, arma::fill::ones);
+    arma::frowvec lcbs(N);
+    arma::frowvec ucbs(N);
+    arma::frowvec numSamples(N, arma::fill::zeros);
+    arma::frowvec exactMask(N, arma::fill::zeros);
+
+    // TODO(@motiwari): #pragma omp parallel for if (this->parallelize)?
+    for (size_t k = 0; k < nMedoids; k++) {
+      // instantiate medoids one-by-one
+      permutationIdx = 0;
+      candidates.fill(1);
+      numSamples.fill(0);
+      exactMask.fill(0);
+      estimates.fill(0);
+      // compute std dev amongst batch of reference points
+      sigma = buildSigma(data, distMat, bestDistances, useAbsolute);
+
+      while (arma::sum(candidates) > precision) {
+        // TODO(@motiwari): Do not need a matrix for this comparison,
+        //  use broadcasting
+        arma::umat compute_exactly =
+                ((numSamples + batchSize) >= N_mat) != exactMask;
+        if (arma::accu(compute_exactly) > 0) {
+          arma::uvec targets = find(compute_exactly);
+          arma::frowvec result = buildTarget(
+                  data,
+                  distMat,
+                  &targets,
+                  &bestDistances,
+                  useAbsolute,
+                  (true ? N > 0 : false));
+          estimates.cols(targets) = result;
+          ucbs.cols(targets) = result;
+          lcbs.cols(targets) = result;
+          exactMask.cols(targets).fill(1);
+          numSamples.cols(targets) += N;
+          candidates.cols(targets).fill(0);
+        }
+        if (arma::sum(candidates) < precision) {
+          break;
+        }
+        arma::uvec targets = arma::find(candidates);
         arma::frowvec result = buildTarget(
-          data,
-          distMat,
-          &targets,
-          &bestDistances,
-          useAbsolute,
-          N);
-        estimates.cols(targets) = result;
-        ucbs.cols(targets) = result;
-        lcbs.cols(targets) = result;
-        exactMask.cols(targets).fill(1);
-        numSamples.cols(targets) += N;
-        candidates.cols(targets).fill(0);
-      }
-      if (arma::sum(candidates) < precision) {
-        break;
+                data,
+                distMat,
+                &targets,
+                &bestDistances,
+                useAbsolute,
+                false);
+        // update the running average
+        estimates.cols(targets) =
+                ((numSamples.cols(targets) % estimates.cols(targets)) +
+                 (result * batchSize)) /
+                (batchSize + numSamples.cols(targets));
+        numSamples.cols(targets) += batchSize;
+        arma::frowvec adjust(targets.n_rows);
+        adjust.fill(p);
+        // Assume buildConfidence is given in logspace
+        adjust = buildConfidence + arma::log(adjust);
+        arma::frowvec confBoundDelta =
+                sigma.cols(targets) %
+                arma::sqrt(adjust / numSamples.cols(targets));
+        ucbs.cols(targets) = estimates.cols(targets) + confBoundDelta;
+        lcbs.cols(targets) = estimates.cols(targets) - confBoundDelta;
+        candidates = (lcbs < ucbs.min()) && (exactMask == 0);
       }
-      arma::uvec targets = arma::find(candidates);
-      arma::frowvec result = buildTarget(
-        data,
-        distMat,
-        &targets,
-        &bestDistances,
-        useAbsolute,
-        0);
-      // update the running average
-      estimates.cols(targets) =
-        ((numSamples.cols(targets) % estimates.cols(targets)) +
-        (result * batchSize)) /
-        (batchSize + numSamples.cols(targets));
-      numSamples.cols(targets) += batchSize;
-      arma::frowvec adjust(targets.n_rows);
-      adjust.fill(p);
-      // Assume buildConfidence is given in logspace
-      adjust = buildConfidence + arma::log(adjust);
-      arma::frowvec confBoundDelta =
-        sigma.cols(targets) %
-        arma::sqrt(adjust / numSamples.cols(targets));
-      ucbs.cols(targets) = estimates.cols(targets) + confBoundDelta;
-      lcbs.cols(targets) = estimates.cols(targets) - confBoundDelta;
-      candidates = (lcbs < ucbs.min()) && (exactMask == 0);
-    }
 
-    medoidIndices->at(k) = lcbs.index_min();
-    medoids->unsafe_col(k) = data.unsafe_col((*medoidIndices)(k));
+      medoidIndices->at(k) = lcbs.index_min();
+      medoids->unsafe_col(k) = data.unsafe_col((*medoidIndices)(k));
 
-    // don't need to do this on final iteration
-    #pragma omp parallel for if (this->parallelize)
-    for (size_t i = 0; i < N; i++) {
+      // don't need to do this on final iteration
+      #pragma omp parallel for if (this->parallelize)
+      for (size_t i = 0; i < N; i++) {
         float cost = KMedoids::cachedLoss(
-          data,
-          distMat,
-          i,
-          (*medoidIndices)(k),
-          0);  // 0 for MISC
+                data,
+                distMat,
+                i,
+                (*medoidIndices)(k),
+                0);  // 0 for MISC
         if (cost < bestDistances(i)) {
-            bestDistances(i) = cost;
+          bestDistances(i) = cost;
         }
+      }
+      // use difference of loss for sigma and sampling, not absolute
+      useAbsolute = false;
     }
-    // use difference of loss for sigma and sampling, not absolute
-    useAbsolute = false;
   }
-}
-
-arma::fmat BanditPAM::swapSigma(
-  const arma::fmat& data,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-  const arma::frowvec* bestDistances,
-  const arma::frowvec* secondBestDistances,
-  const arma::urowvec* assignments) {
-  size_t N = data.n_cols;
-  size_t K = nMedoids;
-  arma::fmat updated_sigma(K, N, arma::fill::zeros);
-  arma::uvec referencePoints;
-  // TODO(@motiwari): Make this wraparound properly
-  //  as last batch_size elements are dropped
-  if (usePerm) {
-    if ((permutationIdx + batchSize - 1) >= N) {
-      permutationIdx = 0;
+
+  arma::fmat BanditPAM::swapSigma(
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const arma::frowvec *bestDistances,
+          const arma::frowvec *secondBestDistances,
+          const arma::urowvec *assignments) {
+    size_t N = data.n_cols;
+    size_t K = nMedoids;
+    arma::fmat updated_sigma(K, N, arma::fill::zeros);
+    arma::uvec referencePoints;
+    // TODO(@motiwari): Make this wraparound properly
+    //  as last batch_size elements are dropped
+    if (usePerm) {
+      if ((permutationIdx + batchSize - 1) >= N) {
+        permutationIdx = 0;
+      }
+      // inclusive of both indices
+      referencePoints = permutation.subvec(
+              permutationIdx,
+              permutationIdx + batchSize - 1);
+      permutationIdx += batchSize;
+    } else {
+      referencePoints = arma::randperm(N, batchSize);
     }
-    // inclusive of both indices
-    referencePoints = permutation.subvec(
-      permutationIdx,
-      permutationIdx + batchSize - 1);
-    permutationIdx += batchSize;
-  } else {
-    referencePoints = arma::randperm(N, batchSize);
-  }
 
-  arma::fvec sample(batchSize);
-  // for each considered swap
-  #pragma omp parallel for if (this->parallelize)
-  for (size_t i = 0; i < K * N; i++) {
-    // extract data point of swap
-    size_t n = i / K;
-    size_t k = i % K;
-
-    // calculate change in loss for some subset of the data
-    for (size_t j = 0; j < batchSize; j++) {
-      // 0 for MISC when estimating sigma
-      float cost =
-          KMedoids::cachedLoss(data, distMat, n, referencePoints(j), 0);
-
-      if (k == (*assignments)(referencePoints(j))) {
-        if (cost < (*secondBestDistances)(referencePoints(j))) {
-          sample(j) = cost;
-        } else {
-          sample(j) = (*secondBestDistances)(referencePoints(j));
-        }
-      } else {
-        if (cost < (*bestDistances)(referencePoints(j))) {
-          sample(j) = cost;
+    arma::fvec sample(batchSize);
+    // for each considered swap
+    #pragma omp parallel for if (this->parallelize)
+    for (size_t i = 0; i < K * N; i++) {
+      // extract data point of swap
+      size_t n = i / K;
+      size_t k = i % K;
+
+      // calculate change in loss for some subset of the data
+      for (size_t j = 0; j < batchSize; j++) {
+        // 0 for MISC when estimating sigma
+        float cost =
+          KMedoids::cachedLoss(data, distMat, n,
+                               referencePoints(j), 0);
+
+        if (k == (*assignments)(referencePoints(j))) {
+          if (cost < (*secondBestDistances)(referencePoints(j))) {
+            sample(j) = cost;
+          } else {
+            sample(j) = (*secondBestDistances)(referencePoints(j));
+          }
         } else {
-          sample(j) = (*bestDistances)(referencePoints(j));
+          if (cost < (*bestDistances)(referencePoints(j))) {
+            sample(j) = cost;
+          } else {
+            sample(j) = (*bestDistances)(referencePoints(j));
+          }
         }
+        sample(j) -= (*bestDistances)(referencePoints(j));
       }
-      sample(j) -= (*bestDistances)(referencePoints(j));
+      updated_sigma(k, n) = arma::stddev(sample);
     }
-    updated_sigma(k, n) = arma::stddev(sample);
-  }
-  return updated_sigma;
-}
-
-arma::fmat BanditPAM::swapTarget(
-  const arma::fmat& data,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-  const arma::urowvec* medoidIndices,
-  const arma::uvec* targets,
-  const arma::frowvec* bestDistances,
-  const arma::frowvec* secondBestDistances,
-  const arma::urowvec* assignments,
-  const size_t exact = 0) {
-  const size_t N = data.n_cols;
-  const size_t T = targets->n_rows;
-  arma::fmat results(nMedoids, T, arma::fill::zeros);
-
-  // Targets should be a list of indices for target CANDIDATE points
-  // Then update all corresponding EXISTING MEDOID indices targets.
-  // If targets is a T-length vector, then the return value should be
-  // a matrix of size K x T. We should perform the appropriate update then
-  // in the swap() function.
-  //
-  // An alternate method to do this would be to pass only the (m, c)
-  // Points under consideration. Then we wouldn't need to update all
-  // k virtual arms for each candidate, just the ones that are passed
-  // However, this would incur a .find() call to find all pairs
-  // (m, c) where c == c', the arm under consideration. I believe this
-  // would be an O(kn) cost. Instead, may need to use another data
-  // structure to avoid this .find() call, like a tree where the top-level
-  // nodes are the candidates and the bottom-level nodes are the corresponding
-  // virtual arms.
-  // A jagged array might also do the trick.
-
-
-  size_t tmpBatchSize = batchSize;
-  if (exact > 0) {
-    tmpBatchSize = N;
+    return updated_sigma;
   }
 
-  arma::uvec referencePoints;
-  // TODO(@motiwari): Make this wraparound properly
-  //  as last batch_size elements are dropped
-  // TODO(@motiwari): Break this duplicated code into a function
-  if (usePerm) {
-    if ((permutationIdx + tmpBatchSize - 1) >= N) {
-      permutationIdx = 0;
+  arma::fmat BanditPAM::swapTarget(
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const arma::urowvec *medoidIndices,
+          const arma::uvec *targets,
+          const arma::frowvec *bestDistances,
+          const arma::frowvec *secondBestDistances,
+          const arma::urowvec *assignments,
+          const bool exact = false) {
+    const size_t N = data.n_cols;
+    const size_t T = targets->n_rows;
+    arma::fmat results(nMedoids, T, arma::fill::zeros);
+
+    // Targets should be a list of indices for target CANDIDATE points
+    // Then update all corresponding EXISTING MEDOID indices targets.
+    // If targets is a T-length vector, then the return value should be
+    // a matrix of size K x T. We should perform the appropriate update then
+    // in the swap() function.
+    //
+    // An alternate method to do this would be to pass only the (m, c)
+    // Points under consideration. Then we wouldn't need to update all
+    // k virtual arms for each candidate, just the ones that are passed
+    // However, this would incur a .find() call to find all pairs
+    // (m, c) where c == c', the arm under consideration. I believe this
+    // would be an O(kn) cost. Instead, may need to use another data
+    // structure to avoid this .find() call, like a tree where the top-level
+    // nodes are the candidates and the bottom-level nodes are the corresponding
+    // virtual arms.
+    // A jagged array might also do the trick.
+
+
+    size_t tmpBatchSize = batchSize;
+    if (exact) {
+      tmpBatchSize = N;
     }
-    // inclusive of both indices
-    referencePoints = permutation.subvec(
-      permutationIdx,
-      permutationIdx + tmpBatchSize - 1);
-    permutationIdx += tmpBatchSize;
-  } else {
-    referencePoints = arma::randperm(N, tmpBatchSize);
-  }
 
-  // TODO(@motiwari): Declare variables outside of loops
-  #pragma omp parallel for if (this->parallelize)
-  for (size_t i = 0; i < T; i++) {
-    // TODO(@motiwari): pragma omp parallel for?
-    for (size_t j = 0; j < tmpBatchSize; j++) {
-      float cost =
-          KMedoids::cachedLoss(
-              data,
-              distMat,
-              (*targets)(i),
-              referencePoints(j),
-              2);  // 2 for SWAP
-      size_t k = (*assignments)(referencePoints(j));
-      if (cost < (*bestDistances)(referencePoints(j))) {
-          // We might be able to change this to .eachrow(every column but k)
+    arma::uvec referencePoints;
+    // TODO(@motiwari): Make this wraparound properly
+    //  as last batch_size elements are dropped
+    // TODO(@motiwari): Break this duplicated code into a function
+    if (usePerm) {
+      if ((permutationIdx + tmpBatchSize - 1) >= N) {
+        permutationIdx = 0;
+      }
+      // inclusive of both indices
+      referencePoints = permutation.subvec(
+              permutationIdx,
+              permutationIdx + tmpBatchSize - 1);
+      permutationIdx += tmpBatchSize;
+    } else {
+      referencePoints = arma::randperm(N, tmpBatchSize);
+    }
+
+    // TODO(@motiwari): Declare variables outside of loops
+    #pragma omp parallel for if (this->parallelize)
+    for (size_t i = 0; i < T; i++) {
+      // TODO(@motiwari): pragma omp parallel for?
+      for (size_t j = 0; j < tmpBatchSize; j++) {
+        float cost =
+                KMedoids::cachedLoss(
+                        data,
+                        distMat,
+                        (*targets)(i),
+                        referencePoints(j),
+                        2);  // 2 for SWAP
+        size_t k = (*assignments)(referencePoints(j));
+        if (cost < (*bestDistances)(referencePoints(j))) {
+          // We might be able to change this to
+          // .eachrow(every column but k)
           // since arma does this in-place and it should not introduce
           // complexity
-          results.col(i) += cost - (*bestDistances)(referencePoints(j));
-      }
+          results.col(i) +=
+                  cost - (*bestDistances)(referencePoints(j));
+        }
 
-      // If cost < bd, this second term will subtract off the "new cost"
-      // added by the all-column call above inside the if
-      results(k, i) +=
-        std::fmin(cost, (*secondBestDistances)(referencePoints(j))) -
-          std::fmin(cost, (*bestDistances)(referencePoints(j)));
+        // If cost < bd, this second term will subtract off the "new cost"
+        // added by the all-column call above inside the if
+        results(k, i) +=
+                std::fmin(cost,
+                          (*secondBestDistances)(referencePoints(j))) -
+                std::fmin(cost, (*bestDistances)(referencePoints(j)));
+      }
     }
+    // TODO(@motiwari): we can probably avoid this division
+    //  if we look at total loss, not average loss
+    results /= tmpBatchSize;
+    return results;
   }
-  // TODO(@motiwari): we can probably avoid this division
-  //  if we look at total loss, not average loss
-  results /= tmpBatchSize;
-  return results;
-}
-
-void BanditPAM::swap(
-  const arma::fmat& data,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-  arma::urowvec* medoidIndices,
-  arma::fmat* medoids,
-  arma::urowvec* assignments) {
-  size_t N = data.n_cols;
-  size_t p = N;
-
-  arma::fmat sigma(nMedoids, N, arma::fill::zeros);
-
-  arma::frowvec bestDistances(N);
-  arma::frowvec secondBestDistances(N);
-  bool swapPerformed = true;
-  arma::umat candidates(nMedoids, N, arma::fill::ones);
-  arma::umat exactMask(nMedoids, N, arma::fill::zeros);
-  arma::fmat estimates(nMedoids, N, arma::fill::zeros);
-  arma::fmat lcbs(nMedoids, N);
-  arma::fmat ucbs(nMedoids, N);
-  arma::umat numSamples(nMedoids, N, arma::fill::zeros);
-
-  // calculate quantities needed for swap, bestDistances and sigma
-  calcBestDistancesSwap(
-    data,
-    distMat,
-    medoidIndices,
-    &bestDistances,
-    &secondBestDistances,
-    assignments,
-    swapPerformed);
-
-  // continue making swaps while loss is decreasing
-  while (swapPerformed && steps < maxIter) {
-    steps++;
-    permutationIdx = 0;
-
-    sigma = swapSigma(
-      data,
-      distMat,
-      &bestDistances,
-      &secondBestDistances,
-      assignments);
-
-    // Reset variables when starting a new swap
-    candidates.fill(1);
-    exactMask.fill(0);
-    estimates.fill(0);
-    numSamples.fill(0);
-
-    // while there is at least one candidate (float comparison issues)
-    while (arma::accu(candidates) > 1.5) {
-      // compute exactly if it's been samples more than N times and
-      // hasn't been computed exactly already
-      arma::umat compute_exactly =
-        ((numSamples + batchSize) >= N) != (exactMask);
-
-      // Get unique candidate medoids from the candidates (second index)
-      // Store all k x T in estimates
-      // TODO(@motiwari): Move this declaration outside loop
-      // Need unique values over second index
-      // Sum the different columns
-      // if any index appears in at least one, compute it exactly
-      // TODO(@motiwari): make sure we're only computing exactly
-      // for the relevant candidates
-      arma::uvec compute_exactly_targets =
-        arma::find(arma::sum(compute_exactly, 0) >= 1);
-
-      if (compute_exactly_targets.size() > 0) {
-        arma::fmat result = swapTarget(
-          data,
-          distMat,
-          medoidIndices,
-          &compute_exactly_targets,
-          &bestDistances,
-          &secondBestDistances,
-          assignments,
-          N);
-
-        // result will be k x T
-        // Now update the correct indices
-        estimates.cols(compute_exactly_targets) = result;
-        ucbs.cols(compute_exactly_targets) = result;
-        lcbs.cols(compute_exactly_targets) = result;
-        exactMask.cols(compute_exactly_targets).fill(1);
-        numSamples.cols(compute_exactly_targets) += N;
-        candidates = (lcbs < ucbs.min()) && (exactMask == 0);
-      }
-      if (arma::accu(candidates) < precision) {
-        break;
+
+  void BanditPAM::swap(
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          arma::urowvec *medoidIndices,
+          arma::fmat *medoids,
+          arma::urowvec *assignments) {
+    size_t N = data.n_cols;
+    size_t p = N;
+
+    arma::fmat sigma(nMedoids, N, arma::fill::zeros);
+
+    arma::frowvec bestDistances(N);
+    arma::frowvec secondBestDistances(N);
+    bool swapPerformed = true;
+    arma::umat candidates(nMedoids, N, arma::fill::ones);
+    arma::umat exactMask(nMedoids, N, arma::fill::zeros);
+    arma::fmat estimates(nMedoids, N, arma::fill::zeros);
+    arma::fmat lcbs(nMedoids, N);
+    arma::fmat ucbs(nMedoids, N);
+    arma::umat numSamples(nMedoids, N, arma::fill::zeros);
+
+    // calculate quantities needed for swap, bestDistances and sigma
+    calcBestDistancesSwap(
+            data,
+            distMat,
+            medoidIndices,
+            &bestDistances,
+            &secondBestDistances,
+            assignments,
+            swapPerformed);
+
+    // continue making swaps while loss is decreasing
+    while (swapPerformed && steps < maxIter) {
+        steps++;
+        permutationIdx = 0;
+
+        sigma = swapSigma(
+                data,
+                distMat,
+                &bestDistances,
+                &secondBestDistances,
+                assignments);
+
+        // Reset variables when starting a new swap
+        candidates.fill(1);
+        exactMask.fill(0);
+        estimates.fill(0);
+        numSamples.fill(0);
+
+        // while there is at least one candidate (float comparison issues)
+        while (arma::accu(candidates) > 1.5) {
+          // compute exactly if it's been samples more than N times and
+          // hasn't been computed exactly already
+          arma::umat compute_exactly =
+                  ((numSamples + batchSize) >= N) != (exactMask);
+
+          // Get unique candidate medoids from the candidates (second index)
+          // Store all k x T in estimates
+          // TODO(@motiwari): Move this declaration outside loop
+          // Need unique values over second index
+          // Sum the different columns
+          // if any index appears in at least one, compute it exactly
+          // TODO(@motiwari): make sure we're only computing exactly
+          // for the relevant candidates
+          arma::uvec compute_exactly_targets =
+                  arma::find(arma::sum(compute_exactly, 0) >= 1);
+
+          if (compute_exactly_targets.size() > 0) {
+            arma::fmat result = swapTarget(
+                    data,
+                    distMat,
+                    medoidIndices,
+                    &compute_exactly_targets,
+                    &bestDistances,
+                    &secondBestDistances,
+                    assignments,
+                    (true ? N > 0 : false));
+
+            // result will be k x T
+            // Now update the correct indices
+            estimates.cols(compute_exactly_targets) = result;
+            ucbs.cols(compute_exactly_targets) = result;
+            lcbs.cols(compute_exactly_targets) = result;
+            exactMask.cols(compute_exactly_targets).fill(1);
+            numSamples.cols(compute_exactly_targets) += N;
+            candidates = (lcbs < ucbs.min()) && (exactMask == 0);
+          }
+          if (arma::accu(candidates) < precision) {
+            break;
+          }
+
+          // candidate_targets should be of size T
+          // Sum the different columns
+          // if any index appears in at least one column, sample it
+          arma::uvec candidate_targets = arma::find(
+                  arma::sum(candidates, 0) >= 1);
+
+          // result will be k x T
+          arma::fmat result = swapTarget(
+                  data,
+                  distMat,
+                  medoidIndices,
+                  &candidate_targets,
+                  &bestDistances,
+                  &secondBestDistances,
+                  assignments,
+                  false);
+
+          // candidate_targets should be of size T, 1
+          estimates.cols(candidate_targets) =
+                  ((numSamples.cols(candidate_targets)
+                    % estimates.cols(candidate_targets))
+                   + (result * batchSize)) / (batchSize +
+                                              numSamples.cols(
+                                                      candidate_targets));
+
+          // numSamples should be k x N
+          // select the T of N columns that are candidates
+          numSamples.cols(candidate_targets) += batchSize;
+
+          arma::fmat adjust(nMedoids, candidate_targets.size());
+          // TODO(@motiwari): Move this ::fill to the previous line
+          adjust.fill(p);
+          // Assume swapConfidence is given in logspace
+          adjust = swapConfidence + arma::log(adjust);
+          arma::fmat confBoundDelta = sigma.cols(candidate_targets) %
+                                      arma::sqrt(adjust / numSamples.cols(
+                                              candidate_targets));
+          ucbs.cols(candidate_targets) = estimates.cols(candidate_targets)
+                                         + confBoundDelta;
+          lcbs.cols(candidate_targets) = estimates.cols(candidate_targets)
+                                         - confBoundDelta;
+
+          candidates = (lcbs < ucbs.min()) && (exactMask == 0);
       }
 
-      // candidate_targets should be of size T
-      // Sum the different columns
-      // if any index appears in at least one column, sample it
-      arma::uvec candidate_targets = arma::find(arma::sum(candidates, 0) >= 1);
-
-      // result will be k x T
-      arma::fmat result = swapTarget(
-        data,
-        distMat,
-        medoidIndices,
-        &candidate_targets,
-        &bestDistances,
-        &secondBestDistances,
-        assignments,
-        0);
-
-      // candidate_targets should be of size T, 1
-      estimates.cols(candidate_targets) =
-        ((numSamples.cols(candidate_targets)
-        % estimates.cols(candidate_targets))
-        + (result * batchSize)) / (batchSize +
-          numSamples.cols(candidate_targets));
-
-      // numSamples should be k x N
-      // select the T of N columns that are candidates
-      numSamples.cols(candidate_targets) += batchSize;
-
-      arma::fmat adjust(nMedoids, candidate_targets.size());
-      // TODO(@motiwari): Move this ::fill to the previous line
-      adjust.fill(p);
-      // Assume swapConfidence is given in logspace
-      adjust = swapConfidence + arma::log(adjust);
-      arma::fmat confBoundDelta = sigma.cols(candidate_targets) %
-        arma::sqrt(adjust / numSamples.cols(candidate_targets));
-      ucbs.cols(candidate_targets) = estimates.cols(candidate_targets)
-        + confBoundDelta;
-      lcbs.cols(candidate_targets) = estimates.cols(candidate_targets)
-        - confBoundDelta;
-
-      candidates = (lcbs < ucbs.min()) && (exactMask == 0);
-    }
+      // Perform the medoid switch
+      arma::uword newMedoid = lcbs.index_min();
+      size_t k = newMedoid % nMedoids;
+      size_t n = newMedoid / nMedoids;
+      swapPerformed = (*medoidIndices)(k) != n;
 
-    // Perform the medoid switch
-    arma::uword newMedoid = lcbs.index_min();
-    size_t k = newMedoid % nMedoids;
-    size_t n = newMedoid / nMedoids;
-    swapPerformed = (*medoidIndices)(k) != n;
+      if (swapPerformed) {
+        (*medoidIndices)(k) = n;
+        medoids->col(k) = data.col((*medoidIndices)(k));
+      }
 
-    if (swapPerformed) {
-      (*medoidIndices)(k) = n;
-      medoids->col(k) = data.col((*medoidIndices)(k));
+      calcBestDistancesSwap(
+              data,
+              distMat,
+              medoidIndices,
+              &bestDistances,
+              &secondBestDistances,
+              assignments,
+              swapPerformed);
     }
-
-    calcBestDistancesSwap(
-      data,
-      distMat,
-      medoidIndices,
-      &bestDistances,
-      &secondBestDistances,
-      assignments,
-      swapPerformed);
   }
-}
 }  // namespace km
diff --git a/src/algorithms/banditpam_orig.cpp b/src/algorithms/banditpam_orig.cpp
index 611d3520..61b76016 100644
--- a/src/algorithms/banditpam_orig.cpp
+++ b/src/algorithms/banditpam_orig.cpp
@@ -13,501 +13,508 @@
 #include <vector>
 
 namespace km {
-void BanditPAM_orig::fitBanditPAM_orig(
-  const arma::fmat& inputData,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat) {
-  data = arma::trans(inputData);
-
-  // Note: even if we are using a distance matrix, we compute the permutation
-  // in the block below because it is used elsewhere in the call stack
-  // TODO(@motiwari): Remove need for data or permutation through when using
-  //  a distance matrix
-  if (this->useCache) {
-    size_t n = data.n_cols;
-    size_t m = fmin(n, cacheWidth);
-    cache = new float[n * m];
-
-    #pragma omp parallel for if (this->parallelize)
-    for (size_t idx = 0; idx < m*n; idx++) {
-      cache[idx] = -1;  // TODO(@motiwari): need better value here
-    }
+  void BanditPAM_orig::fitBanditPAM_orig(
+          const arma::fmat &inputData,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat) {
+    data = arma::trans(inputData);
+
+    // Note: even if we are using a distance matrix,
+    // we compute the permutation
+    // in the block below because it is used elsewhere in the call stack
+    // TODO(@motiwari): Remove need for data or permutation
+    //  through when using a distance matrix
+    if (this->useCache) {
+      size_t n = data.n_cols;
+      size_t m = fmin(n, cacheWidth);
+      cache = new float[n * m];
+
+      #pragma omp parallel for if (this->parallelize)
+      for (size_t idx = 0; idx < m * n; idx++) {
+        cache[idx] = -1;  // TODO(@motiwari): need better value here
+      }
 
-    permutation = arma::randperm(n);
-    permutationIdx = 0;
-    reindex = {};  // TODO(@motiwari): Can this intialization be removed?
-    // TODO(@motiwari): Can we parallelize this?
-    for (size_t counter = 0; counter < m; counter++) {
-      reindex[permutation[counter]] = counter;
-    }
-}
-
-  arma::fmat medoidMatrix(data.n_rows, nMedoids);
-  arma::urowvec medoidIndices(nMedoids);
-  steps = 0;
-  BanditPAM_orig::build(data, distMat, &medoidIndices, &medoidMatrix);
-
-  medoidIndicesBuild = medoidIndices;
-  arma::urowvec assignments(data.n_cols);
-
-  BanditPAM_orig::swap(
-      data,
-      distMat,
-      &medoidIndices,
-      &medoidMatrix,
-      &assignments);
-
-  medoidIndicesFinal = medoidIndices;
-  labels = assignments;
-}
-
-arma::frowvec BanditPAM_orig::buildSigma(
-  const arma::fmat& data,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-  const arma::frowvec& bestDistances,
-  const bool useAbsolute) {
-  size_t N = data.n_cols;
-  arma::uvec referencePoints;
-  // TODO(@motiwari): Make this wraparound properly as
-  //  last batch_size elements are dropped
-  if (usePerm) {
-    if ((permutationIdx + batchSize - 1) >= N) {
+      permutation = arma::randperm(n);
       permutationIdx = 0;
-    }
-    // inclusive of both indices
-    referencePoints = permutation.subvec(
-      permutationIdx,
-      permutationIdx + batchSize - 1);
-    permutationIdx += batchSize;
-  } else {
-    referencePoints = arma::randperm(N, batchSize);
-  }
-
-  arma::fvec sample(batchSize);
-  arma::frowvec updated_sigma(N);
-  #pragma omp parallel for if (this->parallelize)
-  for (size_t i = 0; i < N; i++) {
-    for (size_t j = 0; j < batchSize; j++) {
-      // 0 for MISC
-      float cost =
-          KMedoids::cachedLoss(data, distMat, i, referencePoints(j), 0);
-      if (useAbsolute) {
-        sample(j) = cost;
-      } else {
-        sample(j) = cost < bestDistances(referencePoints(j))
-                          ? cost : bestDistances(referencePoints(j));
-        sample(j) -= bestDistances(referencePoints(j));
+      reindex = {};  // TODO(@motiwari): Can this be removed?
+      // TODO(@motiwari): Can we parallelize this?
+      for (size_t counter = 0; counter < m; counter++) {
+        reindex[permutation[counter]] = counter;
       }
     }
-    updated_sigma(i) = arma::stddev(sample);
-  }
-  return updated_sigma;
-}
-
-arma::frowvec BanditPAM_orig::buildTarget(
-  const arma::fmat& data,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-  const arma::uvec* target,
-  const arma::frowvec* bestDistances,
-  const bool useAbsolute,
-  const size_t exact = 0) {
-  size_t N = data.n_cols;
-  size_t tmpBatchSize = batchSize;
-  if (exact > 0) {
-    tmpBatchSize = N;
-  }
-  arma::frowvec estimates(target->n_rows, arma::fill::zeros);
-  arma::uvec referencePoints;
-  // TODO(@motiwari): Make this wraparound properly
-  //  as last batch_size elements are dropped
-  if (usePerm) {
-    if ((permutationIdx + tmpBatchSize - 1) >= N) {
-      permutationIdx = 0;
-    }
-    // inclusive of both indices
-    referencePoints = permutation.subvec(
-      permutationIdx,
-      permutationIdx + tmpBatchSize - 1);
-    permutationIdx += tmpBatchSize;
-  } else {
-    referencePoints = arma::randperm(N, tmpBatchSize);
-  }
 
-  #pragma omp parallel for if (this->parallelize)
-  for (size_t i = 0; i < target->n_rows; i++) {
-    float total = 0;
-    for (size_t j = 0; j < referencePoints.n_rows; j++) {
-      float cost =
-          KMedoids::cachedLoss(
+    arma::fmat medoidMatrix(data.n_rows, nMedoids);
+    arma::urowvec medoidIndices(nMedoids);
+    steps = 0;
+    BanditPAM_orig::build(data, distMat, &medoidIndices, &medoidMatrix);
+
+    medoidIndicesBuild = medoidIndices;
+    arma::urowvec assignments(data.n_cols);
+
+    if (nMedoids > 1) {
+      BanditPAM_orig::swap(
               data,
               distMat,
-              (*target)(i),
-              referencePoints(j),
-              1);  // 1 for BUILD
-      if (useAbsolute) {
-        total += cost;
-      } else {
-        total += cost < (*bestDistances)(referencePoints(j))
-                      ? cost : (*bestDistances)(referencePoints(j));
-        total -= (*bestDistances)(referencePoints(j));
-      }
+              &medoidIndices,
+              &medoidMatrix,
+              &assignments);
     }
-    estimates(i) = total / tmpBatchSize;
+
+    medoidIndicesFinal = medoidIndices;
+    labels = assignments;
   }
-  return estimates;
-}
-
-void BanditPAM_orig::build(
-  const arma::fmat& data,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-  arma::urowvec* medoidIndices,
-  arma::fmat* medoids) {
-  size_t N = data.n_cols;
-  arma::frowvec N_mat(N);
-  N_mat.fill(N);
-  size_t p = N;
-  bool useAbsolute = true;
-  arma::frowvec estimates(N, arma::fill::zeros);
-  arma::frowvec bestDistances(N);
-  bestDistances.fill(std::numeric_limits<float>::infinity());
-  arma::frowvec sigma(N);
-  arma::urowvec candidates(N, arma::fill::ones);
-  arma::frowvec lcbs(N);
-  arma::frowvec ucbs(N);
-  arma::frowvec numSamples(N, arma::fill::zeros);
-  arma::frowvec exactMask(N, arma::fill::zeros);
-
-  // TODO(@motiwari): #pragma omp parallel for if (this->parallelize)?
-  for (size_t k = 0; k < nMedoids; k++) {
-    // instantiate medoids one-by-one
-    permutationIdx = 0;
-    candidates.fill(1);
-    numSamples.fill(0);
-    exactMask.fill(0);
-    estimates.fill(0);
-    // compute std dev amongst batch of reference points
-    sigma = buildSigma(data, distMat, bestDistances, useAbsolute);
-
-    while (arma::sum(candidates) > precision) {
-      // TODO(@motiwari): Do not need a matrix for this comparison,
-      //  use broadcasting
-      arma::umat compute_exactly =
-          ((numSamples + batchSize) >= N_mat) != exactMask;
-      if (arma::accu(compute_exactly) > 0) {
-        arma::uvec targets = find(compute_exactly);
-        arma::frowvec result = buildTarget(
-          data,
-          distMat,
-          &targets,
-          &bestDistances,
-          useAbsolute,
-          N);
-        estimates.cols(targets) = result;
-        ucbs.cols(targets) = result;
-        lcbs.cols(targets) = result;
-        exactMask.cols(targets).fill(1);
-        numSamples.cols(targets) += N;
-        candidates.cols(targets).fill(0);
-      }
-      if (arma::sum(candidates) < precision) {
-        break;
+
+  arma::frowvec BanditPAM_orig::buildSigma(
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const arma::frowvec &bestDistances,
+          const bool useAbsolute) {
+    size_t N = data.n_cols;
+    arma::uvec referencePoints;
+    // TODO(@motiwari): Make this wraparound properly as
+    //  last batch_size elements are dropped
+    if (usePerm) {
+      if ((permutationIdx + batchSize - 1) >= N) {
+        permutationIdx = 0;
       }
-      arma::uvec targets = arma::find(candidates);
-      arma::frowvec result = buildTarget(
-        data,
-        distMat,
-        &targets,
-        &bestDistances,
-        useAbsolute,
-        0);
-      // update the running average
-      estimates.cols(targets) =
-        ((numSamples.cols(targets) % estimates.cols(targets)) +
-         (result * batchSize)) /
-        (batchSize + numSamples.cols(targets));
-      numSamples.cols(targets) += batchSize;
-      arma::frowvec adjust(targets.n_rows);
-      adjust.fill(p);
-      // Assume buildConfidence is given in logspace
-      adjust = buildConfidence + arma::log(adjust);
-      arma::frowvec confBoundDelta =
-        sigma.cols(targets) %
-        arma::sqrt(adjust / numSamples.cols(targets));
-      ucbs.cols(targets) = estimates.cols(targets) + confBoundDelta;
-      lcbs.cols(targets) = estimates.cols(targets) - confBoundDelta;
-      candidates = (lcbs < ucbs.min()) && (exactMask == 0);
+      // inclusive of both indices
+      referencePoints = permutation.subvec(
+              permutationIdx,
+              permutationIdx + batchSize - 1);
+      permutationIdx += batchSize;
+    } else {
+      referencePoints = arma::randperm(N, batchSize);
     }
 
-    medoidIndices->at(k) = lcbs.index_min();
-    medoids->unsafe_col(k) = data.unsafe_col((*medoidIndices)(k));
-
-    // don't need to do this on final iteration
+    arma::fvec sample(batchSize);
+    arma::frowvec updated_sigma(N);
     #pragma omp parallel for if (this->parallelize)
     for (size_t i = 0; i < N; i++) {
-      float cost = KMedoids::cachedLoss(
-        data,
-        distMat,
-        i,
-        (*medoidIndices)(k),
-        0);  // 0 for MISC
-      if (cost < bestDistances(i)) {
-        bestDistances(i) = cost;
+      for (size_t j = 0; j < batchSize; j++) {
+        // 0 for MISC
+        float cost =
+                KMedoids::cachedLoss(data, distMat, i,
+                                     referencePoints(j), 0);
+        if (useAbsolute) {
+          sample(j) = cost;
+        } else {
+          sample(j) = cost < bestDistances(referencePoints(j))
+                        ? cost : bestDistances(referencePoints(j));
+          sample(j) -= bestDistances(referencePoints(j));
+        }
       }
+      updated_sigma(i) = arma::stddev(sample);
     }
-    // use difference of loss for sigma and sampling, not absolute
-    useAbsolute = false;
+    return updated_sigma;
   }
-}
-
-arma::fmat BanditPAM_orig::swapSigma(
-  const arma::fmat& data,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-  const arma::frowvec* bestDistances,
-  const arma::frowvec* secondBestDistances,
-  const arma::urowvec* assignments) {
-  size_t N = data.n_cols;
-  size_t K = nMedoids;
-  arma::fmat updated_sigma(K, N, arma::fill::zeros);
-  arma::uvec referencePoints;
-  // TODO(@motiwari): Make this wraparound properly
-  //  as last batch_size elements are dropped
-  if (usePerm) {
-    if ((permutationIdx + batchSize - 1) >= N) {
-      permutationIdx = 0;
+
+  arma::frowvec BanditPAM_orig::buildTarget(
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const arma::uvec *target,
+          const arma::frowvec *bestDistances,
+          const bool useAbsolute,
+          const bool exact = false) {
+    size_t N = data.n_cols;
+    size_t tmpBatchSize = batchSize;
+    if (exact) {
+      tmpBatchSize = N;
+    }
+    arma::frowvec estimates(target->n_rows, arma::fill::zeros);
+    arma::uvec referencePoints;
+    // TODO(@motiwari): Make this wraparound properly
+    //  as last batch_size elements are dropped
+    if (usePerm) {
+      if ((permutationIdx + tmpBatchSize - 1) >= N) {
+        permutationIdx = 0;
+      }
+      // inclusive of both indices
+      referencePoints = permutation.subvec(
+              permutationIdx,
+              permutationIdx + tmpBatchSize - 1);
+      permutationIdx += tmpBatchSize;
+    } else {
+      referencePoints = arma::randperm(N, tmpBatchSize);
     }
-    // inclusive of both indices
-    referencePoints = permutation.subvec(
-      permutationIdx,
-      permutationIdx + batchSize - 1);
-    permutationIdx += batchSize;
-  } else {
-    referencePoints = arma::randperm(N, batchSize);
-  }
 
-  arma::fvec sample(batchSize);
-  // for each considered swap
-  #pragma omp parallel for if (this->parallelize)
-  for (size_t i = 0; i < K * N; i++) {
-    // extract data point of swap
-    size_t n = i / K;
-    size_t k = i % K;
-
-    // calculate change in loss for some subset of the data
-    for (size_t j = 0; j < batchSize; j++) {
-      // 0 for MISC when estimating sigma
-      float cost =
-          KMedoids::cachedLoss(data, distMat, n, referencePoints(j), 0);
-
-      if (k == (*assignments)(referencePoints(j))) {
-        if (cost < (*secondBestDistances)(referencePoints(j))) {
-          sample(j) = cost;
-        } else {
-          sample(j) = (*secondBestDistances)(referencePoints(j));
-        }
-      } else {
-        if (cost < (*bestDistances)(referencePoints(j))) {
-          sample(j) = cost;
+    #pragma omp parallel for if (this->parallelize)
+    for (size_t i = 0; i < target->n_rows; i++) {
+      float total = 0;
+      for (size_t j = 0; j < referencePoints.n_rows; j++) {
+        float cost =
+                KMedoids::cachedLoss(
+                        data,
+                        distMat,
+                        (*target)(i),
+                        referencePoints(j),
+                        1);  // 1 for BUILD
+        if (useAbsolute) {
+          total += cost;
         } else {
-          sample(j) = (*bestDistances)(referencePoints(j));
+          total += cost < (*bestDistances)(referencePoints(j))
+                     ? cost : (*bestDistances)(referencePoints(j));
+          total -= (*bestDistances)(referencePoints(j));
         }
       }
-      sample(j) -= (*bestDistances)(referencePoints(j));
+      estimates(i) = total / tmpBatchSize;
     }
-    updated_sigma(k, n) = arma::stddev(sample);
-  }
-  return updated_sigma;
-}
-
-arma::fvec BanditPAM_orig::swapTarget(
-  const arma::fmat& data,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-  const arma::urowvec* medoidIndices,
-  const arma::uvec* targets,
-  const arma::frowvec* bestDistances,
-  const arma::frowvec* secondBestDistances,
-  const arma::urowvec* assignments,
-  const size_t exact = 0) {
-  size_t N = data.n_cols;
-  arma::fvec estimates(targets->n_rows, arma::fill::zeros);
-
-  size_t tmpBatchSize = batchSize;
-  if (exact > 0) {
-    tmpBatchSize = N;
+    return estimates;
   }
 
-  arma::uvec referencePoints;
-  // TODO(@motiwari): Make this wraparound properly
-  //  as last batch_size elements are dropped
-  // TODO(@motiwari): Break this duplicated code into a function
-  if (usePerm) {
-    if ((permutationIdx + tmpBatchSize - 1) >= N) {
+  void BanditPAM_orig::build(
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          arma::urowvec *medoidIndices,
+          arma::fmat *medoids) {
+    size_t N = data.n_cols;
+    arma::frowvec N_mat(N);
+    N_mat.fill(N);
+    size_t p = N;
+    bool useAbsolute = true;
+    arma::frowvec estimates(N, arma::fill::zeros);
+    arma::frowvec bestDistances(N);
+    bestDistances.fill(std::numeric_limits<float>::infinity());
+    arma::frowvec sigma(N);
+    arma::urowvec candidates(N, arma::fill::ones);
+    arma::frowvec lcbs(N);
+    arma::frowvec ucbs(N);
+    arma::frowvec numSamples(N, arma::fill::zeros);
+    arma::frowvec exactMask(N, arma::fill::zeros);
+
+    // TODO(@motiwari): #pragma omp parallel for if (this->parallelize)?
+    for (size_t k = 0; k < nMedoids; k++) {
+      // instantiate medoids one-by-one
       permutationIdx = 0;
+      candidates.fill(1);
+      numSamples.fill(0);
+      exactMask.fill(0);
+      estimates.fill(0);
+      // compute std dev amongst batch of reference points
+      sigma = buildSigma(data, distMat, bestDistances, useAbsolute);
+
+      while (arma::sum(candidates) > precision) {
+        // TODO(@motiwari): Do not need a matrix for this comparison,
+        //  use broadcasting
+        arma::umat compute_exactly =
+                ((numSamples + batchSize) >= N_mat) != exactMask;
+        if (arma::accu(compute_exactly) > 0) {
+          arma::uvec targets = find(compute_exactly);
+          arma::frowvec result = buildTarget(
+                  data,
+                  distMat,
+                  &targets,
+                  &bestDistances,
+                  useAbsolute,
+                  (true ? N > 0 : false));
+          estimates.cols(targets) = result;
+          ucbs.cols(targets) = result;
+          lcbs.cols(targets) = result;
+          exactMask.cols(targets).fill(1);
+          numSamples.cols(targets) += N;
+          candidates.cols(targets).fill(0);
+        }
+        if (arma::sum(candidates) < precision) {
+          break;
+        }
+        arma::uvec targets = arma::find(candidates);
+        arma::frowvec result = buildTarget(
+                data,
+                distMat,
+                &targets,
+                &bestDistances,
+                useAbsolute,
+                false);
+        // update the running average
+        estimates.cols(targets) =
+                ((numSamples.cols(targets) % estimates.cols(targets)) +
+                 (result * batchSize)) /
+                (batchSize + numSamples.cols(targets));
+        numSamples.cols(targets) += batchSize;
+        arma::frowvec adjust(targets.n_rows);
+        adjust.fill(p);
+        // Assume buildConfidence is given in logspace
+        adjust = buildConfidence + arma::log(adjust);
+        arma::frowvec confBoundDelta =
+                sigma.cols(targets) %
+                arma::sqrt(adjust / numSamples.cols(targets));
+        ucbs.cols(targets) = estimates.cols(targets) + confBoundDelta;
+        lcbs.cols(targets) = estimates.cols(targets) - confBoundDelta;
+        candidates = (lcbs < ucbs.min()) && (exactMask == 0);
+      }
+
+      medoidIndices->at(k) = lcbs.index_min();
+      medoids->unsafe_col(k) = data.unsafe_col((*medoidIndices)(k));
+
+      // don't need to do this on final iteration
+      #pragma omp parallel for if (this->parallelize)
+      for (size_t i = 0; i < N; i++) {
+        float cost = KMedoids::cachedLoss(
+                data,
+                distMat,
+                i,
+                (*medoidIndices)(k),
+                0);  // 0 for MISC
+        if (cost < bestDistances(i)) {
+          bestDistances(i) = cost;
+        }
+      }
+      // use difference of loss for sigma and sampling, not absolute
+      useAbsolute = false;
     }
-    // inclusive of both indices
-    referencePoints = permutation.subvec(
-      permutationIdx,
-      permutationIdx + tmpBatchSize - 1);
-    permutationIdx += tmpBatchSize;
-  } else {
-    referencePoints = arma::randperm(N, tmpBatchSize);
   }
 
-  // TODO(@motiwari): Declare variables outside of loops
-  #pragma omp parallel for if (this->parallelize)
-  for (size_t i = 0; i < targets->n_rows; i++) {
-    float total = 0;
-    // extract data point of swap
-    size_t n = (*targets)(i) / medoidIndices->n_cols;
-    size_t k = (*targets)(i) % medoidIndices->n_cols;
-    // calculate total loss for some subset of the data
-    for (size_t j = 0; j < tmpBatchSize; j++) {
-      // 2 for SWAP
-      float cost =
-          KMedoids::cachedLoss(data, distMat, n, referencePoints(j), 2);
-      if (k == (*assignments)(referencePoints(j))) {
-        if (cost < (*secondBestDistances)(referencePoints(j))) {
-          total += cost;
+  arma::fmat BanditPAM_orig::swapSigma(
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const arma::frowvec *bestDistances,
+          const arma::frowvec *secondBestDistances,
+          const arma::urowvec *assignments) {
+    size_t N = data.n_cols;
+    size_t K = nMedoids;
+    arma::fmat updated_sigma(K, N, arma::fill::zeros);
+    arma::uvec referencePoints;
+    // TODO(@motiwari): Make this wraparound properly
+    //  as last batch_size elements are dropped
+    if (usePerm) {
+      if ((permutationIdx + batchSize - 1) >= N) {
+        permutationIdx = 0;
+      }
+      // inclusive of both indices
+      referencePoints = permutation.subvec(
+              permutationIdx,
+              permutationIdx + batchSize - 1);
+      permutationIdx += batchSize;
+    } else {
+      referencePoints = arma::randperm(N, batchSize);
+    }
+
+    arma::fvec sample(batchSize);
+    // for each considered swap
+    #pragma omp parallel for if (this->parallelize)
+    for (size_t i = 0; i < K * N; i++) {
+      // extract data point of swap
+      size_t n = i / K;
+      size_t k = i % K;
+
+      // calculate change in loss for some subset of the data
+      for (size_t j = 0; j < batchSize; j++) {
+        // 0 for MISC when estimating sigma
+        float cost =
+                KMedoids::cachedLoss(data, distMat, n,
+                                     referencePoints(j), 0);
+
+        if (k == (*assignments)(referencePoints(j))) {
+          if (cost < (*secondBestDistances)(referencePoints(j))) {
+            sample(j) = cost;
+          } else {
+            sample(j) = (*secondBestDistances)(referencePoints(j));
+          }
         } else {
-          total += (*secondBestDistances)(referencePoints(j));
+          if (cost < (*bestDistances)(referencePoints(j))) {
+            sample(j) = cost;
+          } else {
+            sample(j) = (*bestDistances)(referencePoints(j));
+          }
         }
-      } else {
-        if (cost < (*bestDistances)(referencePoints(j))) {
-          total += cost;
+        sample(j) -= (*bestDistances)(referencePoints(j));
+      }
+      updated_sigma(k, n) = arma::stddev(sample);
+    }
+    return updated_sigma;
+  }
+
+  arma::fvec BanditPAM_orig::swapTarget(
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const arma::urowvec *medoidIndices,
+          const arma::uvec *targets,
+          const arma::frowvec *bestDistances,
+          const arma::frowvec *secondBestDistances,
+          const arma::urowvec *assignments,
+          const bool exact = false) {
+    size_t N = data.n_cols;
+    arma::fvec estimates(targets->n_rows, arma::fill::zeros);
+
+    size_t tmpBatchSize = batchSize;
+    if (exact) {
+      tmpBatchSize = N;
+    }
+
+    arma::uvec referencePoints;
+    // TODO(@motiwari): Make this wraparound properly
+    //  as last batch_size elements are dropped
+    // TODO(@motiwari): Break this duplicated code into a function
+    if (usePerm) {
+      if ((permutationIdx + tmpBatchSize - 1) >= N) {
+        permutationIdx = 0;
+      }
+      // inclusive of both indices
+      referencePoints = permutation.subvec(
+              permutationIdx,
+              permutationIdx + tmpBatchSize - 1);
+      permutationIdx += tmpBatchSize;
+    } else {
+      referencePoints = arma::randperm(N, tmpBatchSize);
+    }
+
+    // TODO(@motiwari): Declare variables outside of loops
+    #pragma omp parallel for if (this->parallelize)
+    for (size_t i = 0; i < targets->n_rows; i++) {
+      float total = 0;
+      // extract data point of swap
+      size_t n = (*targets)(i) / medoidIndices->n_cols;
+      size_t k = (*targets)(i) % medoidIndices->n_cols;
+      // calculate total loss for some subset of the data
+      for (size_t j = 0; j < tmpBatchSize; j++) {
+        // 2 for SWAP
+        float cost =
+                KMedoids::cachedLoss(data, distMat, n,
+                                     referencePoints(j), 2);
+        if (k == (*assignments)(referencePoints(j))) {
+          if (cost < (*secondBestDistances)(referencePoints(j))) {
+            total += cost;
+          } else {
+            total += (*secondBestDistances)(referencePoints(j));
+          }
         } else {
-          total += (*bestDistances)(referencePoints(j));
+          if (cost < (*bestDistances)(referencePoints(j))) {
+            total += cost;
+          } else {
+            total += (*bestDistances)(referencePoints(j));
+          }
         }
+        total -= (*bestDistances)(referencePoints(j));
       }
-      total -= (*bestDistances)(referencePoints(j));
+      // TODO(@motiwari): we can probably avoid this division
+      //  if we look at total loss, not average loss
+      estimates(i) = total / referencePoints.n_rows;
     }
-    // TODO(@motiwari): we can probably avoid this division
-    //  if we look at total loss, not average loss
-    estimates(i) = total / referencePoints.n_rows;
+    return estimates;
   }
-  return estimates;
-}
-
-void BanditPAM_orig::swap(
-  const arma::fmat& data,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-  arma::urowvec* medoidIndices,
-  arma::fmat* medoids,
-  arma::urowvec* assignments) {
-  size_t N = data.n_cols;
-  size_t p = N;
-
-  arma::fmat sigma(nMedoids, N, arma::fill::zeros);
-
-  arma::frowvec bestDistances(N);
-  arma::frowvec secondBestDistances(N);
-  bool swapPerformed = true;
-  arma::umat candidates(nMedoids, N, arma::fill::ones);
-  arma::umat exactMask(nMedoids, N, arma::fill::zeros);
-  arma::fmat estimates(nMedoids, N, arma::fill::zeros);
-  arma::fmat lcbs(nMedoids, N);
-  arma::fmat ucbs(nMedoids, N);
-  arma::umat numSamples(nMedoids, N, arma::fill::zeros);
-
-  // calculate quantities needed for swap, bestDistances and sigma
-  calcBestDistancesSwap(
-    data,
-    distMat,
-    medoidIndices,
-    &bestDistances,
-    &secondBestDistances,
-    assignments,
-    swapPerformed);
-
-  // continue making swaps while loss is decreasing
-  while (swapPerformed && steps < maxIter) {
-    steps++;
-    permutationIdx = 0;
-
-    sigma = swapSigma(
-      data,
-      distMat,
-      &bestDistances,
-      &secondBestDistances,
-      assignments);
-
-    // Reset variables when starting a new swap
-    candidates.fill(1);
-    exactMask.fill(0);
-    estimates.fill(0);
-    numSamples.fill(0);
-
-    // while there is at least one candidate (float comparison issues)
-    while (arma::accu(candidates) > 1.5) {
-      // compute exactly if it's been samples more than N times and
-      // hasn't been computed exactly already
-      arma::umat compute_exactly =
-          ((numSamples + batchSize) >= N) != (exactMask);
-      arma::uvec targets = arma::find(compute_exactly);
-
-      if (targets.size() > 0) {
+
+  void BanditPAM_orig::swap(
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          arma::urowvec *medoidIndices,
+          arma::fmat *medoids,
+          arma::urowvec *assignments) {
+    size_t N = data.n_cols;
+    size_t p = N;
+
+    arma::fmat sigma(nMedoids, N, arma::fill::zeros);
+
+    arma::frowvec bestDistances(N);
+    arma::frowvec secondBestDistances(N);
+    bool swapPerformed = true;
+    arma::umat candidates(nMedoids, N, arma::fill::ones);
+    arma::umat exactMask(nMedoids, N, arma::fill::zeros);
+    arma::fmat estimates(nMedoids, N, arma::fill::zeros);
+    arma::fmat lcbs(nMedoids, N);
+    arma::fmat ucbs(nMedoids, N);
+    arma::umat numSamples(nMedoids, N, arma::fill::zeros);
+
+    // calculate quantities needed for swap, bestDistances and sigma
+    calcBestDistancesSwap(
+            data,
+            distMat,
+            medoidIndices,
+            &bestDistances,
+            &secondBestDistances,
+            assignments,
+            swapPerformed);
+
+    // continue making swaps while loss is decreasing
+    while (swapPerformed && steps < maxIter) {
+      steps++;
+      permutationIdx = 0;
+
+      sigma = swapSigma(
+              data,
+              distMat,
+              &bestDistances,
+              &secondBestDistances,
+              assignments);
+
+      // Reset variables when starting a new swap
+      candidates.fill(1);
+      exactMask.fill(0);
+      estimates.fill(0);
+      numSamples.fill(0);
+
+      // while there is at least one candidate (float comparison issues)
+      while (arma::accu(candidates) > 1.5) {
+        // compute exactly if it's been samples more than N times and
+        // hasn't been computed exactly already
+        arma::umat compute_exactly =
+                ((numSamples + batchSize) >= N) != (exactMask);
+        arma::uvec targets = arma::find(compute_exactly);
+
+        if (targets.size() > 0) {
+          arma::fvec result = swapTarget(
+                  data,
+                  distMat,
+                  medoidIndices,
+                  &targets,
+                  &bestDistances,
+                  &secondBestDistances,
+                  assignments,
+                  (true ? N > 0 : false));
+          estimates.elem(targets) = result;
+          ucbs.elem(targets) = result;
+          lcbs.elem(targets) = result;
+          exactMask.elem(targets).fill(1);
+          numSamples.elem(targets) += N;
+          candidates = (lcbs < ucbs.min()) && (exactMask == 0);
+        }
+        if (arma::accu(candidates) < precision) {
+          break;
+        }
+        targets = arma::find(candidates);
         arma::fvec result = swapTarget(
-          data,
-          distMat,
-          medoidIndices,
-          &targets,
-          &bestDistances,
-          &secondBestDistances,
-          assignments,
-          N);
-        estimates.elem(targets) = result;
-        ucbs.elem(targets) = result;
-        lcbs.elem(targets) = result;
-        exactMask.elem(targets).fill(1);
-        numSamples.elem(targets) += N;
+                data,
+                distMat,
+                medoidIndices,
+                &targets,
+                &bestDistances,
+                &secondBestDistances,
+                assignments,
+                false);
+        estimates.elem(targets) =
+                ((numSamples.elem(targets) % estimates.elem(targets)) +
+                 (result * batchSize)) /
+                (batchSize + numSamples.elem(targets));
+        numSamples.elem(targets) += batchSize;
+        arma::fvec adjust(targets.n_rows);
+        adjust.fill(p);
+        // Assume swapConfidence is given in logspace
+        adjust = swapConfidence + arma::log(adjust);
+        arma::fvec confBoundDelta = sigma.elem(targets) %
+                                    arma::sqrt(adjust / numSamples.elem(
+                                            targets));
+        ucbs.elem(targets) = estimates.elem(targets) + confBoundDelta;
+        lcbs.elem(targets) = estimates.elem(targets) - confBoundDelta;
         candidates = (lcbs < ucbs.min()) && (exactMask == 0);
       }
-      if (arma::accu(candidates) < precision) {
-        break;
-      }
-      targets = arma::find(candidates);
-      arma::fvec result = swapTarget(
-          data,
-          distMat,
-          medoidIndices,
-          &targets,
-          &bestDistances,
-          &secondBestDistances,
-          assignments,
-          0);
-      estimates.elem(targets) =
-          ((numSamples.elem(targets) % estimates.elem(targets)) +
-           (result * batchSize)) /
-          (batchSize + numSamples.elem(targets));
-      numSamples.elem(targets) += batchSize;
-      arma::fvec adjust(targets.n_rows);
-      adjust.fill(p);
-      // Assume swapConfidence is given in logspace
-      adjust = swapConfidence + arma::log(adjust);
-      arma::fvec confBoundDelta = sigma.elem(targets) %
-                    arma::sqrt(adjust / numSamples.elem(targets));
-      ucbs.elem(targets) = estimates.elem(targets) + confBoundDelta;
-      lcbs.elem(targets) = estimates.elem(targets) - confBoundDelta;
-      candidates = (lcbs < ucbs.min()) && (exactMask == 0);
-    }
 
-    // Perform the medoid switch
-    arma::uword newMedoid = lcbs.index_min();
-    // extract old and new medoids of swap
-    size_t k = newMedoid % medoids->n_cols;
-    size_t n = newMedoid / medoids->n_cols;
-    swapPerformed = (*medoidIndices)(k) != n;
-    steps++;
-
-    if (swapPerformed) {
-      (*medoidIndices)(k) = n;
-      medoids->col(k) = data.col((*medoidIndices)(k));
+      // Perform the medoid switch
+      arma::uword newMedoid = lcbs.index_min();
+      // extract old and new medoids of swap
+      size_t k = newMedoid % medoids->n_cols;
+      size_t n = newMedoid / medoids->n_cols;
+      swapPerformed = (*medoidIndices)(k) != n;
+      steps++;
+
+      if (swapPerformed) {
+        (*medoidIndices)(k) = n;
+        medoids->col(k) = data.col((*medoidIndices)(k));
+      }
+      calcBestDistancesSwap(
+              data,
+              distMat,
+              medoidIndices,
+              &bestDistances,
+              &secondBestDistances,
+              assignments,
+              swapPerformed);
     }
-    calcBestDistancesSwap(
-      data,
-      distMat,
-      medoidIndices,
-      &bestDistances,
-      &secondBestDistances,
-      assignments,
-      swapPerformed);
   }
-}
 }  // namespace km
diff --git a/src/algorithms/fastpam1.cpp b/src/algorithms/fastpam1.cpp
index 48a78975..9fb3f581 100644
--- a/src/algorithms/fastpam1.cpp
+++ b/src/algorithms/fastpam1.cpp
@@ -16,171 +16,193 @@
 #include <unordered_map>
 
 namespace km {
-void FastPAM1::fitFastPAM1(
-  const arma::fmat& inputData,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat) {
-  data = inputData;
-  data = arma::trans(data);
-  arma::urowvec medoidIndices(nMedoids);
-  FastPAM1::buildFastPAM1(data, distMat, &medoidIndices);
-  steps = 0;
-  medoidIndicesBuild = medoidIndices;
-  arma::urowvec assignments(data.n_cols);
-  size_t iter = 0;
-  bool medoidChange = true;
-  while (iter < maxIter && medoidChange) {
-    auto previous{medoidIndices};
-    FastPAM1::swapFastPAM1(data, distMat, &medoidIndices, &assignments);
-    medoidChange = arma::any(medoidIndices != previous);
-    iter++;
+  void FastPAM1::fitFastPAM1(
+          const arma::fmat &inputData,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat) {
+    data = inputData;
+    data = arma::trans(data);
+    arma::urowvec medoidIndices(nMedoids);
+    FastPAM1::buildFastPAM1(data, distMat, &medoidIndices);
+    steps = 0;
+    medoidIndicesBuild = medoidIndices;
+    arma::urowvec assignments(data.n_cols);
+    size_t iter = 0;
+    bool medoidChange = true;
+    while (iter < maxIter && medoidChange) {
+      auto previous{medoidIndices};
+      if (nMedoids > 1) {
+        FastPAM1::swapFastPAM1(data, distMat, &medoidIndices,
+                                 &assignments);
+      }
+
+      medoidChange = arma::any(medoidIndices != previous);
+      iter++;
+    }
+    medoidIndicesFinal = medoidIndices;
+    labels = assignments;
+    steps = iter;
   }
-  medoidIndicesFinal = medoidIndices;
-  labels = assignments;
-  steps = iter;
-}
 
-void FastPAM1::buildFastPAM1(
-  const arma::fmat& data,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-  arma::urowvec* medoidIndices
-) {
-  size_t N = data.n_cols;
-  arma::frowvec estimates(N, arma::fill::zeros);
-  arma::frowvec bestDistances(N);
-  bestDistances.fill(std::numeric_limits<float>::infinity());
-  arma::frowvec sigma(N);
-  float minDistance = std::numeric_limits<float>::infinity();
-  int best = 0;
-  float total = 0;
-  float cost = 0;
+  void FastPAM1::buildFastPAM1(
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          arma::urowvec *medoidIndices
+  ) {
+    size_t N = data.n_cols;
+    arma::frowvec estimates(N, arma::fill::zeros);
+    arma::frowvec bestDistances(N);
+    bestDistances.fill(std::numeric_limits<float>::infinity());
+    arma::frowvec sigma(N);
+    float minDistance = std::numeric_limits<float>::infinity();
+    int best = 0;
+    float total = 0;
+    float cost = 0;
 
-  // TODO(@motiwari): pragma omp parallel for?
-  for (size_t k = 0; k < nMedoids; k++) {
-    minDistance = std::numeric_limits<float>::infinity();
-    best = 0;
-    // fixes a base datapoint
     // TODO(@motiwari): pragma omp parallel for?
-    for (size_t i = 0; i < data.n_cols; i++) {
-      total = 0;
+    for (size_t k = 0; k < nMedoids; k++) {
+      minDistance = std::numeric_limits<float>::infinity();
+      best = 0;
+      // fixes a base datapoint
       // TODO(@motiwari): pragma omp parallel for?
-      for (size_t j = 0; j < data.n_cols; j++) {
-        // computes distance between base and all other points
-        cost = (this->*lossFn)(data, i, j);
-        // compares this with the cached best distance
-        if (bestDistances(j) < cost) {
-          cost = bestDistances(j);
+      for (size_t i = 0; i < data.n_cols; i++) {
+        total = 0;
+        // TODO(@motiwari): pragma omp parallel for?
+        for (size_t j = 0; j < data.n_cols; j++) {
+          // computes distance between base and all other points
+          cost = (this->*lossFn)(data, i, j);
+          // compares this with the cached best distance
+          if (bestDistances(j) < cost) {
+              cost = bestDistances(j);
+          }
+          total += cost;
+        }
+        if (total < minDistance) {
+          minDistance = total;
+          best = i;
         }
-        total += cost;
-      }
-      if (total < minDistance) {
-        minDistance = total;
-        best = i;
       }
-    }
-    (*medoidIndices)(k) = best;
+      (*medoidIndices)(k) = best;
 
-    // update the medoid assignment and best_distance for this datapoint
-    // TODO(@motiwari): pragma omp parallel for?
-    for (size_t l = 0; l < N; l++) {
-      cost = (this->*lossFn)(data, l, (*medoidIndices)(k));
-      if (cost < bestDistances(l)) {
-        bestDistances(l) = cost;
+      // update the medoid assignment and best_distance for this datapoint
+      // TODO(@motiwari): pragma omp parallel for?
+      for (size_t l = 0; l < N; l++) {
+        cost = (this->*lossFn)(data, l, (*medoidIndices)(k));
+        if (cost < bestDistances(l)) {
+          bestDistances(l) = cost;
+        }
       }
     }
   }
-}
-
-void FastPAM1::swapFastPAM1(
-  const arma::fmat& data,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-  arma::urowvec* medoidIndices,
-  arma::urowvec* assignments
-) {
-  float bestChange = 0;
-  size_t swapIn = 0;
-  size_t medoidToSwap = 0;
-  size_t N = data.n_cols;
-  size_t iter = 0;
-  bool swapPerformed = true;
-  arma::fmat sigma(nMedoids, N, arma::fill::zeros);
-  arma::frowvec bestDistances(N);
-  arma::frowvec secondBestDistances(N);
-  arma::frowvec deltaTD(nMedoids, arma::fill::zeros);
 
-  // calculate quantities needed for swap, bestDistances and sigma
-  KMedoids::calcBestDistancesSwap(
-    data,
-    distMat,
-    medoidIndices,
-    &bestDistances,
-    &secondBestDistances,
-    assignments,
-    swapPerformed);
+  void FastPAM1::swapFastPAM1(
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          arma::urowvec *medoidIndices,
+          arma::urowvec *assignments
+  ) {
+    float bestChange = 0;
+    float prevBestChange = 0;
+    size_t swapIn = 0;
+    size_t medoidToSwap = 0;
+    size_t N = data.n_cols;
+    size_t iter = 0;
+    bool swapPerformed = true;
+    arma::fmat sigma(nMedoids, N, arma::fill::zeros);
+    arma::frowvec bestDistances(N);
+    arma::frowvec secondBestDistances(N);
+    arma::frowvec deltaTD(nMedoids, arma::fill::zeros);
 
-  float di = 0;
-  float dij = 0;
+    // calculate quantities needed for swap, bestDistances and sigma
+    KMedoids::calcBestDistancesSwap(
+            data,
+            distMat,
+            medoidIndices,
+            &bestDistances,
+            &secondBestDistances,
+            assignments,
+            swapPerformed);
 
-  while (swapPerformed && iter < maxIter) {
-    iter++;
-    // TODO(@motiwari): pragma omp parallel for?
-    for (size_t i = 0; i < data.n_cols; i++) {
-      di = bestDistances(i);
+    float di = 0;
+    float dij = 0;
 
-      // Consider making point i a medoid.
-      // The total loss then contains at least one term, -di,
-      // because the loss contribution for point i is reduced from di to 0
-      deltaTD.fill(-di);
+    while (swapPerformed && iter < maxIter) {
+      bestChange = 0;
+      deltaTD.zeros();
+      iter++;
       // TODO(@motiwari): pragma omp parallel for?
-      for (size_t j = 0; j < data.n_cols; j++) {
-        if (j != i) {
-          dij = (this->*lossFn)(data, i, j);
-          if (dij < bestDistances(j)) {
-            // Case 1: point i becomes the closest medoid for point j,
-            // regardless of which medoid j was previously assigned to. deltaTD
-            // will be negative across ALL possible medoid indices m
-            deltaTD += (dij -  bestDistances(j));
-          } else if (dij < secondBestDistances(j)) {
-            // Case 2: i. If point i is closer than the second best
-            // medoid but further than the best medoid (enforced by failing
-            // the condition for the above if condition), point i will
-            // become the closest medoid only when we remove its associated
-            // medoid and add point i
-            deltaTD.at((*assignments)(j)) += (dij - bestDistances(j));
-          }  else {
-            // Case 3: dij > secondBestDistances(j). Then the loss for point j
-            // will not change for any medoid swapped out except for its
-            // assignment, in which case it moves to its second nearest medoid
-            deltaTD.at((*assignments)(j)) +=
-              (secondBestDistances(j) - bestDistances(j));
+      for (size_t i = 0; i < data.n_cols; i++) {
+        di = bestDistances(i);
+
+        // Consider making point i a medoid.
+        // The total loss then contains at least one term, -di,
+        // because the loss contribution for point i is
+        // reduced from di to 0
+        deltaTD.fill(-di);
+        // TODO(@motiwari): pragma omp parallel for?
+        for (size_t j = 0; j < data.n_cols; j++) {
+          if (j != i) {
+            dij = (this->*lossFn)(data, i, j);
+            if (dij < bestDistances(j)) {
+              // Case 1: point i becomes the closest
+              // medoid for point j,
+              // regardless of which medoid j was
+              // previously assigned to. deltaTD
+              // will be negative across ALL
+              // possible medoid indices m
+              deltaTD += (dij - bestDistances(j));
+            } else if (dij < secondBestDistances(j)) {
+              // Case 2: i. If point i is closer
+              // than the second best
+              // medoid but further than the best
+              // medoid (enforced by failing
+              // the condition for the above if
+              // condition), point i will
+              // become the closest medoid only
+              // when we remove its associated
+              // medoid and add point i
+              deltaTD.at((*assignments)(j)) += (dij -
+                                                bestDistances(j));
+            } else {
+              // Case 3: dij > secondBestDistances(j).
+              // Then the loss for point j
+              // will not change for any medoid swapped
+              // out except for its
+              // assignment, in which case it moves to
+              // its second nearest medoid
+              deltaTD.at((*assignments)(j)) +=
+                      (secondBestDistances(j) - bestDistances(j));
+            }
           }
         }
-      }
 
-      // Determine the best medoid to swap out
-      arma::uword swapOut = deltaTD.index_min();
-      // If the loss change is better than the best loss change so far,
-      // Update our running best statistics
-      if (deltaTD.min() < bestChange) {
-        bestChange = deltaTD.min();
-        swapIn = i;
-        medoidToSwap = swapOut;
+        // Determine the best medoid to swap out
+        arma::uword swapOut = deltaTD.index_min();
+        // If the loss change is better than
+        // the best loss change so far,
+        // Update our running best statistics
+        if (deltaTD.min() < prevBestChange) {
+          bestChange = deltaTD.min();
+          prevBestChange = bestChange;
+          swapIn = i;
+          medoidToSwap = swapOut;
+        }
       }
-    }
 
-    // Update the loss and perform the swap if the loss would be improved
-    if (bestChange < 0) {
-      (*medoidIndices)(medoidToSwap) = swapIn;
-      calcBestDistancesSwap(
-        data,
-        distMat,
-        medoidIndices,
-        &bestDistances,
-        &secondBestDistances,
-        assignments,
-        swapPerformed);
-    } else {
-      swapPerformed = false;
+      // Update the loss and perform the swap if
+      // the loss would be improved
+      if (bestChange < 0) {
+        (*medoidIndices)(medoidToSwap) = swapIn;
+        calcBestDistancesSwap(
+                data,
+                distMat,
+                medoidIndices,
+                &bestDistances,
+                &secondBestDistances,
+                assignments,
+                swapPerformed);
+      } else {
+        swapPerformed = false;
+      }
     }
   }
-}
 }  // namespace km
diff --git a/src/algorithms/kmedoids_algorithm.cpp b/src/algorithms/kmedoids_algorithm.cpp
index d80bd4c6..ee15fe10 100644
--- a/src/algorithms/kmedoids_algorithm.cpp
+++ b/src/algorithms/kmedoids_algorithm.cpp
@@ -20,422 +20,437 @@ namespace km {
 // NOTE: The order of arguments in this constructor must match that of the
 // arguments in kmedoids_pywrapper.cpp, otherwise undefined behavior can
 // result (variables being initialized with others' values)
-KMedoids::KMedoids(
-  size_t nMedoids,
-  const std::string& algorithm,
-  size_t maxIter,
-  size_t buildConfidence,
-  size_t swapConfidence,
-  bool useCache,
-  bool usePerm,
-  size_t cacheWidth,
-  bool parallelize,
-  size_t seed):
-    nMedoids(nMedoids),
-    algorithm(algorithm),
-    maxIter(maxIter),
-    buildConfidence(buildConfidence),
-    swapConfidence(swapConfidence),
-    useCache(useCache),
-    usePerm(usePerm),
-    cacheWidth(cacheWidth),
-    parallelize(parallelize),
-    seed(seed) {
-  KMedoids::checkAlgorithm(algorithm);
-  // Though we initialize seed from the given parameter,
-  // we need to call setSeed to pass it to arma
-  KMedoids::setSeed(seed);
-}
-
-KMedoids::~KMedoids() {}
-
-void KMedoids::fit(
-  const arma::fmat& inputData,
-  const std::string& loss,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat) {
-  numMiscDistanceComputations = 0;
-  numBuildDistanceComputations = 0;
-  numSwapDistanceComputations = 0;
-  numCacheWrites = 0;
-  numCacheHits = 0;
-  numCacheMisses = 0;
-
-  if (distMat) {  // User has provided a distance matrix
-    if (distMat.value().get().n_cols != distMat.value().get().n_rows) {
-      // TODO(@motiwari): Change this to an assertion that is properly raised
-      throw std::invalid_argument("Malformed distance matrix provided");
+  KMedoids::KMedoids(
+          size_t nMedoids,
+          const std::string &algorithm,
+          size_t maxIter,
+          size_t buildConfidence,
+          size_t swapConfidence,
+          bool useCache,
+          bool usePerm,
+          size_t cacheWidth,
+          bool parallelize,
+          size_t seed) :
+          nMedoids(nMedoids),
+          algorithm(algorithm),
+          maxIter(maxIter),
+          buildConfidence(buildConfidence),
+          swapConfidence(swapConfidence),
+          useCache(useCache),
+          usePerm(usePerm),
+          cacheWidth(cacheWidth),
+          parallelize(parallelize),
+          seed(seed) {
+    KMedoids::checkAlgorithm(algorithm);
+    // Though we initialize seed from the given parameter,
+    // we need to call setSeed to pass it to arma
+    KMedoids::setSeed(seed);
+  }
+
+  KMedoids::~KMedoids() {}
+
+  void KMedoids::fit(
+          const arma::fmat &inputData,
+          const std::string &loss,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat) {
+    numMiscDistanceComputations = 0;
+    numBuildDistanceComputations = 0;
+    numSwapDistanceComputations = 0;
+    numCacheWrites = 0;
+    numCacheHits = 0;
+    numCacheMisses = 0;
+
+    if (distMat) {  // User has provided a distance matrix
+      if (distMat.value().get().n_cols != distMat.value().get().n_rows) {
+          // TODO(@motiwari): Change this to an assertion
+          //  that is properly raised
+          throw std::invalid_argument(
+                  "Malformed distance matrix provided");
+      }
+      useDistMat = true;
+    } else {
+      // In case the user is running a new problem
+      // without a distance matrix
+      // after running a distance matrix problem
+      useDistMat = false;
     }
-    useDistMat = true;
-  } else {
-    // In case the user is running a new problem without a distance matrix
-    // after running a distance matrix problem
-    useDistMat = false;
-  }
-
-  if (inputData.n_rows == 0) {
-    // TODO(@motiwari): Change this to an assertion that is properly raised
-    throw std::invalid_argument("Dataset is empty");
-  }
-  batchSize = fmin(inputData.n_rows, batchSize);
-
-  try {
-    KMedoids::setLossFn(loss);
-    if (algorithm == "PAM") {
-      static_cast<PAM*>(this)->fitPAM(inputData, distMat);
-    } else if (algorithm == "BanditPAM") {
-      static_cast<BanditPAM*>(this)->fitBanditPAM(inputData, distMat);
-    } else if (algorithm == "BanditPAM_orig") {
-      static_cast<BanditPAM_orig*>(this)->fitBanditPAM_orig(inputData, distMat);
-    } else if (algorithm == "FastPAM1") {
-      static_cast<FastPAM1*>(this)->fitFastPAM1(inputData, distMat);
+
+    if (inputData.n_rows == 0) {
+      // TODO(@motiwari): Change this to an assertion
+      //  that is properly raised
+      throw std::invalid_argument("Dataset is empty");
+    }
+    // TODO(@Adarsh321123): assert that the number of medoids is >=
+    //  than the number of points
+    batchSize = fmin(inputData.n_rows, batchSize);
+
+    try {
+      KMedoids::setLossFn(loss);
+      if (algorithm == "PAM") {
+          static_cast<PAM *>(this)->fitPAM(inputData, distMat);
+      } else if (algorithm == "BanditPAM") {
+          static_cast<BanditPAM *>(this)->fitBanditPAM(inputData,
+                                                       distMat);
+      } else if (algorithm == "BanditPAM_orig") {
+          static_cast<BanditPAM_orig *>(this)->fitBanditPAM_orig(
+                  inputData, distMat);
+      } else if (algorithm == "FastPAM1") {
+          static_cast<FastPAM1 *>(this)->fitFastPAM1(inputData, distMat);
+      }
+    } catch (std::invalid_argument &e) {
+      std::cout << e.what() << std::endl;
+      std::cout << "Error: Clustering did not run." << std::endl;
+      throw e;
     }
-  } catch (std::invalid_argument& e) {
-    std::cout << e.what() << std::endl;
-    std::cout << "Error: Clustering did not run." << std::endl;
-    throw e;
   }
-}
 
-arma::urowvec KMedoids::getMedoidsBuild() const {
-  return medoidIndicesBuild;
-}
+  arma::urowvec KMedoids::getMedoidsBuild() const {
+    return medoidIndicesBuild;
+  }
 
-arma::urowvec KMedoids::getMedoidsFinal() const {
-  return medoidIndicesFinal;
-}
+  arma::urowvec KMedoids::getMedoidsFinal() const {
+    return medoidIndicesFinal;
+  }
 
-arma::urowvec KMedoids::getLabels() const {
-  return labels;
-}
+  arma::urowvec KMedoids::getLabels() const {
+    return labels;
+  }
 
-size_t KMedoids::getSteps() const {
-  return steps;
-}
+  size_t KMedoids::getSteps() const {
+    return steps;
+  }
 
-size_t KMedoids::getNMedoids() const {
-  return nMedoids;
-}
+  size_t KMedoids::getNMedoids() const {
+    return nMedoids;
+  }
 
-void KMedoids::setNMedoids(size_t newNMedoids) {
-  nMedoids = newNMedoids;
-}
+  void KMedoids::setNMedoids(size_t newNMedoids) {
+    nMedoids = newNMedoids;
+  }
 
-std::string KMedoids::getAlgorithm() const {
-  return algorithm;
-}
+  std::string KMedoids::getAlgorithm() const {
+    return algorithm;
+  }
 
-void KMedoids::setAlgorithm(const std::string& newAlgorithm) {
-  algorithm = newAlgorithm;
-  KMedoids::checkAlgorithm(algorithm);
-}
+  void KMedoids::setAlgorithm(const std::string &newAlgorithm) {
+    algorithm = newAlgorithm;
+    KMedoids::checkAlgorithm(algorithm);
+  }
 
-size_t KMedoids::getMaxIter() const {
-  return maxIter;
-}
+  size_t KMedoids::getMaxIter() const {
+    return maxIter;
+  }
 
-void KMedoids::setMaxIter(size_t newMaxIter) {
-  maxIter = newMaxIter;
-}
+  void KMedoids::setMaxIter(size_t newMaxIter) {
+    maxIter = newMaxIter;
+  }
 
 
-size_t KMedoids::getBuildConfidence() const {
-  return buildConfidence;
-}
+  size_t KMedoids::getBuildConfidence() const {
+    return buildConfidence;
+  }
 
-void KMedoids::setBuildConfidence(size_t newBuildConfidence) {
-  if (algorithm != "BanditPAM" && algorithm != "BanditPAM_orig") {
-    // TODO(@motiwari): Better error type
-    throw "Cannot set buildConfidence when not using BanditPAM";
+  void KMedoids::setBuildConfidence(size_t newBuildConfidence) {
+    if (algorithm != "BanditPAM" && algorithm != "BanditPAM_orig") {
+        // TODO(@motiwari): Better error type
+        throw "Cannot set buildConfidence when not using BanditPAM";
+    }
+    buildConfidence = newBuildConfidence;
   }
-  buildConfidence = newBuildConfidence;
-}
 
-size_t KMedoids::getSwapConfidence() const {
-  return swapConfidence;
-}
+  size_t KMedoids::getSwapConfidence() const {
+    return swapConfidence;
+  }
 
-void KMedoids::setSwapConfidence(size_t newSwapConfidence) {
-  if (algorithm != "BanditPAM" && algorithm != "BanditPAM_orig") {
-    // TODO(@motiwari): Better error type
-    throw "Cannot set buildConfidence when not using BanditPAM";
+  void KMedoids::setSwapConfidence(size_t newSwapConfidence) {
+    if (algorithm != "BanditPAM" && algorithm != "BanditPAM_orig") {
+        // TODO(@motiwari): Better error type
+        throw "Cannot set buildConfidence when not using BanditPAM";
+    }
+    swapConfidence = newSwapConfidence;
   }
-  swapConfidence = newSwapConfidence;
-}
 
-void KMedoids::setSeed(size_t newSeed) {
-  seed = newSeed;
-  arma::arma_rng::set_seed(seed);
-}
+  void KMedoids::setSeed(size_t newSeed) {
+    seed = newSeed;
+    arma::arma_rng::set_seed(seed);
+  }
 
-size_t KMedoids::getSeed() const {
-  return seed;
-}
+  size_t KMedoids::getSeed() const {
+    return seed;
+  }
 
-bool KMedoids::getUseCache() const {
-  return useCache;
-}
+  bool KMedoids::getUseCache() const {
+    return useCache;
+  }
 
-void KMedoids::setUseCache(bool newUseCache) {
-  // TODO(@motiwari): Throw an error if not using BanditPAM
-  useCache = newUseCache;
-}
+  void KMedoids::setUseCache(bool newUseCache) {
+    // TODO(@motiwari): Throw an error if not using BanditPAM
+    useCache = newUseCache;
+  }
 
-bool KMedoids::getUsePerm() const {
+  bool KMedoids::getUsePerm() const {
     return usePerm;
-}
+  }
 
-void KMedoids::setUsePerm(bool newUsePerm) {
-  // TODO(@motiwari): Throw an error if not using BanditPAM
-  usePerm = newUsePerm;
-}
+  void KMedoids::setUsePerm(bool newUsePerm) {
+    // TODO(@motiwari): Throw an error if not using BanditPAM
+    usePerm = newUsePerm;
+  }
 
-size_t KMedoids::getCacheWidth() const {
-  // TODO(@motiwari): Throw an error if not using BanditPAM
-  return cacheWidth;
-}
+  size_t KMedoids::getCacheWidth() const {
+    // TODO(@motiwari): Throw an error if not using BanditPAM
+    return cacheWidth;
+  }
 
-void KMedoids::setCacheWidth(size_t newCacheWidth) {
-  // TODO(@motiwari): Throw an error if not using BanditPAM
-  cacheWidth = newCacheWidth;
-}
+  void KMedoids::setCacheWidth(size_t newCacheWidth) {
+    // TODO(@motiwari): Throw an error if not using BanditPAM
+    cacheWidth = newCacheWidth;
+  }
 
-bool KMedoids::getParallelize() const {
-  return parallelize;
-}
+  bool KMedoids::getParallelize() const {
+    return parallelize;
+  }
 
 // TODO(@motiwari): Change this to const bool newParallelize
-void KMedoids::setParallelize(bool newParallelize) {
-  parallelize = newParallelize;
-}
+  void KMedoids::setParallelize(bool newParallelize) {
+    parallelize = newParallelize;
+  }
 
-size_t KMedoids::getDistanceComputations(const bool includeMisc) const {
+  size_t KMedoids::getDistanceComputations(const bool includeMisc) const {
     if (includeMisc) {
-        return numMiscDistanceComputations +
-          numBuildDistanceComputations +
-          numSwapDistanceComputations;
+      return numMiscDistanceComputations +
+             numBuildDistanceComputations +
+             numSwapDistanceComputations;
     } else {
-        return numBuildDistanceComputations +
-          numSwapDistanceComputations;
+      return numBuildDistanceComputations +
+             numSwapDistanceComputations;
     }
-}
+  }
 
-size_t KMedoids::getMiscDistanceComputations() const {
+  size_t KMedoids::getMiscDistanceComputations() const {
     return numMiscDistanceComputations;
-}
-
-size_t KMedoids::getBuildDistanceComputations() const {
-  return numBuildDistanceComputations;
-}
-
-size_t KMedoids::getSwapDistanceComputations() const {
-  return numSwapDistanceComputations;
-}
-
-size_t KMedoids::getCacheWrites() const {
-  return numCacheWrites;
-}
-
-size_t KMedoids::getCacheHits() const {
-  return numCacheHits;
-}
-
-size_t KMedoids::getCacheMisses() const {
-  return numCacheMisses;
-}
-
-size_t KMedoids::getTotalSwapTime() const {
-  return totalSwapTime;
-}
-
-float KMedoids::getTimePerSwap() const {
-  return totalSwapTime / steps;
-}
-
-void KMedoids::setLossFn(std::string loss) {
-  // TODO(@motiwari): On setting this, clear the cache and the average loss,
-  // assignments, medoids, etc.
-  std::for_each(loss.begin(), loss.end(), [](char& c){
-    c = ::tolower(c);  // TODO(@motiwari): Put something before ::
-  });
-  // TODO(@motiwari): Change this to a switch
-  if (std::regex_match(loss, std::regex("l\\d*"))) {
-    lossFn = &KMedoids::LP;
-    lp = stoi(loss.substr(1));
-  } else if (loss == "manhattan") {
-    lossFn = &KMedoids::manhattan;
-  } else if (loss == "cos" || loss == "cosine") {
-    lossFn = &KMedoids::cos;
-  } else if (loss == "inf") {
-    lossFn = &KMedoids::LINF;
-  } else if (loss == "euclidean") {
-    lossFn = &KMedoids::LP;
-    lp = 2;
-  } else {
-    throw std::invalid_argument("Error: unrecognized loss function");
-  }
-}
-
-std::string KMedoids::getLossFn() const {
-  // TODO(@motiwari): make the strings constants
-  if (lossFn == &KMedoids::manhattan) {
+  }
+
+  size_t KMedoids::getBuildDistanceComputations() const {
+    return numBuildDistanceComputations;
+  }
+
+  size_t KMedoids::getSwapDistanceComputations() const {
+    return numSwapDistanceComputations;
+  }
+
+  size_t KMedoids::getCacheWrites() const {
+    return numCacheWrites;
+  }
+
+  size_t KMedoids::getCacheHits() const {
+    return numCacheHits;
+  }
+
+  size_t KMedoids::getCacheMisses() const {
+    return numCacheMisses;
+  }
+
+  size_t KMedoids::getTotalSwapTime() const {
+    return totalSwapTime;
+  }
+
+  float KMedoids::getTimePerSwap() const {
+    return totalSwapTime / steps;
+  }
+
+  void KMedoids::setLossFn(std::string loss) {
+    // TODO(@motiwari): On setting this, clear the
+    //  cache and the average loss,
+    // assignments, medoids, etc.
+    std::for_each(loss.begin(), loss.end(), [](char &c) {
+      c = ::tolower(c);  // TODO(@motiwari): Put something before ::
+    });
+    // TODO(@motiwari): Change this to a switch
+    if (std::regex_match(loss, std::regex("l\\d*"))) {
+      lossFn = &KMedoids::LP;
+      lp = stoi(loss.substr(1));
+    } else if (loss == "manhattan") {
+      lossFn = &KMedoids::manhattan;
+    } else if (loss == "cos" || loss == "cosine") {
+      lossFn = &KMedoids::cos;
+    } else if (loss == "inf") {
+      lossFn = &KMedoids::LINF;
+    } else if (loss == "euclidean") {
+      lossFn = &KMedoids::LP;
+      lp = 2;
+    } else {
+      throw std::invalid_argument("Error: unrecognized loss function");
+    }
+  }
+
+  std::string KMedoids::getLossFn() const {
+    // TODO(@motiwari): make the strings constants
+    if (lossFn == &KMedoids::manhattan) {
       return "manhattan";
-  } else if (lossFn == &KMedoids::cos) {
-    return "cosine";
-  } else if (lossFn == &KMedoids::LINF) {
-    return "L-infinity";
-  } else if (lossFn == &KMedoids::LP) {
-    return "L" + std::to_string(lp);
-  } else {
-    throw std::invalid_argument("Error: Loss Function Undefined!");
-  }
-}
-
-
-void KMedoids::calcBestDistancesSwap(
-  const arma::fmat& data,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-  const arma::urowvec* medoidIndices,
-  arma::frowvec* bestDistances,
-  arma::frowvec* secondBestDistances,
-  arma::urowvec* assignments,
-  const bool swapPerformed) {
-  #pragma omp parallel for if (this->parallelize)
-  for (size_t i = 0; i < data.n_cols; i++) {
-    float best = std::numeric_limits<float>::infinity();
-    float second = std::numeric_limits<float>::infinity();
-    for (size_t k = 0; k < medoidIndices->n_cols; k++) {
-      // 0 for MISC
-      float cost =
-          KMedoids::cachedLoss(data, distMat, i, (*medoidIndices)(k), 0);
-      if (cost < best) {
-        (*assignments)(i) = k;
-        second = best;
-        best = cost;
-      } else if (cost < second) {
-        second = cost;
-      }
+    } else if (lossFn == &KMedoids::cos) {
+      return "cosine";
+    } else if (lossFn == &KMedoids::LINF) {
+      return "L-infinity";
+    } else if (lossFn == &KMedoids::LP) {
+      return "L" + std::to_string(lp);
+    } else {
+      throw std::invalid_argument("Error: Loss Function Undefined!");
     }
-    (*bestDistances)(i) = best;
-    (*secondBestDistances)(i) = second;
-  }
-
-  if (!swapPerformed) {
-    // We have converged; update the final loss
-    averageLoss = arma::accu(*bestDistances) / data.n_cols;
-  }
-}
-
-float KMedoids::calcLoss(
-  const arma::fmat& data,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-  const arma::urowvec* medoidIndices) {
-  float total = 0;
-  // TODO(@motiwari): is this parallel loop accumulating properly?
-  #pragma omp parallel for if (this->parallelize)
-  for (size_t i = 0; i < data.n_cols; i++) {
-    float cost = std::numeric_limits<float>::infinity();
-    for (size_t k = 0; k < nMedoids; k++) {
-      float currCost = KMedoids::cachedLoss(
-        data,
-        distMat,
-        i,
-        (*medoidIndices)(k),
-        0);  // 0 for Misc
-      if (currCost < cost) {
-        cost = currCost;
+  }
+
+
+  void KMedoids::calcBestDistancesSwap(
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const arma::urowvec *medoidIndices,
+          arma::frowvec *bestDistances,
+          arma::frowvec *secondBestDistances,
+          arma::urowvec *assignments,
+          const bool swapPerformed) {
+    #pragma omp parallel for if (this->parallelize)
+    for (size_t i = 0; i < data.n_cols; i++) {
+      float best = std::numeric_limits<float>::infinity();
+      float second = std::numeric_limits<float>::infinity();
+      for (size_t k = 0; k < medoidIndices->n_cols; k++) {
+        // 0 for MISC
+        float cost =
+                KMedoids::cachedLoss(data, distMat, i,
+                                     (*medoidIndices)(k), 0);
+        if (cost < best) {
+          (*assignments)(i) = k;
+          second = best;
+          best = cost;
+        } else if (cost < second) {
+          second = cost;
+        }
       }
+      (*bestDistances)(i) = best;
+      (*secondBestDistances)(i) = second;
+    }
+
+    if (!swapPerformed) {
+      // We have converged; update the final loss
+      averageLoss = arma::accu(*bestDistances) / data.n_cols;
     }
-    total += cost;
   }
 
-  // Returns average distance
-  return total/data.n_cols;
-}
+  float KMedoids::calcLoss(
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const arma::urowvec *medoidIndices) {
+    float total = 0;
+    // TODO(@motiwari): is this parallel loop accumulating properly?
+    #pragma omp parallel for if (this->parallelize)
+    for (size_t i = 0; i < data.n_cols; i++) {
+      float cost = std::numeric_limits<float>::infinity();
+      for (size_t k = 0; k < nMedoids; k++) {
+        float currCost = KMedoids::cachedLoss(
+                data,
+                distMat,
+                i,
+                (*medoidIndices)(k),
+                0);  // 0 for Misc
+        if (currCost < cost) {
+          cost = currCost;
+        }
+      }
+      total += cost;
+    }
+
+    // Returns average distance
+    return total / data.n_cols;
+  }
 
-float KMedoids::cachedLoss(
-  const arma::fmat& data,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-  const size_t i,
-  const size_t j,
-  const size_t category,
-  const bool useCacheFunctionOverride
+  float KMedoids::cachedLoss(
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          const size_t i,
+          const size_t j,
+          const size_t category,
+          const bool useCacheFunctionOverride
   ) {
     // TODO(@motiwari): Change category to an enum
     if (category == 0) {  // MISC
-        numMiscDistanceComputations++;
+      numMiscDistanceComputations++;
     } else if (category == 1) {  // BUILD
-        numBuildDistanceComputations++;
+      numBuildDistanceComputations++;
     } else if (category == 2) {  // SWAP
-        numSwapDistanceComputations++;
+      numSwapDistanceComputations++;
     } else {
-        // TODO(@motiwari): Throw exception
+      // TODO(@motiwari): Throw exception
     }
 
-  if (this->useDistMat) {
-    return distMat.value().get().at(i, j);
-  }
+    if (this->useDistMat) {
+      return distMat.value().get().at(i, j);
+    }
+
+    if (!useCache) {
+      return (this->*lossFn)(data, i, j);
+    }
 
-  if (!useCache) {
+    // TODO(@motiwari): Should infer n and m from the size of the cache
+    size_t n = data.n_cols;
+    size_t m = fmin(n, cacheWidth);
+
+    // test this is one of the early points in the permutation
+    if (reindex.find(j) != reindex.end()) {
+      // TODO(@motiwari): Potential race condition with shearing?
+      // T1 begins to write to cache and then T2
+      // access in the middle of write?
+      if (cache[(m * i) + reindex[j]] == -1) {
+        numCacheWrites++;
+        cache[(m * i) + reindex[j]] = (this->*lossFn)(data, i, j);
+      }
+      numCacheHits++;
+      return cache[m * i + reindex[j]];
+    }
+
+    numCacheMisses++;
     return (this->*lossFn)(data, i, j);
   }
 
-  // TODO(@motiwari): Should infer n and m from the size of the cache
-  size_t n = data.n_cols;
-  size_t m = fmin(n, cacheWidth);
-
-  // test this is one of the early points in the permutation
-  if (reindex.find(j) != reindex.end()) {
-    // TODO(@motiwari): Potential race condition with shearing?
-    // T1 begins to write to cache and then T2 access in the middle of write?
-    if (cache[(m*i) + reindex[j]] == -1) {
-      numCacheWrites++;
-      cache[(m*i) + reindex[j]] = (this->*lossFn)(data, i, j);
+  void KMedoids::checkAlgorithm(const std::string &algorithm) const {
+    if ((algorithm != "BanditPAM") &&
+        (algorithm != "BanditPAM_orig") &&
+        (algorithm != "PAM") &&
+        (algorithm != "FastPAM1")) {
+        // TODO(@motiwari): Better error type
+      throw "unrecognized algorithm";
     }
-    numCacheHits++;
-    return cache[m*i + reindex[j]];
-  }
-
-  numCacheMisses++;
-  return (this->*lossFn)(data, i, j);
-}
-
-void KMedoids::checkAlgorithm(const std::string& algorithm) const {
-  if ((algorithm != "BanditPAM") &&
-      (algorithm != "BanditPAM_orig") &&
-      (algorithm != "PAM") &&
-      (algorithm != "FastPAM1")) {
-    // TODO(@motiwari): Better error type
-    throw "unrecognized algorithm";
-  }
-}
-
-float KMedoids::getAverageLoss() const {
-  return averageLoss;
-}
-
-float KMedoids::LP(const arma::fmat& data,
-  const size_t i,
-  const size_t j) const {
-  return arma::norm(data.col(i) - data.col(j), lp);
-}
-
-float KMedoids::LINF(
-  const arma::fmat& data,
-  const size_t i,
-  const size_t j) const {
-  return arma::max(arma::abs(data.col(i) - data.col(j)));
-}
-
-float KMedoids::cos(const arma::fmat& data,
-  const size_t i,
-  const size_t j) const {
-  return 1 - (arma::dot(data.col(i), data.col(j))
-    / (arma::norm(data.col(i)) * arma::norm(data.col(j))));
-}
-
-float KMedoids::manhattan(const arma::fmat& data,
-  const size_t i,
-  const size_t j) const {
-  return arma::accu(arma::abs(data.col(i) - data.col(j)));
-}
+  }
+
+  float KMedoids::getAverageLoss() const {
+    return averageLoss;
+  }
+
+  float KMedoids::getBuildLoss() const {
+    return buildLoss;
+  }
+
+  float KMedoids::LP(const arma::fmat &data,
+                     const size_t i,
+                     const size_t j) const {
+    return arma::norm(data.col(i) - data.col(j), lp);
+  }
+
+  float KMedoids::LINF(
+          const arma::fmat &data,
+          const size_t i,
+          const size_t j) const {
+    return arma::max(arma::abs(data.col(i) - data.col(j)));
+  }
+
+  float KMedoids::cos(const arma::fmat &data,
+                      const size_t i,
+                      const size_t j) const {
+    return 1 - (arma::dot(data.col(i), data.col(j))
+                / (arma::norm(data.col(i)) * arma::norm(data.col(j))));
+  }
+
+  float KMedoids::manhattan(const arma::fmat &data,
+                            const size_t i,
+                            const size_t j) const {
+    return arma::accu(arma::abs(data.col(i) - data.col(j)));
+  }
 }  // namespace km
diff --git a/src/algorithms/pam.cpp b/src/algorithms/pam.cpp
index f6503655..7fe62bd6 100644
--- a/src/algorithms/pam.cpp
+++ b/src/algorithms/pam.cpp
@@ -14,117 +14,121 @@
 #include <unordered_map>
 
 namespace km {
-void PAM::fitPAM(
-  const arma::fmat& inputData,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat) {
-  data = arma::trans(inputData);
-  arma::urowvec medoidIndices(nMedoids);
-  PAM::buildPAM(data, distMat, &medoidIndices);
-  steps = 0;
-  medoidIndicesBuild = medoidIndices;
-  arma::urowvec assignments(data.n_cols);
+  void PAM::fitPAM(
+          const arma::fmat &inputData,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat) {
+    data = arma::trans(inputData);
+    arma::urowvec medoidIndices(nMedoids);
+    PAM::buildPAM(data, distMat, &medoidIndices);
+    steps = 0;
+    medoidIndicesBuild = medoidIndices;
+    arma::urowvec assignments(data.n_cols);
 
-  size_t i = 0;
-  bool medoidChange = true;
-  while (i < maxIter && medoidChange) {
-    auto previous(medoidIndices);
-    PAM::swapPAM(data, distMat, &medoidIndices, &assignments);
-    medoidChange = arma::any(medoidIndices != previous);
-    i++;
+    size_t i = 0;
+    bool medoidChange = true;
+    while (i < maxIter && medoidChange) {
+      auto previous(medoidIndices);
+      if (nMedoids > 1) {
+        PAM::swapPAM(data, distMat, &medoidIndices, &assignments);
+      }
+
+      medoidChange = arma::any(medoidIndices != previous);
+      i++;
+    }
+    medoidIndicesFinal = medoidIndices;
+    labels = assignments;
+    steps = i;
   }
-  medoidIndicesFinal = medoidIndices;
-  labels = assignments;
-  steps = i;
-}
 
-void PAM::buildPAM(
-  const arma::fmat& data,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-  arma::urowvec* medoidIndices) {
-  size_t N = data.n_cols;
-  arma::frowvec estimates(N, arma::fill::zeros);
-  arma::frowvec bestDistances(N);
-  bestDistances.fill(std::numeric_limits<float>::infinity());
-  for (size_t k = 0; k < nMedoids; k++) {
-    float minDistance = std::numeric_limits<float>::infinity();
-    size_t best = 0;
-    for (size_t i = 0; i < data.n_cols; i++) {
-      float total = 0;
-      for (size_t j = 0; j < data.n_cols; j++) {
-        float cost = (this->*lossFn)(data, i, j);
-        // compares this with the cached best distance
-        if (bestDistances(j) < cost) {
-          cost = bestDistances(j);
+  void PAM::buildPAM(
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          arma::urowvec *medoidIndices) {
+    size_t N = data.n_cols;
+    arma::frowvec estimates(N, arma::fill::zeros);
+    arma::frowvec bestDistances(N);
+    bestDistances.fill(std::numeric_limits<float>::infinity());
+    for (size_t k = 0; k < nMedoids; k++) {
+      float minDistance = std::numeric_limits<float>::infinity();
+      size_t best = 0;
+      for (size_t i = 0; i < data.n_cols; i++) {
+        float total = 0;
+        for (size_t j = 0; j < data.n_cols; j++) {
+          float cost = (this->*lossFn)(data, i, j);
+          // compares this with the cached best distance
+          if (bestDistances(j) < cost) {
+              cost = bestDistances(j);
+          }
+          total += cost;
+        }
+        if (total < minDistance) {
+          minDistance = total;
+          best = i;
         }
-        total += cost;
-      }
-      if (total < minDistance) {
-        minDistance = total;
-        best = i;
       }
-    }
-    (*medoidIndices)(k) = best;
+      (*medoidIndices)(k) = best;
 
-    // update the medoid assignment and best_distance for this datapoint
-    for (size_t l = 0; l < N; l++) {
-      float cost = (this->*lossFn)(data, l, (*medoidIndices)(k));
-      if (cost < bestDistances(l)) {
-        bestDistances(l) = cost;
+      // update the medoid assignment and best_distance for this datapoint
+      for (size_t l = 0; l < N; l++) {
+        float cost = (this->*lossFn)(data, l, (*medoidIndices)(k));
+        if (cost < bestDistances(l)) {
+          bestDistances(l) = cost;
+        }
       }
     }
   }
-}
 
-void PAM::swapPAM(
-  const arma::fmat& data,
-  std::optional<std::reference_wrapper<const arma::fmat>> distMat,
-  arma::urowvec* medoidIndices,
-  arma::urowvec* assignments) {
-  float minDistance = std::numeric_limits<float>::infinity();
-  size_t best = 0;
-  size_t medoidToSwap = 0;
-  size_t N = data.n_cols;
-  arma::frowvec bestDistances(N);
-  arma::frowvec secondBestDistances(N);
+  void PAM::swapPAM(
+          const arma::fmat &data,
+          std::optional<std::reference_wrapper<const arma::fmat>> distMat,
+          arma::urowvec *medoidIndices,
+          arma::urowvec *assignments) {
+    float minDistance = std::numeric_limits<float>::infinity();
+    size_t best = 0;
+    size_t medoidToSwap = 0;
+    size_t N = data.n_cols;
+    arma::frowvec bestDistances(N);
+    arma::frowvec secondBestDistances(N);
 
-  KMedoids::calcBestDistancesSwap(
-    data,
-    distMat,
-    medoidIndices,
-    &bestDistances,
-    &secondBestDistances,
-    assignments);
+    KMedoids::calcBestDistancesSwap(
+            data,
+            distMat,
+            medoidIndices,
+            &bestDistances,
+            &secondBestDistances,
+            assignments);
 
-  for (size_t k = 0; k < nMedoids; k++) {
-    for (size_t i = 0; i < data.n_cols; i++) {
-      float total = 0;
-      for (size_t j = 0; j < data.n_cols; j++) {
-        // compute distance between base point and every other datapoint
-        float cost = (this->*lossFn)(data, i, j);
-        // if x_j is NOT assigned to k: compares this with
-        //   the cached best distance
-        // if x_j is assigned to k: compares this with
-        //   the cached second best distance
-        if ((*assignments)(j) != k) {
-          if (bestDistances(j) < cost) {
-            cost = bestDistances(j);
-          }
-        } else {
-          if (secondBestDistances(j) < cost) {
-            cost = secondBestDistances(j);
+    for (size_t k = 0; k < nMedoids; k++) {
+      for (size_t i = 0; i < data.n_cols; i++) {
+        float total = 0;
+        for (size_t j = 0; j < data.n_cols; j++) {
+          // compute distance between base point
+          // and every other datapoint
+          float cost = (this->*lossFn)(data, i, j);
+          // if x_j is NOT assigned to k: compares this with
+          //   the cached best distance
+          // if x_j is assigned to k: compares this with
+          //   the cached second best distance
+          if ((*assignments)(j) != k) {
+            if (bestDistances(j) < cost) {
+              cost = bestDistances(j);
+            }
+          } else {
+            if (secondBestDistances(j) < cost) {
+              cost = secondBestDistances(j);
+            }
           }
+          total += cost;
+        }
+        // if total distance for new base point is better than
+        // that of the medoid, update the best index identified so far
+        if (total < minDistance) {
+          minDistance = total;
+          best = i;
+          medoidToSwap = k;
         }
-        total += cost;
-      }
-      // if total distance for new base point is better than
-      // that of the medoid, update the best index identified so far
-      if (total < minDistance) {
-        minDistance = total;
-        best = i;
-        medoidToSwap = k;
       }
     }
+    (*medoidIndices)(medoidToSwap) = best;
   }
-  (*medoidIndices)(medoidToSwap) = best;
-}
 }  // namespace km
diff --git a/src/main.cpp b/src/main.cpp
index c33ef9dd..0a65609d 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -16,114 +16,114 @@
 
 #include "kmedoids_algorithm.hpp"
 
-int main(int argc, char* argv[]) {
-    std::string input_name;
-    size_t k;
-    int opt;
-    int prev_ind;
-    size_t maxIter = 100;
-    size_t buildConfidence = 3;
-    size_t swapConfidence = 5;
-    std::string loss = "L2";
-    bool f_flag = false;
-    bool k_flag = false;
-    const int ARGUMENT_ERROR_CODE = 1;
-    bool useCache = true;
-    bool usePerm = true;
-    int seed = 0;
-    int num_data = 0;
-    bool parallelize = false;
+int main(int argc, char *argv[]) {
+  std::string input_name;
+  size_t k;
+  int opt;
+  int prev_ind;
+  size_t maxIter = 100;
+  size_t buildConfidence = 3;
+  size_t swapConfidence = 5;
+  std::string loss = "L2";
+  bool f_flag = false;
+  bool k_flag = false;
+  const int ARGUMENT_ERROR_CODE = 1;
+  bool useCache = true;
+  bool usePerm = true;
+  int seed = 0;
+  int num_data = 0;
+  bool parallelize = false;
 
-    // TODO(@motiwari): Use a variadic function signature for this instead
-    while (prev_ind = optind,
-        (opt = getopt(argc, argv, "f:l:k:v:s:cpnw:")) != -1) {
-        if ( optind == prev_ind + 2 && *optarg == '-' ) {
-        opt = ':';
-        --optind;
-        }
+  // TODO(@motiwari): Use a variadic function signature for this instead
+  while (prev_ind = optind,
+          (opt = getopt(argc, argv, "f:l:k:v:s:cpnw:")) != -1) {
+    if (optind == prev_ind + 2 && *optarg == '-') {
+      opt = ':';
+      --optind;
+    }
 
-        switch (opt) {
-            // path to the data file to be read in
-            case 'f':
-                input_name = optarg;
-                f_flag = true;
-                break;
-            // number of clusters to create
-            case 'k':
-                k = std::stoi(optarg);
-                k_flag = true;
-                break;
-            // type of loss/distance function to use
-            case 'l':
-                loss = optarg;
-                break;
-            case 's':
-                seed = std::stoi(optarg);
-                break;
-            case ':':
-                printf("option needs a value\n");
-                return ARGUMENT_ERROR_CODE;
-            case 'c':
-                useCache = true;
-                break;
-            case 'p':
-                usePerm = true;
-                break;
-            case 'n':
-                num_data = std::stoi(optarg);
-                break;
-            case 'w':
-                parallelize = true;
-                break;
-            case '?':
-                printf("unknown option: %c\n", optopt);
-                return ARGUMENT_ERROR_CODE;
-        }
+    switch (opt) {
+      // path to the data file to be read in
+      case 'f':
+        input_name = optarg;
+        f_flag = true;
+        break;
+        // number of clusters to create
+      case 'k':
+        k = std::stoi(optarg);
+        k_flag = true;
+        break;
+        // type of loss/distance function to use
+      case 'l':
+        loss = optarg;
+        break;
+      case 's':
+        seed = std::stoi(optarg);
+        break;
+      case ':':
+        printf("option needs a value\n");
+        return ARGUMENT_ERROR_CODE;
+      case 'c':
+        useCache = true;
+        break;
+      case 'p':
+        usePerm = true;
+        break;
+      case 'n':
+        num_data = std::stoi(optarg);
+        break;
+      case 'w':
+        parallelize = true;
+        break;
+      case '?':
+        printf("unknown option: %c\n", optopt);
+        return ARGUMENT_ERROR_CODE;
     }
+  }
 
-    try {
-      if (!f_flag) {
-        throw std::invalid_argument(
-          "Error: Must specify input file via -f flag");
-      } else if (!k_flag) {
-        throw std::invalid_argument(
-          "Error: Must specify number of clusters via -k flag");
-      } else if (!std::filesystem::exists(input_name)) {
-        throw std::invalid_argument(
-          "Error: The file does not exist");
-      }
-    } catch (std::invalid_argument& e) {
-      std::cout << e.what() << std::endl;
-      return ARGUMENT_ERROR_CODE;
+  try {
+    if (!f_flag) {
+      throw std::invalid_argument(
+              "Error: Must specify input file via -f flag");
+    } else if (!k_flag) {
+      throw std::invalid_argument(
+              "Error: Must specify number of clusters via -k flag");
+    } else if (!std::filesystem::exists(input_name)) {
+      throw std::invalid_argument(
+              "Error: The file does not exist");
     }
+  } catch (std::invalid_argument &e) {
+    std::cout << e.what() << std::endl;
+    return ARGUMENT_ERROR_CODE;
+  }
 
-    arma::fmat data;
-    data.load(input_name);
+  arma::fmat data;
+  data.load(input_name);
 
-    if (num_data < 0) {
-        throw std::invalid_argument(
-        "Error: num_data passed was less than 0");
-    } else if (num_data != 0) {
-        data.resize(num_data, data.n_cols);
-    }
+  if (num_data < 0) {
+    throw std::invalid_argument(
+            "Error: num_data passed was less than 0");
+  } else if (num_data != 0) {
+    data.resize(num_data, data.n_cols);
+  }
 
-    km::KMedoids kmed(
-      k,
-      "BanditPAM",
-      maxIter,
-      buildConfidence,
-      swapConfidence,
-      useCache,
-      usePerm,
-      1000,  // Cache Width
-      parallelize,
-      seed);
-    kmed.fit(data, loss, {});
-    for (auto medoid : kmed.getMedoidsFinal()) {
-      std::cout << medoid << ",";
-    }
-    std::cout << "Final loss: " << kmed.getAverageLoss() << "\n";
-    std::cout << "Num Swap Steps: " << kmed.getSteps() << "\n";
-    std::cout << "Total Swap Milliseconds: " << kmed.getTotalSwapTime() << "\n";
-    std::cout << "Average Swap Milliseconds: " << kmed.getTimePerSwap() << "\n";
+  km::KMedoids kmed(
+          k,
+          "BanditPAM",
+          maxIter,
+          buildConfidence,
+          swapConfidence,
+          useCache,
+          usePerm,
+          1000,  // Cache Width
+          parallelize,
+          seed);
+  kmed.fit(data, loss, {});
+  for (auto medoid : kmed.getMedoidsFinal()) {
+    std::cout << medoid << ",";
+  }
+  std::cout << "Final loss: " << kmed.getAverageLoss() << "\n";
+  std::cout << "Num Swap Steps: " << kmed.getSteps() << "\n";
+  std::cout << "Total Swap Milliseconds: " << kmed.getTotalSwapTime() << "\n";
+  std::cout << "Average Swap Milliseconds: " << kmed.getTimePerSwap() << "\n";
 }
diff --git a/src/python_bindings/build_medoids_python.cpp b/src/python_bindings/build_medoids_python.cpp
index 2d123894..0453597e 100644
--- a/src/python_bindings/build_medoids_python.cpp
+++ b/src/python_bindings/build_medoids_python.cpp
@@ -14,17 +14,18 @@
 #include "kmedoids_pywrapper.hpp"
 
 namespace km {
-pybind11::array_t<arma::uword> km::KMedoidsWrapper::getMedoidsBuildPython() {
-  if (KMedoids::getMedoidsBuild().size() > 1) {
-    return
-      carma::row_to_arr<arma::uword>(KMedoids::getMedoidsBuild()).squeeze();
-  } else {
+  pybind11::array_t <arma::uword> km::KMedoidsWrapper::getMedoidsBuildPython() {
+    if (KMedoids::getMedoidsBuild().size() > 1) {
+      return
+              carma::row_to_arr<arma::uword>(
+                      KMedoids::getMedoidsBuild()).squeeze();
+    } else {
       return carma::row_to_arr<arma::uword>(KMedoids::getMedoidsBuild());
+    }
   }
-}
 
-void build_medoids_python(pybind11::class_<KMedoidsWrapper> *cls) {
-  cls->def_property_readonly("build_medoids",
-    &KMedoidsWrapper::getMedoidsBuildPython);
-}
+  void build_medoids_python(pybind11::class_ <KMedoidsWrapper> *cls) {
+    cls->def_property_readonly("build_medoids",
+                               &KMedoidsWrapper::getMedoidsBuildPython);
+  }
 }  // namespace km
diff --git a/src/python_bindings/cache_python.cpp b/src/python_bindings/cache_python.cpp
index 582aa0a5..7dd13af2 100644
--- a/src/python_bindings/cache_python.cpp
+++ b/src/python_bindings/cache_python.cpp
@@ -13,76 +13,77 @@
 #include "kmedoids_pywrapper.hpp"
 
 namespace km {
-size_t km::KMedoidsWrapper::getDistanceComputationsPython(
-  const bool includeMisc) {
-  return KMedoids::getDistanceComputations(includeMisc);
-}
+  size_t km::KMedoidsWrapper::getDistanceComputationsPython(
+          const bool includeMisc) {
+    return KMedoids::getDistanceComputations(includeMisc);
+  }
 
-size_t km::KMedoidsWrapper::getMiscDistanceComputationsPython() {
+  size_t km::KMedoidsWrapper::getMiscDistanceComputationsPython() {
     return KMedoids::getMiscDistanceComputations();
-}
+  }
 
-size_t km::KMedoidsWrapper::getBuildDistanceComputationsPython() {
+  size_t km::KMedoidsWrapper::getBuildDistanceComputationsPython() {
     return KMedoids::getBuildDistanceComputations();
-}
+  }
 
-size_t km::KMedoidsWrapper::getSwapDistanceComputationsPython() {
+  size_t km::KMedoidsWrapper::getSwapDistanceComputationsPython() {
     return KMedoids::getSwapDistanceComputations();
-}
+  }
 
-size_t km::KMedoidsWrapper::getCacheWritesPython() {
-  return KMedoids::getCacheWrites();
-}
+  size_t km::KMedoidsWrapper::getCacheWritesPython() {
+    return KMedoids::getCacheWrites();
+  }
 
-size_t km::KMedoidsWrapper::getCacheHitsPython() {
-  return KMedoids::getCacheHits();
-}
+  size_t km::KMedoidsWrapper::getCacheHitsPython() {
+    return KMedoids::getCacheHits();
+  }
 
-size_t km::KMedoidsWrapper::getCacheMissesPython() {
-  return KMedoids::getCacheMisses();
-}
+  size_t km::KMedoidsWrapper::getCacheMissesPython() {
+    return KMedoids::getCacheMisses();
+  }
 
-void distance_computations_python(pybind11::class_<km::KMedoidsWrapper> *cls) {
-  cls->def("getDistanceComputations",
-           &KMedoidsWrapper::getDistanceComputationsPython);
-}
+  void distance_computations_python(
+          pybind11::class_ <km::KMedoidsWrapper> *cls) {
+    cls->def("getDistanceComputations",
+             &KMedoidsWrapper::getDistanceComputationsPython);
+  }
 
-void misc_distance_computations_python(
-    pybind11::class_<km::KMedoidsWrapper> *cls) {
+  void misc_distance_computations_python(
+          pybind11::class_ <km::KMedoidsWrapper> *cls) {
     cls->def_property_readonly(
-      "misc_distance_computations",
-      &KMedoidsWrapper::getMiscDistanceComputationsPython);
-}
+            "misc_distance_computations",
+            &KMedoidsWrapper::getMiscDistanceComputationsPython);
+  }
 
-void build_distance_computations_python(
-    pybind11::class_<km::KMedoidsWrapper> *cls) {
+  void build_distance_computations_python(
+          pybind11::class_ <km::KMedoidsWrapper> *cls) {
     cls->def_property_readonly(
-      "build_distance_computations",
-      &KMedoidsWrapper::getBuildDistanceComputationsPython);
-}
+            "build_distance_computations",
+            &KMedoidsWrapper::getBuildDistanceComputationsPython);
+  }
 
-void swap_distance_computations_python(
-    pybind11::class_<km::KMedoidsWrapper> *cls) {
+  void swap_distance_computations_python(
+          pybind11::class_ <km::KMedoidsWrapper> *cls) {
     cls->def_property_readonly(
-      "swap_distance_computations",
-      &KMedoidsWrapper::getSwapDistanceComputationsPython);
-}
-
-void cache_writes_python(pybind11::class_<km::KMedoidsWrapper> *cls) {
-  cls->def_property_readonly(
-      "cache_writes",
-      &KMedoidsWrapper::getCacheWritesPython);
-}
-
-void cache_hits_python(pybind11::class_<km::KMedoidsWrapper> *cls) {
-  cls->def_property_readonly(
-      "cache_hits",
-      &KMedoidsWrapper::getCacheHitsPython);
-}
-
-void cache_misses_python(pybind11::class_<km::KMedoidsWrapper> *cls) {
-  cls->def_property_readonly(
-      "cache_misses",
-      &KMedoidsWrapper::getCacheMisses);
-}
+            "swap_distance_computations",
+            &KMedoidsWrapper::getSwapDistanceComputationsPython);
+  }
+
+  void cache_writes_python(pybind11::class_ <km::KMedoidsWrapper> *cls) {
+    cls->def_property_readonly(
+            "cache_writes",
+            &KMedoidsWrapper::getCacheWritesPython);
+  }
+
+  void cache_hits_python(pybind11::class_ <km::KMedoidsWrapper> *cls) {
+    cls->def_property_readonly(
+            "cache_hits",
+            &KMedoidsWrapper::getCacheHitsPython);
+  }
+
+  void cache_misses_python(pybind11::class_ <km::KMedoidsWrapper> *cls) {
+    cls->def_property_readonly(
+            "cache_misses",
+            &KMedoidsWrapper::getCacheMisses);
+  }
 }  // namespace km
diff --git a/src/python_bindings/fit_python.cpp b/src/python_bindings/fit_python.cpp
index e5666903..366374d8 100644
--- a/src/python_bindings/fit_python.cpp
+++ b/src/python_bindings/fit_python.cpp
@@ -16,44 +16,48 @@
 #include "kmedoids_pywrapper.hpp"
 
 namespace km {
-void km::KMedoidsWrapper::fitPython(
-  const pybind11::array_t<float>& inputData,
-  const std::string& loss,
-  pybind11::kwargs kw) {
-  // throw an error if the number of medoids is not specified in either
-  // the KMedoids object or the fitPython function
-  try {
-    if (KMedoids::getNMedoids() == 0) {  // Check for 0 as NULL
-      if (kw.size() == 0) {
-        throw pybind11::value_error("Error: must specify number of medoids.");
+  void km::KMedoidsWrapper::fitPython(
+          const pybind11::array_t<float> &inputData,
+          const std::string &loss,
+          pybind11::kwargs kw) {
+    // throw an error if the number of medoids is not specified in either
+    // the KMedoids object or the fitPython function
+    try {
+      if (KMedoids::getNMedoids() == 0) {  // Check for 0 as NULL
+        if (kw.size() == 0) {
+          throw pybind11::value_error(
+                  "Error: must specify number of medoids.");
+        }
       }
+    } catch (pybind11::value_error &e) {
+      // Throw it again (pybind11 will raise ValueError)
+      // TODO(@motiwari): Make this more informative
+      // TODO(@motiwari): Is this throw redundant with the above?
+      throw;
+    }
+    // if k is specified here, then
+    // we set the number of medoids as k and override previous value
+    if ((kw.size() != 0) && (kw.contains("k"))) {
+      KMedoids::setNMedoids(pybind11::cast<int>(kw["k"]));
     }
-  } catch (pybind11::value_error &e) {
-    // Throw it again (pybind11 will raise ValueError)
-    // TODO(@motiwari): Make this more informative
-    // TODO(@motiwari): Is this throw redundant with the above?
-    throw;
-  }
-  // if k is specified here, then
-  // we set the number of medoids as k and override previous value
-  if ((kw.size() != 0) && (kw.contains("k"))) {
-    KMedoids::setNMedoids(pybind11::cast<int>(kw["k"]));
-  }
 
-  // TODO(@motiwari): Add a comment here explaining the massaging, cleanup
-  if ((kw.size() != 0) && (kw.contains("dist_mat"))) {
-    const pybind11::array_t<float>& distMat =
-      pybind11::cast<const pybind11::array_t<float>>(kw["dist_mat"]);
-    const arma::fmat& tmp_array = carma::arr_to_mat<float>(distMat);
-    KMedoids::fit(carma::arr_to_mat<float>(inputData), loss,
-      std::make_optional<std::reference_wrapper<const arma::fmat>>(tmp_array));
-  } else {
-    // TODO(@motiwari): change std::nullopt to nullopt?
-    KMedoids::fit(carma::arr_to_mat<float>(inputData), loss, std::nullopt);
+    // TODO(@motiwari): Add a comment here explaining the massaging, cleanup
+    if ((kw.size() != 0) && (kw.contains("dist_mat"))) {
+      const pybind11::array_t<float> &distMat =
+              pybind11::cast<const pybind11::array_t<float>>(
+                      kw["dist_mat"]);
+      const arma::fmat &tmp_array = carma::arr_to_mat<float>(distMat);
+      KMedoids::fit(carma::arr_to_mat<float>(inputData), loss,
+            std::make_optional<std::reference_wrapper<const arma::fmat>>(
+                    tmp_array));
+    } else {
+      // TODO(@motiwari): change std::nullopt to nullopt?
+      KMedoids::fit(carma::arr_to_mat<float>(inputData), loss,
+                    std::nullopt);
+    }
   }
-}
 
-void fit_python(pybind11::class_<KMedoidsWrapper> *cls) {
-  cls->def("fit", &KMedoidsWrapper::fitPython);
-}
+  void fit_python(pybind11::class_ <KMedoidsWrapper> *cls) {
+    cls->def("fit", &KMedoidsWrapper::fitPython);
+  }
 }  // namespace km
diff --git a/src/python_bindings/kmedoids_pywrapper.cpp b/src/python_bindings/kmedoids_pywrapper.cpp
index 8b4dc660..c336909d 100644
--- a/src/python_bindings/kmedoids_pywrapper.cpp
+++ b/src/python_bindings/kmedoids_pywrapper.cpp
@@ -16,25 +16,26 @@
 #include "kmedoids_pywrapper.hpp"
 
 namespace km {
-PYBIND11_MODULE(banditpam, m) {
-  // Module functions
-  m.doc() = "BanditPAM Python library, implemented in C++";
-  m.attr("__version__") = VERSION_INFO;
-  m.def("get_max_threads",  // TODO(@motiwari): change to get_num_threads
+  PYBIND11_MODULE(banditpam, m) {
+    // Module functions
+    m.doc() = "BanditPAM Python library, implemented in C++";
+
+    m.attr("__version__") = VERSION_INFO;
+    m.def("get_max_threads",  // TODO(@motiwari): change to get_num_threads
     &omp_get_max_threads, "Returns max number of threads");
-  m.def("set_num_threads",
+    m.def("set_num_threads",
     &omp_set_num_threads, "Set the maximum number of threads");
 
-  // Class functions
-  pybind11::class_<KMedoidsWrapper> cls(m, "KMedoids");
+    // Class functions
+    pybind11::class_ <KMedoidsWrapper> cls(m, "KMedoids");
 
-  // Constructor
-  // NOTE: The order of these matters! Otherwise you can get warnings about
-  //  variables not being intialized in the same order as their declarations,
-  //  which can lead to undefined behavior (variables being intialized with
-  //  each others' values). The order here much also match that of the
-  //  constructor in kmedoids_algorithm.*pp
-  cls.def(
+    // Constructor
+    // NOTE: The order of these matters! Otherwise you can get warnings about
+    //  variables not being intialized in the same order as their declarations,
+    //  which can lead to undefined behavior (variables being intialized with
+    //  each others' values). The order here much also match that of the
+    //  constructor in kmedoids_algorithm.*pp
+    cls.def(
     pybind11::init<int, std::string, int, int, int, bool, bool, int, bool>(),
     pybind11::arg("n_medoids") = 5,
     pybind11::arg("algorithm") = "BanditPAM",
@@ -47,49 +48,50 @@ PYBIND11_MODULE(banditpam, m) {
     pybind11::arg("cache_width") = 1000,
     pybind11::arg("parallelize") = true);
 
-  // Properties
-  cls.def_property("n_medoids",
+    // Properties
+    cls.def_property("n_medoids",
     &KMedoidsWrapper::getNMedoids, &KMedoidsWrapper::setNMedoids);
-  cls.def_property("algorithm",
+    cls.def_property("algorithm",
     &KMedoidsWrapper::getAlgorithm, &KMedoidsWrapper::setAlgorithm);
-  cls.def_property("max_iter",
+    cls.def_property("max_iter",
     &KMedoidsWrapper::getMaxIter, &KMedoidsWrapper::setMaxIter);
-  cls.def_property("build_confidence",
+    cls.def_property("build_confidence",
     &KMedoidsWrapper::getBuildConfidence, &KMedoidsWrapper::setBuildConfidence);
-  cls.def_property("swap_confidence",
+    cls.def_property("swap_confidence",
     &KMedoidsWrapper::getSwapConfidence, &KMedoidsWrapper::setSwapConfidence);
-  cls.def_property("use_cache",
+    cls.def_property("use_cache",
     &KMedoidsWrapper::getUseCache, &KMedoidsWrapper::setUseCache);
-  cls.def_property("use_perm",
+    cls.def_property("use_perm",
     &KMedoidsWrapper::getUsePerm, &KMedoidsWrapper::setUsePerm);
-  cls.def_property("cache_width",
+    cls.def_property("cache_width",
     &KMedoidsWrapper::getCacheWidth, &KMedoidsWrapper::setCacheWidth);
-  cls.def_property("parallelize",
+    cls.def_property("parallelize",
     &KMedoidsWrapper::getParallelize, &KMedoidsWrapper::setParallelize);
-  cls.def_property("loss_function",
+    cls.def_property("loss_function",
     &KMedoidsWrapper::getLossFn, &KMedoidsWrapper::setLossFn);
-  cls.def_property("seed",
+    cls.def_property("seed",
     &KMedoidsWrapper::getSeed, &KMedoidsWrapper::setSeed);
 
-  // Other functions
-  medoids_python(&cls);
-  build_medoids_python(&cls);
-  labels_python(&cls);
-  steps_python(&cls);
-  fit_python(&cls);
-  loss_python(&cls);
+    // Other functions
+    medoids_python(&cls);
+    build_medoids_python(&cls);
+    labels_python(&cls);
+    steps_python(&cls);
+    fit_python(&cls);
+    loss_python(&cls);
+    build_loss_python(&cls);
 
-  // Cache functions
-  distance_computations_python(&cls);
-  misc_distance_computations_python(&cls);
-  build_distance_computations_python(&cls);
-  swap_distance_computations_python(&cls);
-  cache_writes_python(&cls);
-  cache_hits_python(&cls);
-  cache_misses_python(&cls);
+    // Cache functions
+    distance_computations_python(&cls);
+    misc_distance_computations_python(&cls);
+    build_distance_computations_python(&cls);
+    swap_distance_computations_python(&cls);
+    cache_writes_python(&cls);
+    cache_hits_python(&cls);
+    cache_misses_python(&cls);
 
-  // Swap timing functions
-  time_per_swap_python(&cls);
-  total_swap_time_python(&cls);
-}
+    // Swap timing functions
+    time_per_swap_python(&cls);
+    total_swap_time_python(&cls);
+  }
 }  // namespace km
diff --git a/src/python_bindings/labels_python.cpp b/src/python_bindings/labels_python.cpp
index ff21d7a6..7b465883 100644
--- a/src/python_bindings/labels_python.cpp
+++ b/src/python_bindings/labels_python.cpp
@@ -14,15 +14,17 @@
 #include "kmedoids_pywrapper.hpp"
 
 namespace km {
-pybind11::array_t<arma::uword> km::KMedoidsWrapper::getLabelsPython() {
-  if (KMedoids::getLabels().size() > 1) {
-    return carma::row_to_arr<arma::uword>(KMedoids::getLabels()).squeeze();
-  } else {
-    return carma::row_to_arr<arma::uword>(KMedoids::getLabels());
+  pybind11::array_t <arma::uword> km::KMedoidsWrapper::getLabelsPython() {
+    if (KMedoids::getLabels().size() > 1) {
+      return carma::row_to_arr<arma::uword>(
+                KMedoids::getLabels()).squeeze();
+    } else {
+      return carma::row_to_arr<arma::uword>(KMedoids::getLabels());
+    }
   }
-}
 
-void labels_python(pybind11::class_<km::KMedoidsWrapper> *cls) {
-  cls->def_property_readonly("labels", &km::KMedoidsWrapper::getLabelsPython);
-}
+  void labels_python(pybind11::class_ <km::KMedoidsWrapper> *cls) {
+    cls->def_property_readonly("labels",
+                               &km::KMedoidsWrapper::getLabelsPython);
+  }
 }  // namespace km
diff --git a/src/python_bindings/loss_python.cpp b/src/python_bindings/loss_python.cpp
index 4ec64e43..0454be69 100644
--- a/src/python_bindings/loss_python.cpp
+++ b/src/python_bindings/loss_python.cpp
@@ -2,7 +2,7 @@
  * @file loss_python.cpp
  * @date 2021-08-16
  *
- * Defines the function getLossPython in KMedoidsWrapper class 
+ * Defines the function getLossPython and getBuildLossPython in KMedoidsWrapper class
  * which is used in Python bindings.
  */
 
@@ -14,11 +14,21 @@
 #include "kmedoids_pywrapper.hpp"
 
 namespace km {
-float km::KMedoidsWrapper::getLossPython() {
-  return KMedoids::getAverageLoss();
-}
+  float km::KMedoidsWrapper::getLossPython() {
+    return KMedoids::getAverageLoss();
+  }
 
-void loss_python(pybind11::class_<KMedoidsWrapper> *cls) {
-  cls->def_property_readonly("average_loss", &KMedoidsWrapper::getLossPython);
-}
+  float km::KMedoidsWrapper::getBuildLossPython() {
+    return KMedoids::getBuildLoss();
+  }
+
+  void loss_python(pybind11::class_ <KMedoidsWrapper> *cls) {
+    cls->def_property_readonly("average_loss",
+                               &KMedoidsWrapper::getLossPython);
+  }
+
+  void build_loss_python(pybind11::class_ <KMedoidsWrapper> *cls) {
+    cls->def_property_readonly("build_loss",
+                               &KMedoidsWrapper::getBuildLossPython);
+  }
 }  // namespace km
diff --git a/src/python_bindings/medoids_python.cpp b/src/python_bindings/medoids_python.cpp
index 155032dc..08b2f164 100644
--- a/src/python_bindings/medoids_python.cpp
+++ b/src/python_bindings/medoids_python.cpp
@@ -14,17 +14,18 @@
 #include "kmedoids_pywrapper.hpp"
 
 namespace km {
-pybind11::array_t<arma::uword> km::KMedoidsWrapper::getMedoidsFinalPython() {
-  if (KMedoids::getMedoidsFinal().size() > 1) {
-    return carma::row_to_arr<arma::uword>(
-      KMedoids::getMedoidsFinal()).squeeze();
-  } else {
-    return carma::row_to_arr<arma::uword>(KMedoids::getMedoidsFinal());
+  pybind11::array_t <arma::uword>
+  km::KMedoidsWrapper::getMedoidsFinalPython() {
+    if (KMedoids::getMedoidsFinal().size() > 1) {
+      return carma::row_to_arr<arma::uword>(
+              KMedoids::getMedoidsFinal()).squeeze();
+    } else {
+      return carma::row_to_arr<arma::uword>(KMedoids::getMedoidsFinal());
+    }
   }
-}
 
-void medoids_python(pybind11::class_<km::KMedoidsWrapper> *cls) {
-  cls->def_property_readonly("medoids",
-    &KMedoidsWrapper::getMedoidsFinalPython);
-}
+  void medoids_python(pybind11::class_ <km::KMedoidsWrapper> *cls) {
+    cls->def_property_readonly("medoids",
+                               &KMedoidsWrapper::getMedoidsFinalPython);
+  }
 }  // namespace km
diff --git a/src/python_bindings/steps_python.cpp b/src/python_bindings/steps_python.cpp
index 22e3ae77..fda4f512 100644
--- a/src/python_bindings/steps_python.cpp
+++ b/src/python_bindings/steps_python.cpp
@@ -14,11 +14,11 @@
 #include "kmedoids_pywrapper.hpp"
 
 namespace km {
-int km::KMedoidsWrapper::getStepsPython() {
-  return KMedoids::getSteps();
-}
+  int km::KMedoidsWrapper::getStepsPython() {
+    return KMedoids::getSteps();
+  }
 
-void steps_python(pybind11::class_<KMedoidsWrapper> *cls) {
-  cls->def_property_readonly("steps", &KMedoidsWrapper::getStepsPython);
-}
+  void steps_python(pybind11::class_ <KMedoidsWrapper> *cls) {
+    cls->def_property_readonly("steps", &KMedoidsWrapper::getStepsPython);
+  }
 }  // namespace km
diff --git a/src/python_bindings/swap_times_python.cpp b/src/python_bindings/swap_times_python.cpp
index 3901a04b..a2dbf9c4 100644
--- a/src/python_bindings/swap_times_python.cpp
+++ b/src/python_bindings/swap_times_python.cpp
@@ -13,23 +13,23 @@
 #include "kmedoids_pywrapper.hpp"
 
 namespace km {
-size_t km::KMedoidsWrapper::getTotalSwapTimePython() {
-  return KMedoids::getTotalSwapTime();
-}
+  size_t km::KMedoidsWrapper::getTotalSwapTimePython() {
+    return KMedoids::getTotalSwapTime();
+  }
 
-float km::KMedoidsWrapper::getTimePerSwapPython() {
-  return KMedoids::getTimePerSwap();
-}
+  float km::KMedoidsWrapper::getTimePerSwapPython() {
+    return KMedoids::getTimePerSwap();
+  }
 
-void total_swap_time_python(pybind11::class_<km::KMedoidsWrapper> *cls) {
-  cls->def_property_readonly(
-    "total_swap_time",
-    &KMedoidsWrapper::getTotalSwapTimePython);
-}
+  void total_swap_time_python(pybind11::class_ <km::KMedoidsWrapper> *cls) {
+    cls->def_property_readonly(
+            "total_swap_time",
+            &KMedoidsWrapper::getTotalSwapTimePython);
+  }
 
-void time_per_swap_python(pybind11::class_<km::KMedoidsWrapper> *cls) {
-  cls->def_property_readonly(
-    "time_per_swap",
-    &KMedoidsWrapper::getTimePerSwapPython);
-}
+  void time_per_swap_python(pybind11::class_ <km::KMedoidsWrapper> *cls) {
+    cls->def_property_readonly(
+            "time_per_swap",
+            &KMedoidsWrapper::getTimePerSwapPython);
+  }
 }  // namespace km


Test patch:
diff --git a/.github/workflows/run_linux_tests.yml b/.github/workflows/run_linux_tests.yml
index cfb9b9b1..01091846 100644
--- a/.github/workflows/run_linux_tests.yml
+++ b/.github/workflows/run_linux_tests.yml
@@ -24,8 +24,9 @@ jobs:
       - name: Install Python dependencies 
         run: |
           python -m pip install --upgrade pip
-          pip install pytest
-          pip install -r requirements.txt
+          python -m pip install pytest
+          python -m pip install -r requirements.txt
+          python -m pip install --upgrade pip setuptools wheel
 
       - name: Install Armadillo 10.7.5+
         run: |
@@ -50,7 +51,7 @@ jobs:
         run: |
           # The flags are necessary to ignore the pyproject.toml
           # See https://github.com/pypa/pip/issues/9738
-          pip install --no-use-pep517 --no-build-isolation -vvv -e .
+          python -m pip install --no-use-pep517 --no-build-isolation -vvv -e .
         env:
           # The default compiler on the Github Ubuntu runners is gcc
           # Would need to make a respective include change for clang 
@@ -76,9 +77,16 @@ jobs:
       #   run : |
       #     pytest tests/test_larger.py
 
+      - name: Clone for PR
+        if: ${{ github.event_name == 'pull_request' }}
+        run: git clone -b $GITHUB_HEAD_REF https://github.com/motiwari/BanditPAM
+
+      - name: Clone for Push
+        if: ${{ github.event_name != 'pull_request' }}
+        run: git clone -b ${GITHUB_REF#refs/heads/} https://github.com/motiwari/BanditPAM
+
       - name: Verify that the C++ executable compiles and runs
         run : |
-          git clone -b ${GITHUB_REF#refs/heads/} https://github.com/motiwari/BanditPAM
           cd BanditPAM
           mkdir build
           cd build
diff --git a/.github/workflows/run_mac_tests.yml b/.github/workflows/run_mac_tests.yml
index d7ab864f..fda3f5d1 100644
--- a/.github/workflows/run_mac_tests.yml
+++ b/.github/workflows/run_mac_tests.yml
@@ -25,8 +25,9 @@ jobs:
       - name: Install Python dependencies
         run: |
           python -m pip install --upgrade pip
-          pip install pytest
-          pip install -r requirements.txt
+          python -m pip install pytest
+          python -m pip install -r requirements.txt
+          python -m pip install --upgrade pip setuptools wheel
 
       - name: Install Armadillo 10.7.5+
         run: |
@@ -51,7 +52,7 @@ jobs:
         run: |
           # The flags are necessary to ignore the pyproject.toml
           # See https://github.com/pypa/pip/issues/9738
-          pip install --no-use-pep517 --no-build-isolation -vvv -e .
+          python -m pip install --no-use-pep517 --no-build-isolation -vvv -e .
         env:
           # The default compiler on the Github Ubuntu runners is gcc
           # Would need to make a respective include change for clang
diff --git a/.github/workflows/run_windows_tests.yml b/.github/workflows/run_windows_tests.yml
index a3f41979..b09d693d 100644
--- a/.github/workflows/run_windows_tests.yml
+++ b/.github/workflows/run_windows_tests.yml
@@ -25,8 +25,8 @@ jobs:
       - name: Install Python dependencies 
         run: |
           python -m pip install --upgrade pip
-          pip install pytest
-          pip install -r requirements.txt
+          python -m pip install pytest
+          python -m pip install -r requirements.txt
       - name: Install Armadillo 10.7.5+
         run: |
           cd ~
@@ -45,7 +45,7 @@ jobs:
           cd ~
       - name: Install BanditPAM package
         run: |
-          pip install -vvv -e .
+          python -m pip install -vvv -e .
         env:
           # The default compiler on the Github Ubuntu runners is gcc
           # Would need to make a respective include change for clang 
diff --git a/tests/test_initialization.py b/tests/test_initialization.py
index ad48d60b..3689d787 100644
--- a/tests/test_initialization.py
+++ b/tests/test_initialization.py
@@ -10,7 +10,7 @@ def test_default_initialization(self):
         assert kmed.algorithm == "BanditPAM", "algorithm does not match"
         assert kmed.max_iter == 100, "max_iter does not match"
         assert kmed.build_confidence == 3, "build_confidence does not match"
-        assert kmed.swap_confidence == 4, "swap_confidence does not match"
+        assert kmed.swap_confidence == 5, "swap_confidence does not match"
         assert kmed.use_cache, "use_cache does not match"
         assert kmed.use_perm, "use_perm does not match"
         assert kmed.cache_width == 1000, "cache_width does not match"
@@ -78,8 +78,8 @@ def test_initialization_3(self):
         assert kmed.n_medoids == 100, "n_medoids does not match"
         assert kmed.algorithm == "PAM", "algorithm does not match"
         assert kmed.max_iter == 10000, "max_iter does not match"
-        assert kmed.build_confidence == 11000, "build_confidence does not" \
-                                               " match"
+        assert kmed.build_confidence == 11000, "build_confidence does not " \
+                                               "match"
         assert kmed.swap_confidence == 2029, "swap_confidence does not match"
         assert kmed.use_cache, "use_cache does not match"
         assert not kmed.use_perm, "use_perm does not match"