From df5e3a94893a1bd557321d20410354947feee2ad Mon Sep 17 00:00:00 2001
From: hmacdope <hugomacdermott@gmail.com>
Date: Fri, 6 Sep 2024 16:02:17 +1000
Subject: [PATCH 01/12] add

---
 distopia/_distopia.pyx | 93 +++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 92 insertions(+), 1 deletion(-)

diff --git a/distopia/_distopia.pyx b/distopia/_distopia.pyx
index 5b0d077b..dbed1d61 100644
--- a/distopia/_distopia.pyx
+++ b/distopia/_distopia.pyx
@@ -129,6 +129,95 @@ cdef extern from "distopia.h" namespace "distopia" nogil:
         const T *box,
         T *out,
     )
+    void CalcBondsNoBoxIdx[T](
+        const T *coords,
+        const int* a_idx,
+        const int* b_idx,
+        size_t n,
+        T *out,
+    )
+
+    void CalcBondsOrthoIdx[T](
+        const T *coords,
+        const int* a_idx,
+        const int* b_idx,
+        size_t n,
+        const T *box,
+        T *out,
+    )
+
+    void CalcBondsTriclinicIdx[T](
+        const T *coords,
+        const int* a_idx,
+        const int* b_idx,
+        size_t n,
+        const T *box,
+        T *out,
+    )
+
+    void CalcAnglesNoBoxIdx[T](
+        const T *coords,
+        const int* a_idx,
+        const int* b_idx,
+        const int* c_idx,
+        size_t n,
+        T *out,
+    )
+
+    void CalcAnglesOrthoIdx[T](
+        const T *coords,
+        const int* a_idx,
+        const int* b_idx,
+        const int* c_idx,
+        size_t n,
+        const T *box,
+        T *out,
+    )
+
+    void CalcAnglesTriclinicIdx[T](
+        const T *coords,
+        const int* a_idx,
+        const int* b_idx,
+        const int* c_idx,
+        size_t n,
+        const T *box,
+        T *out,
+    )
+
+    void CalcDihedralsNoBoxIdx[T](
+        const T *coords,
+        const int* a_idx,
+        const int* b_idx,
+        const int* c_idx,
+        const int* d_idx,
+        size_t n,
+        T *out,
+    )
+
+    void CalcDihedralsOrthoIdx[T](
+        const T *coords,
+        const int* a_idx,
+        const int* b_idx,
+        const int* c_idx,
+        const int* d_idx,
+        size_t n,
+        const T *box,
+        T *out,
+    )
+
+
+    void CalcDihedralsTriclinicIdx[T](
+        const T *coords,
+        const int* a_idx,
+        const int* b_idx,
+        const int* c_idx,
+        const int* d_idx,
+        size_t n,
+        const T *box,
+        T *out,
+    )
+
+
 
 def get_n_float_lanes():
     """The number of floats per register distopia will handle on this system"""
@@ -890,4 +979,6 @@ def calc_self_distance_array_triclinic(
                                &box[0][0],
                                &results_view[0])
 
-    return np.array(results).reshape(coords0.shape[0], coords0.shape[0])
\ No newline at end of file
+    return np.array(results).reshape(coords0.shape[0], coords0.shape[0])
+
+

From 9ec1dd280f4b9af61de2a482d1e3a86f33b121d0 Mon Sep 17 00:00:00 2001
From: hmacdope <hugomacdermott@gmail.com>
Date: Fri, 6 Sep 2024 16:17:06 +1000
Subject: [PATCH 02/12] add

---
 distopia/__init__.py            |   6 ++
 distopia/_distopia.pyx          | 144 ++++++++++++++++++++++++++++++++
 distopia/tests/test_distopia.py |  35 +++++++-
 3 files changed, 184 insertions(+), 1 deletion(-)

diff --git a/distopia/__init__.py b/distopia/__init__.py
index 6e07abf3..8a062d6d 100644
--- a/distopia/__init__.py
+++ b/distopia/__init__.py
@@ -16,6 +16,9 @@
     calc_self_distance_array_no_box,
     calc_self_distance_array_ortho,
     calc_self_distance_array_triclinic,
+    calc_bonds_no_box_idx,
+    calc_bonds_ortho_idx,
+    calc_bonds_triclinic_idx,
 )
 
 __all__ = [
@@ -34,5 +37,8 @@
     'calc_self_distance_array_no_box',
     'calc_self_distance_array_ortho',
     'calc_self_distance_array_triclinic',
+    'calc_bonds_no_box_idx',
+    'calc_bonds_ortho_idx',
+    'calc_bonds_triclinic_idx',
     '__version__',
 ]
\ No newline at end of file
diff --git a/distopia/_distopia.pyx b/distopia/_distopia.pyx
index dbed1d61..b40a6c51 100644
--- a/distopia/_distopia.pyx
+++ b/distopia/_distopia.pyx
@@ -982,3 +982,147 @@ def calc_self_distance_array_triclinic(
     return np.array(results).reshape(coords0.shape[0], coords0.shape[0])
 
 
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def calc_bonds_no_box_idx(
+     floating[:, ::1] coords,
+     int[::1] a_idx,
+     int[::1] b_idx,
+     floating[::1] results=None):
+    """Calculate pairwise distances between coords[a_idx] and coords[b_idx] with no periodic boundary conditions
+
+    Parameters
+    ----------
+    coords : float32 or float64 array
+      must be same length and dtype
+    a_idx, b_idx : int array
+      must be same length and dtype
+    results: float32 or float64 array (optional)
+      array to store results in, must be same size and dtype as a_idx/b_idx
+
+    Returns
+    -------
+    distances : float32 or float64 array
+      same length and dtype as a_idx/b_idx
+    """
+    cdef floating[::1] results_view
+    cdef size_t nvals = a_idx.shape[0]
+    cdef cnp.npy_intp[1] dims
+    dims[0] = <ssize_t > nvals  # FIXME truncation?
+
+    _check_shapes(a_idx, b_idx)
+
+    if results is None:
+        if floating is float:
+            results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0)
+        else:
+            results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT64, 0)
+
+    else:
+        _check_results(results, nvals)
+
+    results_view = results
+
+    CalcBondsNoBoxIdx(& coords[0][0], & a_idx[0], & b_idx[0], nvals, & results_view[0])
+
+    return np.array(results)
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def calc_bonds_ortho_idx(
+     floating[:, ::1] coords,
+     int[::1] a_idx,
+     int[::1] b_idx,
+     floating[::1] box,
+     floating[::1] results=None):
+    """Calculate pairwise distances between coords[a_idx] and coords[b_idx] under orthorhombic boundary conditions
+
+    Parameters
+    ----------
+    coords : float32 or float64 array
+      must be same length and dtype
+    a_idx, b_idx : int array
+      must be same length and dtype
+    box : float32 or float64 array
+      orthorhombic periodic boundary dimensions in [L, L, L] format
+    results: float32 or float64 array (optional)
+      array to store results in, must be same size and dtype as a_idx/b_idx
+
+    Returns
+    -------
+    distances : float32 or float64 array
+      same length and dtype as a_idx/b_idx
+    """
+    cdef floating[::1] results_view
+    cdef size_t nvals = a_idx.shape[0]
+    cdef cnp.npy_intp[1] dims
+    dims[0] = <ssize_t > nvals  # FIXME truncation?
+
+    _check_shapes(a_idx, b_idx)
+
+    if results is None:
+        if floating is float:
+            results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0)
+        else:
+            results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT64, 0)
+
+    else:
+        _check_results(results, nvals)
+
+    results_view = results
+
+    CalcBondsOrthoIdx(& coords[0][0], & a_idx[0], & b_idx[0], nvals, & box[0], & results_view[0])
+
+    return np.array(results)
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def calc_bonds_triclinic_idx(
+     floating[:, ::1] coords,
+     int[::1] a_idx,
+     int[::1] b_idx,
+     floating[:, ::1] box,
+     floating[::1] results=None):
+    """Calculate pairwise distances between coords[a_idx] and coords[b_idx] under triclinic boundary conditions
+
+    Parameters
+    ----------
+    coords : float32 or float64 array
+      must be same length and dtype
+    a_idx, b_idx : int array
+      must be same length and dtype
+    box : float32 or float64 array
+      periodic boundary dimensions, in 3x3 format
+    results: float32 or float64 array (optional)
+      array to store results in, must be same size and dtype as a_idx/b_idx
+
+    Returns
+    -------
+    distances : float32 or float64 array
+      same length and dtype as a_idx/b_idx
+    """
+    cdef floating[::1] results_view
+    cdef size_t nvals = a_idx.shape[0]
+    cdef cnp.npy_intp[1] dims
+    dims[0] = <ssize_t > nvals  # FIXME truncation?
+
+    _check_shapes(a_idx, b_idx)
+
+    if results is None:
+        if floating is float:
+            results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0)
+        else:
+            results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT64, 0)
+
+    else:
+        _check_results(results, nvals)
+
+    results_view = results
+
+    CalcBondsTriclinicIdx(& coords[0][0], & a_idx[0], & b_idx[0], nvals, & box[0][0], & results_view[0])
+
+    return np.array(results)
+
diff --git a/distopia/tests/test_distopia.py b/distopia/tests/test_distopia.py
index 798780d8..37824b8b 100644
--- a/distopia/tests/test_distopia.py
+++ b/distopia/tests/test_distopia.py
@@ -240,6 +240,35 @@ def test_triclinic_bad_result(self):
             distopia.calc_self_distance_array_triclinic(c0, box, results=np.empty(1, dtype=np.float32))
 
 
+
+class TestDistancesIdx:
+    def arange_input(self, N, dtype):
+        return np.arange(3 * N, dtype=dtype).reshape(N, 3)
+
+    def result_shim(self, make_arr, N, dtype):
+        if not make_arr:
+            return None
+        else:
+            return np.empty(N, dtype=dtype)
+
+
+    @pytest.mark.parametrize("dtype", (np.float32, np.float64))
+    @pytest.mark.parametrize("N", (0, 10, 1000, 10000))
+    @pytest.mark.parametrize("use_result_buffer", (True, False))
+    def test_calc_bonds_no_box_idx(self, N, use_result_buffer, dtype):
+        c = self.arange_input(N, dtype)
+        a_idx = np.asarray([0, 1, 2, 3, 4, 5])
+        b_idx = np.asarray([6, 7, 8, 9, 10, 11])
+        result_buffer = self.result_shim(use_result_buffer, N, dtype)
+        result = distopia.calc_bonds_no_box_idx(
+            c, a_idx, b_idx, results=result_buffer
+        )
+        assert_allclose(result, np.zeros(N))
+        # check dtype of result
+        assert result.dtype == dtype
+    
+
+
 class TestMDA:
     """
     Copy of some of the tests from MDAnalysisTests repo
@@ -381,4 +410,8 @@ def test_periodic_angles(self, box_bonds, positions_angles, dtype):
         test3 = distopia.calc_angles_ortho(a, b, c2, box=box)
 
         for val in [test1, test2, test3]:
-            assert_almost_equal(ref, val, self.prec, err_msg="Min image in angle calculation failed")
\ No newline at end of file
+            assert_almost_equal(ref, val, self.prec, err_msg="Min image in angle calculation failed")
+
+
+
+

From e13a092d57b226b874b7dee99cb655b457ffb018 Mon Sep 17 00:00:00 2001
From: hmacdope <hugomacdermott@gmail.com>
Date: Fri, 6 Sep 2024 16:30:11 +1000
Subject: [PATCH 03/12] use integral

---
 distopia/_distopia.pyx | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/distopia/_distopia.pyx b/distopia/_distopia.pyx
index b40a6c51..5988972f 100644
--- a/distopia/_distopia.pyx
+++ b/distopia/_distopia.pyx
@@ -13,7 +13,7 @@ The python functions for distopia
 import numpy as np
 cimport cython
 cimport numpy as cnp
-from cython cimport floating
+from cython cimport floating, integral
 cnp.import_array()
 
 cdef extern from "distopia.h" namespace "distopia" nogil:
@@ -987,8 +987,8 @@ def calc_self_distance_array_triclinic(
 @cython.wraparound(False)
 def calc_bonds_no_box_idx(
      floating[:, ::1] coords,
-     int[::1] a_idx,
-     int[::1] b_idx,
+     integral[::1] a_idx,
+     integral[::1] b_idx,
      floating[::1] results=None):
     """Calculate pairwise distances between coords[a_idx] and coords[b_idx] with no periodic boundary conditions
 
@@ -1033,8 +1033,8 @@ def calc_bonds_no_box_idx(
 @cython.wraparound(False)
 def calc_bonds_ortho_idx(
      floating[:, ::1] coords,
-     int[::1] a_idx,
-     int[::1] b_idx,
+     integral[::1] a_idx,
+     integral[::1] b_idx,
      floating[::1] box,
      floating[::1] results=None):
     """Calculate pairwise distances between coords[a_idx] and coords[b_idx] under orthorhombic boundary conditions
@@ -1082,8 +1082,8 @@ def calc_bonds_ortho_idx(
 @cython.wraparound(False)
 def calc_bonds_triclinic_idx(
      floating[:, ::1] coords,
-     int[::1] a_idx,
-     int[::1] b_idx,
+     integral[::1] a_idx,
+     integral[::1] b_idx,
      floating[:, ::1] box,
      floating[::1] results=None):
     """Calculate pairwise distances between coords[a_idx] and coords[b_idx] under triclinic boundary conditions

From 2fe6b65554bdec4945804c1a90d94747dcbca5a6 Mon Sep 17 00:00:00 2001
From: hmacdope <hugomacdermott@gmail.com>
Date: Fri, 6 Sep 2024 16:48:59 +1000
Subject: [PATCH 04/12] segfaulting test

---
 distopia/_distopia.pyx          | 15 +++++++++------
 distopia/tests/test_distopia.py |  8 ++++----
 2 files changed, 13 insertions(+), 10 deletions(-)

diff --git a/distopia/_distopia.pyx b/distopia/_distopia.pyx
index 5988972f..7ae4bac6 100644
--- a/distopia/_distopia.pyx
+++ b/distopia/_distopia.pyx
@@ -987,8 +987,8 @@ def calc_self_distance_array_triclinic(
 @cython.wraparound(False)
 def calc_bonds_no_box_idx(
      floating[:, ::1] coords,
-     integral[::1] a_idx,
-     integral[::1] b_idx,
+     int[::1] a_idx,
+     int[::1] b_idx,
      floating[::1] results=None):
     """Calculate pairwise distances between coords[a_idx] and coords[b_idx] with no periodic boundary conditions
 
@@ -1013,6 +1013,7 @@ def calc_bonds_no_box_idx(
 
     _check_shapes(a_idx, b_idx)
 
+
     if results is None:
         if floating is float:
             results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0)
@@ -1033,8 +1034,8 @@ def calc_bonds_no_box_idx(
 @cython.wraparound(False)
 def calc_bonds_ortho_idx(
      floating[:, ::1] coords,
-     integral[::1] a_idx,
-     integral[::1] b_idx,
+     int[::1] a_idx,
+     int[::1] b_idx,
      floating[::1] box,
      floating[::1] results=None):
     """Calculate pairwise distances between coords[a_idx] and coords[b_idx] under orthorhombic boundary conditions
@@ -1062,6 +1063,7 @@ def calc_bonds_ortho_idx(
 
     _check_shapes(a_idx, b_idx)
 
+
     if results is None:
         if floating is float:
             results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0)
@@ -1082,8 +1084,8 @@ def calc_bonds_ortho_idx(
 @cython.wraparound(False)
 def calc_bonds_triclinic_idx(
      floating[:, ::1] coords,
-     integral[::1] a_idx,
-     integral[::1] b_idx,
+     int[::1] a_idx,
+     int[::1] b_idx,
      floating[:, ::1] box,
      floating[::1] results=None):
     """Calculate pairwise distances between coords[a_idx] and coords[b_idx] under triclinic boundary conditions
@@ -1111,6 +1113,7 @@ def calc_bonds_triclinic_idx(
 
     _check_shapes(a_idx, b_idx)
 
+
     if results is None:
         if floating is float:
             results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0)
diff --git a/distopia/tests/test_distopia.py b/distopia/tests/test_distopia.py
index 37824b8b..d2a2c696 100644
--- a/distopia/tests/test_distopia.py
+++ b/distopia/tests/test_distopia.py
@@ -253,13 +253,13 @@ def result_shim(self, make_arr, N, dtype):
 
 
     @pytest.mark.parametrize("dtype", (np.float32, np.float64))
-    @pytest.mark.parametrize("N", (0, 10, 1000, 10000))
+    @pytest.mark.parametrize("N", (10, 1000, 10000))
     @pytest.mark.parametrize("use_result_buffer", (True, False))
     def test_calc_bonds_no_box_idx(self, N, use_result_buffer, dtype):
         c = self.arange_input(N, dtype)
-        a_idx = np.asarray([0, 1, 2, 3, 4, 5])
-        b_idx = np.asarray([6, 7, 8, 9, 10, 11])
-        result_buffer = self.result_shim(use_result_buffer, N, dtype)
+        a_idx = np.asarray(np.arange(N -3)).astype(np.int32)
+        b_idx = np.asarray(np.arange(N -3)).astype(np.int32)
+        result_buffer = self.result_shim(use_result_buffer, len(a_idx), dtype)
         result = distopia.calc_bonds_no_box_idx(
             c, a_idx, b_idx, results=result_buffer
         )

From ee0d160b917b7b48e436e6478a7f089cd283b7fb Mon Sep 17 00:00:00 2001
From: hmacdope <hugomacdermott@gmail.com>
Date: Fri, 6 Sep 2024 18:00:23 +1000
Subject: [PATCH 05/12] everything

---
 distopia/_distopia.pyx              |  1 -
 distopia/tests/test_distopia.py     |  3 ++-
 libdistopia/test/test_mda_match.cpp | 19 +++++++++++++++++++
 3 files changed, 21 insertions(+), 2 deletions(-)

diff --git a/distopia/_distopia.pyx b/distopia/_distopia.pyx
index 7ae4bac6..78519701 100644
--- a/distopia/_distopia.pyx
+++ b/distopia/_distopia.pyx
@@ -1013,7 +1013,6 @@ def calc_bonds_no_box_idx(
 
     _check_shapes(a_idx, b_idx)
 
-
     if results is None:
         if floating is float:
             results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0)
diff --git a/distopia/tests/test_distopia.py b/distopia/tests/test_distopia.py
index d2a2c696..11683d24 100644
--- a/distopia/tests/test_distopia.py
+++ b/distopia/tests/test_distopia.py
@@ -253,10 +253,11 @@ def result_shim(self, make_arr, N, dtype):
 
 
     @pytest.mark.parametrize("dtype", (np.float32, np.float64))
-    @pytest.mark.parametrize("N", (10, 1000, 10000))
+    @pytest.mark.parametrize("N", (3, 6, 10, 1000, 10000))
     @pytest.mark.parametrize("use_result_buffer", (True, False))
     def test_calc_bonds_no_box_idx(self, N, use_result_buffer, dtype):
         c = self.arange_input(N, dtype)
+
         a_idx = np.asarray(np.arange(N -3)).astype(np.int32)
         b_idx = np.asarray(np.arange(N -3)).astype(np.int32)
         result_buffer = self.result_shim(use_result_buffer, len(a_idx), dtype)
diff --git a/libdistopia/test/test_mda_match.cpp b/libdistopia/test/test_mda_match.cpp
index e226ec6c..023366fd 100644
--- a/libdistopia/test/test_mda_match.cpp
+++ b/libdistopia/test/test_mda_match.cpp
@@ -1,6 +1,7 @@
 #include <cmath>
 #include <iostream>
 #include <numeric>
+#include <stdlib.h>     
 
 #include "gtest/gtest.h"
 #include "distopia.h"
@@ -502,6 +503,24 @@ TYPED_TEST(CoordinatesIdx, CalcBondsNoBoxIdx) {
     }
 }
 
+TYPED_TEST(CoordinatesIdx, CalcBondsNoBoxIdxSmall) {
+    int ncoords = 250;
+    // 3 is always a wierd remainder case due to 
+    int nidx = 3;
+
+    this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
+
+    // reference result
+    distopia::CalcBondsNoBox(this->a_coords_contig, this->b_coords_contig, this->nidx, this->ref_results);
+    std::cout << "ref_results: " << std::endl;
+    // test the idx
+    distopia::CalcBondsNoBoxIdx(this->coords, this->a_idx, this->b_idx, this->nidx, this->results);
+
+    for (std::size_t i=0; i<this->nidx; i++) {
+        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+    }
+}
+
 
 TYPED_TEST(CoordinatesIdx, CalcBondsOrthoIdx) {
     int ncoords = 250;

From 3d3267eaeb9fe8aa8f832a8b6a9ff7428ad40447 Mon Sep 17 00:00:00 2001
From: hmacdope <hugomacdermott@gmail.com>
Date: Fri, 6 Sep 2024 20:22:58 +1000
Subject: [PATCH 06/12] dispatch

---
 libdistopia/src/distopia.cpp | 82 ++++++++++++++++++++++++++++++++++--
 1 file changed, 79 insertions(+), 3 deletions(-)

diff --git a/libdistopia/src/distopia.cpp b/libdistopia/src/distopia.cpp
index b1187fa0..7bff6f1a 100644
--- a/libdistopia/src/distopia.cpp
+++ b/libdistopia/src/distopia.cpp
@@ -1607,74 +1607,132 @@ namespace distopia {
     }
     HWY_DLLEXPORT template <> void CalcBondsNoBoxIdx(const float *coords, const int *a_idx, const int *b_idx,
                                                      int n, float *out) {
+
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcBondsNoBoxIdxSingle(coords, a_idx, b_idx, n, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcBondsNoBoxIdxSingle)(coords, a_idx, b_idx, n, out);
     }
     HWY_DLLEXPORT template <> void CalcBondsNoBoxIdx(const double *coords, const int *a_idx, const int *b_idx,
                                                      int n, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcBondsNoBoxIdxDouble(coords, a_idx, b_idx, n, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcBondsNoBoxIdxDouble)(coords, a_idx, b_idx, n, out);
     }
     HWY_DLLEXPORT template <> void CalcBondsOrthoIdx(const float *coords, const int *a_idx, const int *b_idx,
                                                      int n, const float *box, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcBondsOrthoIdxSingle(coords, a_idx, b_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcBondsOrthoIdxSingle)(coords, a_idx, b_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcBondsOrthoIdx(const double *coords, const int *a_idx, const int *b_idx,
                                                      int n, const double *box, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcBondsOrthoIdxDouble(coords, a_idx, b_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcBondsOrthoIdxDouble)(coords, a_idx, b_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcBondsTriclinicIdx(const float *coords, const int *a_idx, const int *b_idx,
                                                          int n, const float *box, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcBondsTriclinicIdxSingle(coords, a_idx, b_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcBondsTriclinicIdxSingle)(coords, a_idx, b_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcBondsTriclinicIdx(const double *coords, const int *a_idx, const int *b_idx,
                                                          int n, const double *box, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcBondsTriclinicIdxDouble(coords, a_idx, b_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcBondsTriclinicIdxDouble)(coords, a_idx, b_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcAnglesNoBoxIdx(const float *coords, const int *a_idx, const int *b_idx, const int *c_idx,
                                                       int n, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcAnglesNoBoxIdxSingle(coords, a_idx, b_idx, c_idx, n, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcAnglesNoBoxIdxSingle)(coords, a_idx, b_idx, c_idx, n, out);
     }
     HWY_DLLEXPORT template <> void CalcAnglesNoBoxIdx(const double *coords, const int *a_idx, const int *b_idx, const int *c_idx,
                                                       int n, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcAnglesNoBoxIdxDouble(coords, a_idx, b_idx, c_idx, n, out);
+        }
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcAnglesNoBoxIdxDouble(coords, a_idx, b_idx, c_idx, n, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcAnglesNoBoxIdxDouble)(coords, a_idx, b_idx, c_idx, n, out);
     }
     HWY_DLLEXPORT template <> void CalcAnglesOrthoIdx(const float *coords, const int *a_idx, const int *b_idx, const int *c_idx,
                                                       int n, const float *box, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcAnglesOrthoIdxSingle(coords, a_idx, b_idx, c_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcAnglesOrthoIdxSingle)(coords, a_idx, b_idx, c_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcAnglesOrthoIdx(const double *coords, const int *a_idx, const int *b_idx, const int *c_idx,
                                                       int n, const double *box, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcAnglesOrthoIdxDouble(coords, a_idx, b_idx, c_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcAnglesOrthoIdxDouble)(coords, a_idx, b_idx, c_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcAnglesTriclinicIdx(const float *coords, const int *a_idx, const int *b_idx, const int *c_idx,
                                                           int n, const float *box, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcAnglesTriclinicIdxSingle(coords, a_idx, b_idx, c_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcAnglesTriclinicIdxSingle)(coords, a_idx, b_idx, c_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcAnglesTriclinicIdx(const double *coords, const int *a_idx, const int *b_idx, const int *c_idx,
                                                           int n, const double *box, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcAnglesTriclinicIdxDouble(coords, a_idx, b_idx, c_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcAnglesTriclinicIdxDouble)(coords, a_idx, b_idx, c_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcDihedralsNoBoxIdx(const float *coords, const int *a_idx, const int *b_idx, const int *c_idx, const int *d_idx,
                                                          int n, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcDihedralsNoBoxIdxSingle(coords, a_idx, b_idx, c_idx, d_idx, n, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcDihedralsNoBoxIdxSingle)(coords, a_idx, b_idx, c_idx, d_idx, n, out);
     }
     HWY_DLLEXPORT template <> void CalcDihedralsNoBoxIdx(const double *coords, const int *a_idx, const int *b_idx, const int *c_idx, const int *d_idx,
                                                          int n, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcDihedralsNoBoxIdxDouble(coords, a_idx, b_idx, c_idx, d_idx, n, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcDihedralsNoBoxIdxDouble)(coords, a_idx, b_idx, c_idx, d_idx, n, out);
     }
     HWY_DLLEXPORT template <> void CalcDihedralsOrthoIdx(const float *coords, const int *a_idx, const int *b_idx, const int *c_idx, const int *d_idx,
                                                          int n, const float *box, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcDihedralsOrthoIdxSingle(coords, a_idx, b_idx, c_idx, d_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcDihedralsOrthoIdxSingle)(coords, a_idx, b_idx, c_idx, d_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcDihedralsOrthoIdx(const double *coords, const int *a_idx, const int *b_idx, const int *c_idx, const int *d_idx,
                                                          int n, const double *box, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcDihedralsOrthoIdxDouble(coords, a_idx, b_idx, c_idx, d_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcDihedralsOrthoIdxDouble)(coords, a_idx, b_idx, c_idx, d_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcDihedralsTriclinicIdx(const float *coords, const int *a_idx, const int *b_idx, const int *c_idx, const int *d_idx,
                                                              int n, const float *box, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcDihedralsTriclinicIdxSingle(coords, a_idx, b_idx, c_idx, d_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcDihedralsTriclinicIdxSingle)(coords, a_idx, b_idx, c_idx, d_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcDihedralsTriclinicIdx(const double *coords, const int *a_idx, const int *b_idx, const int *c_idx, const int *d_idx,
                                                              int n, const double *box, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcDihedralsTriclinicIdxDouble(coords, a_idx, b_idx, c_idx, d_idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcDihedralsTriclinicIdxDouble)(coords, a_idx, b_idx, c_idx, d_idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcDistanceArrayNoBoxIdx(const float *coords, const int *a_idx, const int *b_idx, int na, int nb, float *out) {
@@ -1684,7 +1742,7 @@ namespace distopia {
         return HWY_DYNAMIC_DISPATCH(CalcDistanceArrayNoBoxIdxSingle)(coords, a_idx, b_idx, na, nb, out);
     }
     HWY_DLLEXPORT template <> void CalcDistanceArrayNoBoxIdx(const double *coords, const int *a_idx, const int *b_idx, int na, int nb, double *out) {
-        if (nb < GetNFloatLanes()) {
+        if (nb < GetNDoubleLanes()) {
             return distopia::N_SCALAR::CalcDistanceArrayNoBoxIdxDouble(coords, a_idx, b_idx, na, nb, out);
         }
         return HWY_DYNAMIC_DISPATCH(CalcDistanceArrayNoBoxIdxDouble)(coords, a_idx, b_idx, na, nb, out);
@@ -1696,7 +1754,7 @@ namespace distopia {
         return HWY_DYNAMIC_DISPATCH(CalcDistanceArrayOrthoIdxSingle)(coords, a_idx, b_idx, na, nb, box, out);
     }
     HWY_DLLEXPORT template <> void CalcDistanceArrayOrthoIdx(const double *coords, const int *a_idx, const int *b_idx, int na, int nb, const double *box, double *out) {
-        if (nb < GetNFloatLanes()) {
+        if (nb < GetNDoubleLanes()) {
             return distopia::N_SCALAR::CalcDistanceArrayOrthoIdxDouble(coords, a_idx, b_idx, na, nb, box, out);
         }
         return HWY_DYNAMIC_DISPATCH(CalcDistanceArrayOrthoIdxDouble)(coords, a_idx, b_idx, na, nb, box, out);
@@ -1708,27 +1766,45 @@ namespace distopia {
         return HWY_DYNAMIC_DISPATCH(CalcDistanceArrayTriclinicIdxSingle)(coords, a_idx, b_idx, na, nb, box, out);
     }
     HWY_DLLEXPORT template <> void CalcDistanceArrayTriclinicIdx(const double *coords, const int *a_idx, const int *b_idx, int na, int nb, const double *box, double *out) {
-        if (nb < GetNFloatLanes()) {
+        if (nb < GetNDoubleLanes()) {
             return distopia::N_SCALAR::CalcDistanceArrayTriclinicIdxDouble(coords, a_idx, b_idx, na, nb, box, out);
         }
         return HWY_DYNAMIC_DISPATCH(CalcDistanceArrayTriclinicIdxDouble)(coords, a_idx, b_idx, na, nb, box, out);
     }
     HWY_DLLEXPORT template <> void CalcSelfDistanceArrayNoBoxIdx(const float *coords, const int *idx, int n, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcSelfDistanceArrayNoBoxIdxSingle(coords, idx, n, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcSelfDistanceArrayNoBoxIdxSingle)(coords, idx, n, out);
     }
     HWY_DLLEXPORT template <> void CalcSelfDistanceArrayNoBoxIdx(const double *coords, const int *idx, int n, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcSelfDistanceArrayNoBoxIdxDouble(coords, idx, n, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcSelfDistanceArrayNoBoxIdxDouble)(coords, idx, n, out);
     }
     HWY_DLLEXPORT template <> void CalcSelfDistanceArrayOrthoIdx(const float *coords, const int *idx, int n, const float *box, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcSelfDistanceArrayOrthoIdxSingle(coords, idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcSelfDistanceArrayOrthoIdxSingle)(coords, idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcSelfDistanceArrayOrthoIdx(const double *coords, const int *idx, int n, const double *box, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcSelfDistanceArrayOrthoIdxDouble(coords, idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcSelfDistanceArrayOrthoIdxDouble)(coords, idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcSelfDistanceArrayTriclinicIdx(const float *coords, const int *idx, int n, const float *box, float *out) {
+        if (n < GetNFloatLanes()) {
+            return distopia::N_SCALAR::CalcSelfDistanceArrayTriclinicIdxSingle(coords, idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcSelfDistanceArrayTriclinicIdxSingle)(coords, idx, n, box, out);
     }
     HWY_DLLEXPORT template <> void CalcSelfDistanceArrayTriclinicIdx(const double *coords, const int *idx, int n, const double *box, double *out) {
+        if (n < GetNDoubleLanes()) {
+            return distopia::N_SCALAR::CalcSelfDistanceArrayTriclinicIdxDouble(coords, idx, n, box, out);
+        }
         return HWY_DYNAMIC_DISPATCH(CalcSelfDistanceArrayTriclinicIdxDouble)(coords, idx, n, box, out);
     }
 

From ced356bf45f5c6709e74091822c0eeb83a7eab0d Mon Sep 17 00:00:00 2001
From: hmacdope <hugomacdermott@gmail.com>
Date: Sat, 14 Sep 2024 10:38:52 +1000
Subject: [PATCH 07/12] progress

---
 libdistopia/CMakeLists.txt          |  5 +++
 libdistopia/src/distopia.cpp        |  5 +++
 libdistopia/test/test_mda_match.cpp | 62 +++++++++++++++--------------
 libdistopia/test/test_utils.h       |  2 +
 4 files changed, 44 insertions(+), 30 deletions(-)

diff --git a/libdistopia/CMakeLists.txt b/libdistopia/CMakeLists.txt
index 23b9a45b..ef05c6d0 100644
--- a/libdistopia/CMakeLists.txt
+++ b/libdistopia/CMakeLists.txt
@@ -41,10 +41,14 @@ Include(GoogleTest)
 add_subdirectory("googletest")
 enable_testing()
 include_directories(${gtest_SOURCE_DIR}/include ${gtest_SOURCE_DIR})
+include_directories(${gmock_SOURCE_DIR}/include ${gmock_SOURCE_DIR})
+
 
 add_executable(test)
 target_sources(test PRIVATE "test/test.cpp")
 target_link_libraries(test PUBLIC gtest gtest_main)
+target_link_libraries(test PUBLIC gmock gmock_main)
+
 target_link_libraries(test PUBLIC libdistopia)
 target_include_directories(test PUBLIC  ${CMAKE_SOURCE_DIR})
 gtest_discover_tests(test WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY})
@@ -58,6 +62,7 @@ target_include_directories(targets PUBLIC  ${CMAKE_SOURCE_DIR})
 add_executable(test_mda_match)
 target_sources(test_mda_match PRIVATE "test/test_mda_match.cpp")
 target_link_libraries(test_mda_match PUBLIC gtest gtest_main)
+target_link_libraries(test PUBLIC gmock gmock_main)
 target_link_libraries(test_mda_match PUBLIC libdistopia)
 target_include_directories(test_mda_match PUBLIC  ${CMAKE_SOURCE_DIR})
 gtest_discover_tests(test_mda_match WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY})
diff --git a/libdistopia/src/distopia.cpp b/libdistopia/src/distopia.cpp
index 7bff6f1a..dfb5e6fc 100644
--- a/libdistopia/src/distopia.cpp
+++ b/libdistopia/src/distopia.cpp
@@ -759,6 +759,10 @@ namespace distopia {
             auto b_y = hn::Undefined(d);
             auto b_z = hn::Undefined(d);
 
+            std::cout << "nlanes: " << nlanes << std::endl;
+            std::cout << "n: " << n << std::endl;
+            std::cout << "n - nlanes: " << n - nlanes << std::endl;
+
             for (size_t i=0; i <= n - nlanes; i += nlanes) {
                 // load N indices of each source
                 // interleaved gather these indices
@@ -769,6 +773,7 @@ namespace distopia {
 
                 hn::StoreU(result, d, dst + i);
             }
+            std::cout << "n % nlanes: " << n % nlanes << std::endl;
             size_t rem = n % nlanes;
             if (rem) {  // if we had a non-multiple of nlanes, do final nlanes values again
                 LoadInterleavedIdx(a_idx + n - nlanes, coords, a_x, a_y, a_z);
diff --git a/libdistopia/test/test_mda_match.cpp b/libdistopia/test/test_mda_match.cpp
index 023366fd..6c95ee03 100644
--- a/libdistopia/test/test_mda_match.cpp
+++ b/libdistopia/test/test_mda_match.cpp
@@ -511,10 +511,10 @@ TYPED_TEST(CoordinatesIdx, CalcBondsNoBoxIdxSmall) {
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
     // reference result
-    distopia::CalcBondsNoBox(this->a_coords_contig, this->b_coords_contig, this->nidx, this->ref_results);
+    distopia::CalcBondsNoBox<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->nidx, this->ref_results);
     std::cout << "ref_results: " << std::endl;
     // test the idx
-    distopia::CalcBondsNoBoxIdx(this->coords, this->a_idx, this->b_idx, this->nidx, this->results);
+    distopia::CalcBondsNoBoxIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
         EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
@@ -529,10 +529,11 @@ TYPED_TEST(CoordinatesIdx, CalcBondsOrthoIdx) {
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
     // reference result
-    distopia::CalcBondsOrtho(this->a_coords_contig, this->b_coords_contig, this->nidx, this->box, this->ref_results);
+    distopia::CalcBondsOrtho<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->nidx, this->box, this->ref_results);
+
 
     // test the idx
-    distopia::CalcBondsOrthoIdx(this->coords, this->a_idx, this->b_idx, this->nidx, this->box, this->results);
+    distopia::CalcBondsOrthoIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->box, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
         EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
@@ -540,6 +541,7 @@ TYPED_TEST(CoordinatesIdx, CalcBondsOrthoIdx) {
 }
 
 
+
 TYPED_TEST(CoordinatesIdx, CalcBondsTriclinicIdx) {
     int ncoords = 250;
     int nidx = 100;
@@ -547,10 +549,10 @@ TYPED_TEST(CoordinatesIdx, CalcBondsTriclinicIdx) {
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
     // reference result
-    distopia::CalcBondsTriclinic(this->a_coords_contig, this->b_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
+    distopia::CalcBondsTriclinic<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
 
     // test the idx
-    distopia::CalcBondsTriclinicIdx(this->coords, this->a_idx, this->b_idx, this->nidx, this->triclinic_box, this->results);
+    distopia::CalcBondsTriclinicIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->triclinic_box, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
         EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
@@ -565,10 +567,10 @@ TYPED_TEST(CoordinatesIdx, CalcAnglesNoBoxIdx) {
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
     // reference result
-    distopia::CalcAnglesNoBox(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->nidx, this->ref_results);
+    distopia::CalcAnglesNoBox<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->nidx, this->ref_results);
 
     // test the idx
-    distopia::CalcAnglesNoBoxIdx(this->coords, this->a_idx, this->b_idx, this->c_idx, this->nidx, this->results);
+    distopia::CalcAnglesNoBoxIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->nidx, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
         EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
@@ -583,10 +585,10 @@ TYPED_TEST(CoordinatesIdx, CalcAnglesOrthoIdx) {
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
     // reference result
-    distopia::CalcAnglesOrtho(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->nidx, this->box, this->ref_results);
+    distopia::CalcAnglesOrtho<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->nidx, this->box, this->ref_results);
 
     // test the idx
-    distopia::CalcAnglesOrthoIdx(this->coords, this->a_idx, this->b_idx, this->c_idx, this->nidx, this->box, this->results);
+    distopia::CalcAnglesOrthoIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->nidx, this->box, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
         EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
@@ -601,10 +603,10 @@ TYPED_TEST(CoordinatesIdx, CalcAnglesTriclinicIdx) {
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
     // reference result
-    distopia::CalcAnglesTriclinic(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
+    distopia::CalcAnglesTriclinic<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
 
     // test the idx
-    distopia::CalcAnglesTriclinicIdx(this->coords, this->a_idx, this->b_idx, this->c_idx, this->nidx, this->triclinic_box, this->results);
+    distopia::CalcAnglesTriclinicIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->nidx, this->triclinic_box, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
         EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
@@ -619,10 +621,10 @@ TYPED_TEST(CoordinatesIdx, CalcDihedralsNoBoxIdx) {
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
     // reference result
-    distopia::CalcDihedralsNoBox(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->d_coords_contig, this->nidx, this->ref_results);
+    distopia::CalcDihedralsNoBox<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->d_coords_contig, this->nidx, this->ref_results);
 
     // test the idx
-    distopia::CalcDihedralsNoBoxIdx(this->coords, this->a_idx, this->b_idx, this->c_idx, this->d_idx, this->nidx, this->results);
+    distopia::CalcDihedralsNoBoxIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->d_idx, this->nidx, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
         EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
@@ -637,10 +639,10 @@ TYPED_TEST(CoordinatesIdx, CalcDihedralsOrthoIdx) {
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
     // reference result
-    distopia::CalcDihedralsOrtho(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->d_coords_contig, this->nidx, this->box, this->ref_results);
+    distopia::CalcDihedralsOrtho<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->d_coords_contig, this->nidx, this->box, this->ref_results);
 
     // test the idx
-    distopia::CalcDihedralsOrthoIdx(this->coords, this->a_idx, this->b_idx, this->c_idx, this->d_idx, this->nidx, this->box, this->results);
+    distopia::CalcDihedralsOrthoIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->d_idx, this->nidx, this->box, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
         EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
@@ -655,10 +657,10 @@ TYPED_TEST(CoordinatesIdx, CalcDihedralsTriclinicIdx) {
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
     // reference result
-    distopia::CalcDihedralsTriclinic(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->d_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
+    distopia::CalcDihedralsTriclinic<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->d_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
 
     // test the idx
-    distopia::CalcDihedralsTriclinicIdx(this->coords, this->a_idx, this->b_idx, this->c_idx, this->d_idx, this->nidx, this->triclinic_box, this->results);
+    distopia::CalcDihedralsTriclinicIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->d_idx, this->nidx, this->triclinic_box, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
         EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
@@ -673,10 +675,10 @@ TYPED_TEST(CoordinatesIdx, DistanceArrayNoBoxIdx) {
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
     // ref
-    distopia::CalcDistanceArrayNoBox(this->a_coords_contig, this->b_coords_contig, this->nidx, this->nidx, this->ref_results);
+    distopia::CalcDistanceArrayNoBox<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->nidx, this->nidx, this->ref_results);
 
     // test
-    distopia::CalcDistanceArrayNoBoxIdx(this->coords, this->a_idx, this->b_idx, this->nidx, this->nidx, this->results);
+    distopia::CalcDistanceArrayNoBoxIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->nidx, this->results);
 
     size_t n_results = nidx * nidx;
     for (std::size_t i=0; i<n_results; i++) {
@@ -693,10 +695,10 @@ TYPED_TEST(CoordinatesIdx, DistanceArrayOrthoIdx) {
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
     // ref
-    distopia::CalcDistanceArrayOrtho(this->a_coords_contig, this->b_coords_contig, this->nidx, this->nidx, this->box, this->ref_results);
+    distopia::CalcDistanceArrayOrtho<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->nidx, this->nidx, this->box, this->ref_results);
 
     // test
-    distopia::CalcDistanceArrayOrthoIdx(this->coords, this->a_idx, this->b_idx, this->nidx, this->nidx, this->box, this->results);
+    distopia::CalcDistanceArrayOrthoIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->nidx, this->box, this->results);
 
     size_t n_results = nidx * nidx;
     for (std::size_t i=0; i<n_results; i++) {
@@ -713,10 +715,10 @@ TYPED_TEST(CoordinatesIdx, DistanceArrayTriclinicIdx) {
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
     // ref
-    distopia::CalcDistanceArrayTriclinic(this->a_coords_contig, this->b_coords_contig, this->nidx, this->nidx, this->triclinic_box, this->ref_results);
+    distopia::CalcDistanceArrayTriclinic<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->nidx, this->nidx, this->triclinic_box, this->ref_results);
 
     // test
-    distopia::CalcDistanceArrayTriclinicIdx(this->coords, this->a_idx, this->b_idx, this->nidx, this->nidx, this->triclinic_box, this->results);
+    distopia::CalcDistanceArrayTriclinicIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->nidx, this->triclinic_box, this->results);
 
     size_t n_results = nidx * nidx;
     for (std::size_t i=0; i<n_results; i++) {
@@ -733,10 +735,10 @@ TYPED_TEST(CoordinatesIdx, SelfDistanceArrayNoBoxIdx) {
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
     // ref
-    distopia::CalcSelfDistanceArrayNoBox(this->a_coords_contig, this->nidx, this->ref_results);
+    distopia::CalcSelfDistanceArrayNoBox<TypeParam>(this->a_coords_contig, this->nidx, this->ref_results);
 
     // test
-    distopia::CalcSelfDistanceArrayNoBoxIdx(this->coords, this->a_idx, this->nidx, this->results);
+    distopia::CalcSelfDistanceArrayNoBoxIdx<TypeParam>(this->coords, this->a_idx, this->nidx, this->results);
 
 
     size_t n_results = nidx * (nidx-1) / 2;
@@ -753,10 +755,10 @@ TYPED_TEST(CoordinatesIdx, SelfDistanceArrayOrthoIdx) {
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
     // ref
-    distopia::CalcSelfDistanceArrayOrtho(this->a_coords_contig, this->nidx, this->box, this->ref_results);
+    distopia::CalcSelfDistanceArrayOrtho<TypeParam>(this->a_coords_contig, this->nidx, this->box, this->ref_results);
 
     // test
-    distopia::CalcSelfDistanceArrayOrthoIdx(this->coords, this->a_idx, this->nidx, this->box, this->results);
+    distopia::CalcSelfDistanceArrayOrthoIdx<TypeParam>(this->coords, this->a_idx, this->nidx, this->box, this->results);
 
 
     size_t n_results = nidx * (nidx-1) / 2;
@@ -773,10 +775,10 @@ TYPED_TEST(CoordinatesIdx, SelfDistanceArrayTriclinicIdx) {
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
     // ref
-    distopia::CalcSelfDistanceArrayTriclinic(this->a_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
+    distopia::CalcSelfDistanceArrayTriclinic<TypeParam>(this->a_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
 
     // test
-    distopia::CalcSelfDistanceArrayTriclinicIdx(this->coords, this->a_idx, this->nidx, this->triclinic_box, this->results);
+    distopia::CalcSelfDistanceArrayTriclinicIdx<TypeParam>(this->coords, this->a_idx, this->nidx, this->triclinic_box, this->results);
 
 
     size_t n_results = nidx * (nidx-1) / 2;
diff --git a/libdistopia/test/test_utils.h b/libdistopia/test/test_utils.h
index 469da52f..900f745a 100644
--- a/libdistopia/test/test_utils.h
+++ b/libdistopia/test/test_utils.h
@@ -4,6 +4,8 @@
 #include <random>
 
 #include <gtest/gtest.h>
+#include <gmock/gmock.h>
+
 // creates nrandom floating points between pos and neg limit
 template <typename T>
 void RandomFloatingPoint(T *target, const int nrandom, const int neglimit,

From 632102017787e565fbabdf265075c9b3da77f5bb Mon Sep 17 00:00:00 2001
From: hmacdope <hugomacdermott@gmail.com>
Date: Sat, 14 Sep 2024 10:40:59 +1000
Subject: [PATCH 08/12] Nan sensitive testing

---
 libdistopia/test/test_utils.h | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/libdistopia/test/test_utils.h b/libdistopia/test/test_utils.h
index 900f745a..5ba61063 100644
--- a/libdistopia/test/test_utils.h
+++ b/libdistopia/test/test_utils.h
@@ -35,21 +35,21 @@ void RandomInt(std::size_t *target, const int nrandom, const int neglimit,
 
 inline void EXPECT_SCALAR_EQ(float result, float ref)
 {
-    EXPECT_FLOAT_EQ(result, ref);
+    EXPECT_THAT(result, ::testing::NanSensitiveFloatEq(ref));
 }
 
 inline void EXPECT_SCALAR_EQ(double result, double ref)
 {
-    EXPECT_DOUBLE_EQ(result, ref);
+    EXPECT_THAT(result, ::testing::NanSensitiveDoubleEq(ref));
 }
 
 inline void EXPECT_SCALAR_NEAR(float result, float ref, float tol)
 {
-    EXPECT_NEAR(result, ref, tol);
+    EXPECT_THAT(result, ::testing::NanSensitiveFloatNear(ref, tol));
 }
 inline void EXPECT_SCALAR_NEAR(double result, double ref, float tol)
 {
-    EXPECT_NEAR(result, ref, tol);
+    EXPECT_THAT(result, ::testing::NanSensitiveDoubleNear(ref, tol));
 }
 
 

From 7b9c2a300ca9eb211f21b38813ad5c6e6a4800bd Mon Sep 17 00:00:00 2001
From: hmacdope <hugomacdermott@gmail.com>
Date: Sat, 14 Sep 2024 11:13:47 +1000
Subject: [PATCH 09/12] fix comp

---
 libdistopia/test/test_mda_match.cpp | 76 ++++++++++++++---------------
 libdistopia/test/test_utils.h       |  7 +++
 2 files changed, 45 insertions(+), 38 deletions(-)

diff --git a/libdistopia/test/test_mda_match.cpp b/libdistopia/test/test_mda_match.cpp
index 6c95ee03..bcba6423 100644
--- a/libdistopia/test/test_mda_match.cpp
+++ b/libdistopia/test/test_mda_match.cpp
@@ -40,7 +40,7 @@ TYPED_TEST(CoordinatesTest, CalcBondsOrthoMatchesMDA)
 
     for (std::size_t i = 0; i < this->nresults; i++)
     {
-        EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
     }
 }
 
@@ -57,7 +57,7 @@ TYPED_TEST(CoordinatesTest, CalcBondsNoBoxMatchesMDA)
 
   for (std::size_t i = 0; i < this->nresults; i++)
   {
-    EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+    EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
   }
 }
 
@@ -93,7 +93,7 @@ TYPED_TEST(CoordinatesTest, CalcBondsTriclinicMatchesMDA)
 
   for (std::size_t i = 0; i < this->nresults; i++)
   {
-    EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+    EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
   }
 }
 
@@ -112,7 +112,7 @@ TYPED_TEST(CoordinatesTest, CalcAnglesOrthoMatchesMDA)
 
   for (std::size_t i = 0; i < this->nresults; i++)
   {
-    EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+    EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
   }
 }
 
@@ -129,7 +129,7 @@ TYPED_TEST(CoordinatesTest, CalcAnglesNoBoxMatchesMDA)
 
   for (std::size_t i = 0; i < this->nresults; i++)
   {
-    EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+    EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
   }
 }
 
@@ -158,7 +158,7 @@ TYPED_TEST(CoordinatesTest, CalcAnglesTriclinicMatchesMDA)
 
     for (std::size_t i = 0; i < this->nresults; i++)
     {
-        EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
     }
 }
 
@@ -175,7 +175,7 @@ TYPED_TEST(CoordinatesTest, CalcDihedralsOrthoMatchesMDA)
 
   for (std::size_t i = 0; i < this->nresults; i++)
   {
-    EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+    EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
   }
 }
 
@@ -192,7 +192,7 @@ TYPED_TEST(CoordinatesTest, CalcDihedralsNoBoxMatchesMDA)
 
   for (std::size_t i = 0; i < this->nresults; i++)
   {
-    EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+    EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
   }
 }
 
@@ -222,7 +222,7 @@ TYPED_TEST(CoordinatesTest, CalcDihedralsTriclinicMatchesMDA)
 
     for (std::size_t i = 0; i < this->nresults; i++)
     {
-        EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
     }
 }
 
@@ -242,7 +242,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayOrthoMatchesMDA) {
 
     for (std::size_t i = 0; i < this->nresults; i++)
     {
-        EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
     }
 
 }
@@ -262,7 +262,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayOrthoMatchesMDAWide) {
 
     for (std::size_t i = 0; i < this->nresults; i++)
     {
-        EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
     }
 
 }
@@ -283,7 +283,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayOrthoMatchesMDAScalarPath)
 
     for (std::size_t i = 0; i < this->nresults; i++)
     {
-        EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
     }
 }
 
@@ -301,7 +301,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayNoBoxMatchesMDA) {
 
     for (std::size_t i = 0; i < this->nresults; i++)
     {
-        EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
     }
 }
 
@@ -321,7 +321,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayNoBoxMatchesMDAScalarPath)
 
     for (std::size_t i = 0; i <  this->nresults; i++)
     {
-        EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
     }
 }
 
@@ -339,7 +339,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayTriclinicMatchesMDA) {
 
     for (std::size_t i = 0;  i < this->nresults; i++)
     {
-        EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
     }
 }
 
@@ -359,7 +359,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayTriclinicMatchesMDAScalarP
 
     for (std::size_t i = 0; i < this->nresults; i++)
     {
-        EXPECT_NEAR(this->results[i], this->ref[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->results[i], this->ref[i], abs_err);
     }
 }
 
@@ -378,7 +378,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcSelfDistanceArrayNoBox) {
     //nresults >>= 2;
 
     for (std::size_t i=0; i<nresults; ++i) {
-        EXPECT_NEAR(this->ref[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref[i], this->results[i], abs_err);
     }
 }
 
@@ -398,7 +398,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcSelfDistanceArrayOrtho) {
     //nresults >>= 2;
 
     for (std::size_t i=0; i<nresults; ++i) {
-        EXPECT_NEAR(this->ref[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref[i], this->results[i], abs_err);
     }
 }
 
@@ -418,7 +418,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcSelfDistanceArrayTriclinic) {
     //nresults >>= 2;
 
     for (std::size_t i=0; i<nresults; ++i) {
-        EXPECT_NEAR(this->ref[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref[i], this->results[i], abs_err);
     }
 }
 
@@ -438,7 +438,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcSelfDistanceArrayNoBoxScalar) {
     //nresults >>= 2;
 
     for (std::size_t i=0; i<nresults; ++i) {
-        EXPECT_NEAR(this->ref[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref[i], this->results[i], abs_err);
     }
 }
 
@@ -458,7 +458,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcSelfDistanceArrayOrthoScalar) {
     //nresults >>= 2;
 
     for (std::size_t i=0; i<nresults; ++i) {
-        EXPECT_NEAR(this->ref[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref[i], this->results[i], abs_err);
     }
 }
 
@@ -478,7 +478,7 @@ TYPED_TEST(DistanceArrayCoordinates, CalcSelfDistanceArrayTriclinicScalar) {
     //nresults >>= 2;
 
     for (std::size_t i=0; i<nresults; ++i) {
-        EXPECT_NEAR(this->ref[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref[i], this->results[i], abs_err);
     }
 }
 
@@ -499,7 +499,7 @@ TYPED_TEST(CoordinatesIdx, CalcBondsNoBoxIdx) {
     distopia::CalcBondsNoBoxIdx(this->coords, this->a_idx, this->b_idx, this->nidx, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
@@ -517,7 +517,7 @@ TYPED_TEST(CoordinatesIdx, CalcBondsNoBoxIdxSmall) {
     distopia::CalcBondsNoBoxIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
@@ -536,7 +536,7 @@ TYPED_TEST(CoordinatesIdx, CalcBondsOrthoIdx) {
     distopia::CalcBondsOrthoIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->box, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
@@ -555,7 +555,7 @@ TYPED_TEST(CoordinatesIdx, CalcBondsTriclinicIdx) {
     distopia::CalcBondsTriclinicIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->triclinic_box, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
@@ -573,7 +573,7 @@ TYPED_TEST(CoordinatesIdx, CalcAnglesNoBoxIdx) {
     distopia::CalcAnglesNoBoxIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->nidx, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
@@ -591,7 +591,7 @@ TYPED_TEST(CoordinatesIdx, CalcAnglesOrthoIdx) {
     distopia::CalcAnglesOrthoIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->nidx, this->box, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
@@ -609,7 +609,7 @@ TYPED_TEST(CoordinatesIdx, CalcAnglesTriclinicIdx) {
     distopia::CalcAnglesTriclinicIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->nidx, this->triclinic_box, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
@@ -627,7 +627,7 @@ TYPED_TEST(CoordinatesIdx, CalcDihedralsNoBoxIdx) {
     distopia::CalcDihedralsNoBoxIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->d_idx, this->nidx, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
@@ -645,7 +645,7 @@ TYPED_TEST(CoordinatesIdx, CalcDihedralsOrthoIdx) {
     distopia::CalcDihedralsOrthoIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->d_idx, this->nidx, this->box, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
@@ -663,7 +663,7 @@ TYPED_TEST(CoordinatesIdx, CalcDihedralsTriclinicIdx) {
     distopia::CalcDihedralsTriclinicIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->d_idx, this->nidx, this->triclinic_box, this->results);
 
     for (std::size_t i=0; i<this->nidx; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
@@ -682,7 +682,7 @@ TYPED_TEST(CoordinatesIdx, DistanceArrayNoBoxIdx) {
 
     size_t n_results = nidx * nidx;
     for (std::size_t i=0; i<n_results; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 
 }
@@ -702,7 +702,7 @@ TYPED_TEST(CoordinatesIdx, DistanceArrayOrthoIdx) {
 
     size_t n_results = nidx * nidx;
     for (std::size_t i=0; i<n_results; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 
 }
@@ -722,7 +722,7 @@ TYPED_TEST(CoordinatesIdx, DistanceArrayTriclinicIdx) {
 
     size_t n_results = nidx * nidx;
     for (std::size_t i=0; i<n_results; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 
 }
@@ -743,7 +743,7 @@ TYPED_TEST(CoordinatesIdx, SelfDistanceArrayNoBoxIdx) {
 
     size_t n_results = nidx * (nidx-1) / 2;
     for (std::size_t i=0; i<n_results; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
@@ -763,7 +763,7 @@ TYPED_TEST(CoordinatesIdx, SelfDistanceArrayOrthoIdx) {
 
     size_t n_results = nidx * (nidx-1) / 2;
     for (std::size_t i=0; i<n_results; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
 
@@ -783,6 +783,6 @@ TYPED_TEST(CoordinatesIdx, SelfDistanceArrayTriclinicIdx) {
 
     size_t n_results = nidx * (nidx-1) / 2;
     for (std::size_t i=0; i<n_results; i++) {
-        EXPECT_NEAR(this->ref_results[i], this->results[i], abs_err);
+        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 }
\ No newline at end of file
diff --git a/libdistopia/test/test_utils.h b/libdistopia/test/test_utils.h
index 5ba61063..842432b6 100644
--- a/libdistopia/test/test_utils.h
+++ b/libdistopia/test/test_utils.h
@@ -6,6 +6,7 @@
 #include <gtest/gtest.h>
 #include <gmock/gmock.h>
 
+
 // creates nrandom floating points between pos and neg limit
 template <typename T>
 void RandomFloatingPoint(T *target, const int nrandom, const int neglimit,
@@ -35,20 +36,24 @@ void RandomInt(std::size_t *target, const int nrandom, const int neglimit,
 
 inline void EXPECT_SCALAR_EQ(float result, float ref)
 {
+
     EXPECT_THAT(result, ::testing::NanSensitiveFloatEq(ref));
 }
 
 inline void EXPECT_SCALAR_EQ(double result, double ref)
 {
+
     EXPECT_THAT(result, ::testing::NanSensitiveDoubleEq(ref));
 }
 
 inline void EXPECT_SCALAR_NEAR(float result, float ref, float tol)
 {
+
     EXPECT_THAT(result, ::testing::NanSensitiveFloatNear(ref, tol));
 }
 inline void EXPECT_SCALAR_NEAR(double result, double ref, float tol)
 {
+
     EXPECT_THAT(result, ::testing::NanSensitiveDoubleNear(ref, tol));
 }
 
@@ -72,4 +77,6 @@ void pretty_print_matrix(T *matrix, int rows, int cols)
 
 
 
+
+
 #endif // DISTOPIA_TEST_UTILS_H
\ No newline at end of file

From 55e68d12120825599d5e45ba0ac1beab1ec546f8 Mon Sep 17 00:00:00 2001
From: hmacdope <hugomacdermott@gmail.com>
Date: Sat, 14 Sep 2024 11:26:02 +1000
Subject: [PATCH 10/12] remove prints

---
 libdistopia/src/distopia.cpp | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/libdistopia/src/distopia.cpp b/libdistopia/src/distopia.cpp
index dfb5e6fc..7bff6f1a 100644
--- a/libdistopia/src/distopia.cpp
+++ b/libdistopia/src/distopia.cpp
@@ -759,10 +759,6 @@ namespace distopia {
             auto b_y = hn::Undefined(d);
             auto b_z = hn::Undefined(d);
 
-            std::cout << "nlanes: " << nlanes << std::endl;
-            std::cout << "n: " << n << std::endl;
-            std::cout << "n - nlanes: " << n - nlanes << std::endl;
-
             for (size_t i=0; i <= n - nlanes; i += nlanes) {
                 // load N indices of each source
                 // interleaved gather these indices
@@ -773,7 +769,6 @@ namespace distopia {
 
                 hn::StoreU(result, d, dst + i);
             }
-            std::cout << "n % nlanes: " << n % nlanes << std::endl;
             size_t rem = n % nlanes;
             if (rem) {  // if we had a non-multiple of nlanes, do final nlanes values again
                 LoadInterleavedIdx(a_idx + n - nlanes, coords, a_x, a_y, a_z);

From 2bb5b157457810aba2d3d18a0adc71fb0d10ab4e Mon Sep 17 00:00:00 2001
From: hmacdope <hugomacdermott@gmail.com>
Date: Sat, 14 Sep 2024 12:04:58 +1000
Subject: [PATCH 11/12] temp silence test noise

---
 libdistopia/test/test_mda_match.cpp | 203 ++++++++++++++--------------
 1 file changed, 101 insertions(+), 102 deletions(-)

diff --git a/libdistopia/test/test_mda_match.cpp b/libdistopia/test/test_mda_match.cpp
index bcba6423..3e0ada92 100644
--- a/libdistopia/test/test_mda_match.cpp
+++ b/libdistopia/test/test_mda_match.cpp
@@ -493,7 +493,10 @@ TYPED_TEST(CoordinatesIdx, CalcBondsNoBoxIdx) {
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
     // reference result
-    distopia::CalcBondsNoBox(this->a_coords_contig, this->b_coords_contig, this->nidx, this->ref_results);
+    // distopia::CalcBondsNoBox(this->a_coords_contig, this->b_coords_contig, this->nidx, this->ref_results);
+    using ctype = ScalarToCoordinateT<TypeParam>;
+
+    _calc_bond_distance((ctype*)this->a_coords_contig, (ctype*)this->b_coords_contig, this->nidx, this->ref_results);
 
     // test the idx
     distopia::CalcBondsNoBoxIdx(this->coords, this->a_idx, this->b_idx, this->nidx, this->results);
@@ -501,6 +504,7 @@ TYPED_TEST(CoordinatesIdx, CalcBondsNoBoxIdx) {
     for (std::size_t i=0; i<this->nidx; i++) {
         EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
+
 }
 
 TYPED_TEST(CoordinatesIdx, CalcBondsNoBoxIdxSmall) {
@@ -509,10 +513,10 @@ TYPED_TEST(CoordinatesIdx, CalcBondsNoBoxIdxSmall) {
     int nidx = 3;
 
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
+    using ctype = ScalarToCoordinateT<TypeParam>;
 
     // reference result
-    distopia::CalcBondsNoBox<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->nidx, this->ref_results);
-    std::cout << "ref_results: " << std::endl;
+    _calc_bond_distance((ctype*)this->a_coords_contig, (ctype*)this->b_coords_contig, this->nidx, this->ref_results);
     // test the idx
     distopia::CalcBondsNoBoxIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->results);
 
@@ -527,10 +531,10 @@ TYPED_TEST(CoordinatesIdx, CalcBondsOrthoIdx) {
     int nidx = 100;
 
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
+    using ctype = ScalarToCoordinateT<TypeParam>;
 
     // reference result
-    distopia::CalcBondsOrtho<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->nidx, this->box, this->ref_results);
-
+    _calc_bond_distance_ortho((ctype*)this->a_coords_contig, (ctype*)this->b_coords_contig, this->nidx, this->box, this->ref_results);
 
     // test the idx
     distopia::CalcBondsOrthoIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->box, this->results);
@@ -547,10 +551,10 @@ TYPED_TEST(CoordinatesIdx, CalcBondsTriclinicIdx) {
     int nidx = 100;
 
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
+    using ctype = ScalarToCoordinateT<TypeParam>;
 
     // reference result
-    distopia::CalcBondsTriclinic<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
-
+    _calc_bond_distance_triclinic((ctype*)this->a_coords_contig, (ctype*)this->b_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
     // test the idx
     distopia::CalcBondsTriclinicIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->triclinic_box, this->results);
 
@@ -565,10 +569,9 @@ TYPED_TEST(CoordinatesIdx, CalcAnglesNoBoxIdx) {
     int nidx = 100;
 
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
-
+    using ctype = ScalarToCoordinateT<TypeParam>;
     // reference result
-    distopia::CalcAnglesNoBox<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->nidx, this->ref_results);
-
+    _calc_angle((ctype*)this->a_coords_contig, (ctype*)this->b_coords_contig, (ctype*)this->c_coords_contig, this->nidx, this->ref_results);
     // test the idx
     distopia::CalcAnglesNoBoxIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->nidx, this->results);
 
@@ -583,10 +586,9 @@ TYPED_TEST(CoordinatesIdx, CalcAnglesOrthoIdx) {
     int nidx = 100;
 
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
-
+    using ctype = ScalarToCoordinateT<TypeParam>;
     // reference result
-    distopia::CalcAnglesOrtho<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->nidx, this->box, this->ref_results);
-
+    _calc_angle_ortho((ctype*)this->a_coords_contig, (ctype*)this->b_coords_contig, (ctype*)this->c_coords_contig, this->nidx, this->box, this->ref_results);
     // test the idx
     distopia::CalcAnglesOrthoIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->nidx, this->box, this->results);
 
@@ -599,11 +601,11 @@ TYPED_TEST(CoordinatesIdx, CalcAnglesOrthoIdx) {
 TYPED_TEST(CoordinatesIdx, CalcAnglesTriclinicIdx) {
     int ncoords = 250;
     int nidx = 100;
-
+    using ctype = ScalarToCoordinateT<TypeParam>;
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
     // reference result
-    distopia::CalcAnglesTriclinic<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
+    _calc_angle_triclinic((ctype*)this->a_coords_contig, (ctype*)this->b_coords_contig, (ctype*)this->c_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
 
     // test the idx
     distopia::CalcAnglesTriclinicIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->nidx, this->triclinic_box, this->results);
@@ -619,10 +621,9 @@ TYPED_TEST(CoordinatesIdx, CalcDihedralsNoBoxIdx) {
     int nidx = 100;
 
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
-
+    using ctype = ScalarToCoordinateT<TypeParam>;
     // reference result
-    distopia::CalcDihedralsNoBox<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->d_coords_contig, this->nidx, this->ref_results);
-
+    _calc_dihedral((ctype*)this->a_coords_contig, (ctype*)this->b_coords_contig, (ctype*)this->c_coords_contig, (ctype*)this->d_coords_contig, this->nidx, this->ref_results);
     // test the idx
     distopia::CalcDihedralsNoBoxIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->d_idx, this->nidx, this->results);
 
@@ -637,10 +638,9 @@ TYPED_TEST(CoordinatesIdx, CalcDihedralsOrthoIdx) {
     int nidx = 100;
 
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
-
+    using ctype = ScalarToCoordinateT<TypeParam>;
     // reference result
-    distopia::CalcDihedralsOrtho<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->d_coords_contig, this->nidx, this->box, this->ref_results);
-
+    _calc_dihedral_ortho((ctype*)this->a_coords_contig, (ctype*)this->b_coords_contig, (ctype*)this->c_coords_contig, (ctype*)this->d_coords_contig, this->nidx, this->box, this->ref_results);
     // test the idx
     distopia::CalcDihedralsOrthoIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->d_idx, this->nidx, this->box, this->results);
 
@@ -655,10 +655,9 @@ TYPED_TEST(CoordinatesIdx, CalcDihedralsTriclinicIdx) {
     int nidx = 100;
 
     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
-
+    using ctype = ScalarToCoordinateT<TypeParam>;
     // reference result
-    distopia::CalcDihedralsTriclinic<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->c_coords_contig, this->d_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
-
+    _calc_dihedral_triclinic((ctype*)this->a_coords_contig, (ctype*)this->b_coords_contig, (ctype*)this->c_coords_contig, (ctype*)this->d_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
     // test the idx
     distopia::CalcDihedralsTriclinicIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->c_idx, this->d_idx, this->nidx, this->triclinic_box, this->results);
 
@@ -668,121 +667,121 @@ TYPED_TEST(CoordinatesIdx, CalcDihedralsTriclinicIdx) {
 }
 
 
-TYPED_TEST(CoordinatesIdx, DistanceArrayNoBoxIdx) {
-    int ncoords = 250;
-    int nidx = 100;
+// TYPED_TEST(CoordinatesIdx, DistanceArrayNoBoxIdx) {
+//     int ncoords = 250;
+//     int nidx = 100;
 
-    this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
+//     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
-    // ref
-    distopia::CalcDistanceArrayNoBox<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->nidx, this->nidx, this->ref_results);
+//     // ref
+//     distopia::CalcDistanceArrayNoBox<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->nidx, this->nidx, this->ref_results);
 
-    // test
-    distopia::CalcDistanceArrayNoBoxIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->nidx, this->results);
+//     // test
+//     distopia::CalcDistanceArrayNoBoxIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->nidx, this->results);
 
-    size_t n_results = nidx * nidx;
-    for (std::size_t i=0; i<n_results; i++) {
-        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
-    }
+//     size_t n_results = nidx * nidx;
+//     for (std::size_t i=0; i<n_results; i++) {
+//         EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
+//     }
 
-}
+// }
 
 
-TYPED_TEST(CoordinatesIdx, DistanceArrayOrthoIdx) {
-    int ncoords = 250;
-    int nidx = 100;
+// TYPED_TEST(CoordinatesIdx, DistanceArrayOrthoIdx) {
+//     int ncoords = 250;
+//     int nidx = 100;
 
-    this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
+//     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
-    // ref
-    distopia::CalcDistanceArrayOrtho<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->nidx, this->nidx, this->box, this->ref_results);
+//     // ref
+//     distopia::CalcDistanceArrayOrtho<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->nidx, this->nidx, this->box, this->ref_results);
 
-    // test
-    distopia::CalcDistanceArrayOrthoIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->nidx, this->box, this->results);
+//     // test
+//     distopia::CalcDistanceArrayOrthoIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->nidx, this->box, this->results);
 
-    size_t n_results = nidx * nidx;
-    for (std::size_t i=0; i<n_results; i++) {
-        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
-    }
+//     size_t n_results = nidx * nidx;
+//     for (std::size_t i=0; i<n_results; i++) {
+//         EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
+//     }
 
-}
+// }
 
 
-TYPED_TEST(CoordinatesIdx, DistanceArrayTriclinicIdx) {
-    int ncoords = 250;
-    int nidx = 100;
+// TYPED_TEST(CoordinatesIdx, DistanceArrayTriclinicIdx) {
+//     int ncoords = 250;
+//     int nidx = 100;
 
-    this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
+//     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
-    // ref
-    distopia::CalcDistanceArrayTriclinic<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->nidx, this->nidx, this->triclinic_box, this->ref_results);
+//     // ref
+//     distopia::CalcDistanceArrayTriclinic<TypeParam>(this->a_coords_contig, this->b_coords_contig, this->nidx, this->nidx, this->triclinic_box, this->ref_results);
 
-    // test
-    distopia::CalcDistanceArrayTriclinicIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->nidx, this->triclinic_box, this->results);
+//     // test
+//     distopia::CalcDistanceArrayTriclinicIdx<TypeParam>(this->coords, this->a_idx, this->b_idx, this->nidx, this->nidx, this->triclinic_box, this->results);
 
-    size_t n_results = nidx * nidx;
-    for (std::size_t i=0; i<n_results; i++) {
-        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
-    }
+//     size_t n_results = nidx * nidx;
+//     for (std::size_t i=0; i<n_results; i++) {
+//         EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
+//     }
 
-}
+// }
 
 
-TYPED_TEST(CoordinatesIdx, SelfDistanceArrayNoBoxIdx) {
-    int ncoords = 250;
-    int nidx = 100;
+// TYPED_TEST(CoordinatesIdx, SelfDistanceArrayNoBoxIdx) {
+//     int ncoords = 250;
+//     int nidx = 100;
 
-    this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
+//     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
-    // ref
-    distopia::CalcSelfDistanceArrayNoBox<TypeParam>(this->a_coords_contig, this->nidx, this->ref_results);
+//     // ref
+//     distopia::CalcSelfDistanceArrayNoBox<TypeParam>(this->a_coords_contig, this->nidx, this->ref_results);
 
-    // test
-    distopia::CalcSelfDistanceArrayNoBoxIdx<TypeParam>(this->coords, this->a_idx, this->nidx, this->results);
+//     // test
+//     distopia::CalcSelfDistanceArrayNoBoxIdx<TypeParam>(this->coords, this->a_idx, this->nidx, this->results);
 
 
-    size_t n_results = nidx * (nidx-1) / 2;
-    for (std::size_t i=0; i<n_results; i++) {
-        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
-    }
-}
+//     size_t n_results = nidx * (nidx-1) / 2;
+//     for (std::size_t i=0; i<n_results; i++) {
+//         EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
+//     }
+// }
 
 
-TYPED_TEST(CoordinatesIdx, SelfDistanceArrayOrthoIdx) {
-    int ncoords = 250;
-    int nidx = 100;
+// TYPED_TEST(CoordinatesIdx, SelfDistanceArrayOrthoIdx) {
+//     int ncoords = 250;
+//     int nidx = 100;
 
-    this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
+//     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
-    // ref
-    distopia::CalcSelfDistanceArrayOrtho<TypeParam>(this->a_coords_contig, this->nidx, this->box, this->ref_results);
+//     // ref
+//     distopia::CalcSelfDistanceArrayOrtho<TypeParam>(this->a_coords_contig, this->nidx, this->box, this->ref_results);
 
-    // test
-    distopia::CalcSelfDistanceArrayOrthoIdx<TypeParam>(this->coords, this->a_idx, this->nidx, this->box, this->results);
+//     // test
+//     distopia::CalcSelfDistanceArrayOrthoIdx<TypeParam>(this->coords, this->a_idx, this->nidx, this->box, this->results);
 
 
-    size_t n_results = nidx * (nidx-1) / 2;
-    for (std::size_t i=0; i<n_results; i++) {
-        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
-    }
-}
+//     size_t n_results = nidx * (nidx-1) / 2;
+//     for (std::size_t i=0; i<n_results; i++) {
+//         EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
+//     }
+// }
 
 
-TYPED_TEST(CoordinatesIdx, SelfDistanceArrayTriclinicIdx) {
-    int ncoords = 250;
-    int nidx = 100;
+// TYPED_TEST(CoordinatesIdx, SelfDistanceArrayTriclinicIdx) {
+//     int ncoords = 250;
+//     int nidx = 100;
 
-    this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
+//     this->SetUp(ncoords, nidx, BOXSIZE, 3 * BOXSIZE);
 
-    // ref
-    distopia::CalcSelfDistanceArrayTriclinic<TypeParam>(this->a_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
+//     // ref
+//     distopia::CalcSelfDistanceArrayTriclinic<TypeParam>(this->a_coords_contig, this->nidx, this->triclinic_box, this->ref_results);
 
-    // test
-    distopia::CalcSelfDistanceArrayTriclinicIdx<TypeParam>(this->coords, this->a_idx, this->nidx, this->triclinic_box, this->results);
+//     // test
+//     distopia::CalcSelfDistanceArrayTriclinicIdx<TypeParam>(this->coords, this->a_idx, this->nidx, this->triclinic_box, this->results);
 
 
-    size_t n_results = nidx * (nidx-1) / 2;
-    for (std::size_t i=0; i<n_results; i++) {
-        EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
-    }
-}
\ No newline at end of file
+//     size_t n_results = nidx * (nidx-1) / 2;
+//     for (std::size_t i=0; i<n_results; i++) {
+//         EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
+//     }
+// }
\ No newline at end of file

From 24293e6bf8e08b693681787863046fa4b324b00f Mon Sep 17 00:00:00 2001
From: hmacdope <hugomacdermott@gmail.com>
Date: Sun, 13 Oct 2024 11:35:16 +1100
Subject: [PATCH 12/12] wip?

---
 libdistopia/src/distopia.cpp              | 4 ++--
 libdistopia/test/compare/calc_distances.h | 3 +++
 libdistopia/test/test_mda_match.cpp       | 9 +++++++--
 libdistopia/test/test_utils.h             | 8 ++++++++
 4 files changed, 20 insertions(+), 4 deletions(-)

diff --git a/libdistopia/src/distopia.cpp b/libdistopia/src/distopia.cpp
index 7bff6f1a..212bfb22 100644
--- a/libdistopia/src/distopia.cpp
+++ b/libdistopia/src/distopia.cpp
@@ -54,7 +54,7 @@ namespace distopia {
                 acc = hn::MulAdd(dx, dx, acc);
                 acc = hn::MulAdd(dy, dy, acc);
                 acc = hn::MulAdd(dz, dz, acc);
-
+                hn::Print(d, "rsq", acc);
                 return hn::Sqrt(acc);
             }
         };
@@ -492,7 +492,7 @@ namespace distopia {
 
             y = y / vb_norm;
 
-            return   hn::Neg(hn::Atan2(d, y, x));
+            return hn::Neg(hn::Atan2(d, y, x));
         }
 
         template <typename T, typename B>
diff --git a/libdistopia/test/compare/calc_distances.h b/libdistopia/test/compare/calc_distances.h
index 4be27d57..330b0e30 100644
--- a/libdistopia/test/compare/calc_distances.h
+++ b/libdistopia/test/compare/calc_distances.h
@@ -607,7 +607,10 @@ static void _calc_bond_distance(ScalarToCoordinateT<T>* atom1, ScalarToCoordinat
     dx[1] = atom1[i][1] - atom2[i][1];
     dx[2] = atom1[i][2] - atom2[i][2];
     T rsq = (dx[0]*dx[0])+(dx[1]*dx[1])+(dx[2]*dx[2]);
+    std::cout << "rsq = " << rsq << std::endl;
+
     *(distances+i) = sqrt(rsq);
+    std::cout << "distances[" << i << "] = " << sqrt(rsq) << std::endl;
   }
 }
 
diff --git a/libdistopia/test/test_mda_match.cpp b/libdistopia/test/test_mda_match.cpp
index 3e0ada92..c71cd234 100644
--- a/libdistopia/test/test_mda_match.cpp
+++ b/libdistopia/test/test_mda_match.cpp
@@ -15,8 +15,8 @@ typedef Types<float, double> ScalarTypes;
 
 // constants
 constexpr int BOXSIZE = 30;
-constexpr int NRESULTS = 10;
-constexpr int NINDICIES = 5;
+constexpr int NRESULTS = 250;
+constexpr int NINDICIES = 100;
 constexpr double abs_err = 1.0e-4;
 
 
@@ -505,6 +505,11 @@ TYPED_TEST(CoordinatesIdx, CalcBondsNoBoxIdx) {
         EXPECT_SCALAR_NEAR(this->ref_results[i], this->results[i], abs_err);
     }
 
+    // print 
+    for (std::size_t i=0; i<this->nidx; i++) {
+        std::cout << "i: " << i << " res: " << this->results[i] << " ressq: " <<  std::pow(this->results[i], 2) << " ref: " << this->ref_results[i] << " refsq: " << std::pow(this->ref_results[i], 2) << " near: " << isnear(this->results[i], this->ref_results[i], abs_err) << std::endl;
+    }
+
 }
 
 TYPED_TEST(CoordinatesIdx, CalcBondsNoBoxIdxSmall) {
diff --git a/libdistopia/test/test_utils.h b/libdistopia/test/test_utils.h
index 842432b6..2057686b 100644
--- a/libdistopia/test/test_utils.h
+++ b/libdistopia/test/test_utils.h
@@ -58,6 +58,14 @@ inline void EXPECT_SCALAR_NEAR(double result, double ref, float tol)
 }
 
 
+// isnear 
+
+bool isnear(float a, float b, float tol)
+{
+    return std::abs(a - b) < tol;
+}
+
+
 template <typename T>
 void pretty_print_matrix(T *matrix, int rows, int cols)
 {