From beeb9931de0db80baca803ae374d270c6c1e6ba6 Mon Sep 17 00:00:00 2001 From: Pardhav Maradani Date: Wed, 18 Feb 2026 12:06:54 +0530 Subject: [PATCH 1/5] Add support for MDA minimize_vectors --- distopia/__init__.py | 4 + distopia/pydistopia.pyx | 129 +++++++++++++++++++++++++++++++- distopia/tests/test_distopia.py | 73 ++++++++++++++++++ docs/docs/api.md | 8 +- docs/docs/examples.md | 35 +++++++++ libdistopia/include/distopia.h | 2 + libdistopia/src/distopia.cpp | 74 ++++++++++++++++++ 7 files changed, 323 insertions(+), 2 deletions(-) diff --git a/distopia/__init__.py b/distopia/__init__.py index 5740293..742e5e1 100644 --- a/distopia/__init__.py +++ b/distopia/__init__.py @@ -15,6 +15,8 @@ self_distance_array_no_box, self_distance_array_ortho, self_distance_array_triclinic, + minimize_vectors_ortho, + minimize_vectors_triclinic, ) __all__ = [ @@ -33,6 +35,8 @@ 'self_distance_array_no_box', 'self_distance_array_ortho', 'self_distance_array_triclinic', + 'minimize_vectors_ortho', + 'minimize_vectors_triclinic', ] from importlib.metadata import version diff --git a/distopia/pydistopia.pyx b/distopia/pydistopia.pyx index e62f2cb..68de81a 100644 --- a/distopia/pydistopia.pyx +++ b/distopia/pydistopia.pyx @@ -129,6 +129,19 @@ cdef extern from "distopia.h" namespace "distopia" nogil: const T *box, T *out, ) + void MinimizeVectorsOrtho[T]( + const T *vectors, + size_t n, + const T *box, + T *output, + ) + void MinimizeVectorsTriclinic[T]( + const T *vectors, + size_t n, + const T *box, + T *output, + ) + def get_n_float_lanes(): """The number of floats per register distopia will handle on this system""" @@ -902,4 +915,118 @@ def self_distance_array_triclinic( &box[0][0], &results_view[0]) - return np.array(results) \ No newline at end of file + return np.array(results) + + +@cython.boundscheck(False) +@cython.wraparound(False) +def minimize_vectors_ortho( + floating[:, ::1] vectors, + floating[::1] box, + floating[:, ::1] output=None): + """Apply minimum image convention to an array of vectors under orthorhombic boundary conditions + + Parameters + ---------- + vectors : float32 or float64 array + Vector array of shape ``(n, 3)``, either float32 or float64. These + represent many vectors (such as between two particles). + box : float32 or float64 array + periodic boundary dimensions, in 3x3 format + output: float32 or float64 array + Same shape and dtype as input. The vectors from the input, but + minimized according to the size of the box. + + Returns + ------- + minimized_vectors : np.array + Same shape and dtype as input. The vectors from the input, but + minimized according to the size of the box. + """ + cdef floating[:, ::1] output_view + cdef size_t nvals0 = vectors.shape[0] + + cdef cnp.npy_intp[1] dims + + cdef ssize_t final_size = nvals0 * (nvals0 -1) // 2 + + dims[0] = final_size + + # return early, will seg + if nvals0 == 0: + return np.array([]) + + if output is None: + # cnp.PyArray_NewLikeArray is not directly accessible + output = np.empty_like(vectors) + else: + if output.ndim != vectors.ndim: + raise ValueError("output dimensions not same as vectors") + if output.shape[0] != nvals0: + raise ValueError(f"output must be the same length as vectors ({nvals0}), you provided {output.shape[0]}") + + output_view = output + + MinimizeVectorsOrtho(&vectors[0][0], + nvals0, + &box[0], + &output_view[0][0]) + + return np.array(output) + + +@cython.boundscheck(False) +@cython.wraparound(False) +def minimize_vectors_triclinic( + floating[:, ::1] vectors, + floating[:, ::1] box, + floating[:, ::1] output=None): + """Apply minimum image convention to an array of vectors under triclinic boundary conditions + + Parameters + ---------- + vectors : float32 or float64 array + Vector array of shape ``(n, 3)``, either float32 or float64. These + represent many vectors (such as between two particles). + box : float32 or float64 array + periodic boundary dimensions, in 3x3 format + output: float32 or float64 array + Same shape and dtype as input. The vectors from the input, but + minimized according to the size of the box. + + Returns + ------- + minimized_vectors : np.array + Same shape and dtype as input. The vectors from the input, but + minimized according to the size of the box. + """ + cdef floating[:, ::1] output_view + cdef size_t nvals0 = vectors.shape[0] + + cdef cnp.npy_intp[1] dims + + cdef ssize_t final_size = nvals0 * (nvals0 -1) // 2 + + dims[0] = final_size + + # return early, will seg + if nvals0 == 0: + return np.array([]) + + if output is None: + # cnp.PyArray_NewLikeArray is not directly accessible + output = np.empty_like(vectors) + else: + if output.ndim != vectors.ndim: + raise ValueError("output dimensions not same as vectors") + if output.shape[0] != nvals0: + raise ValueError(f"output must be the same length as vectors ({nvals0}), you provided {output.shape[0]}") + + output_view = output + + MinimizeVectorsTriclinic(&vectors[0][0], + nvals0, + &box[0][0], + &output_view[0][0]) + + return np.array(output) diff --git a/distopia/tests/test_distopia.py b/distopia/tests/test_distopia.py index e9621a7..0022be6 100644 --- a/distopia/tests/test_distopia.py +++ b/distopia/tests/test_distopia.py @@ -1,3 +1,4 @@ +import itertools import pytest import numpy as np import distopia @@ -479,3 +480,75 @@ def test_periodic_angles(self, box_bonds, positions_angles, dtype): for val in [test1, test2, test3]: assert_almost_equal(ref, val, self.prec, err_msg="Min image in angle calculation failed") + + # try shifts of -2 to +2 in each dimension, and all combinations of shifts + @pytest.mark.parametrize( + "shift", itertools.product(range(-2, 3), range(-2, 3), range(-2, 3)) + ) + @pytest.mark.parametrize("dtype", (np.float32, np.float64)) + def test_minimize_vectors_ortho(self, shift, dtype): + # test vectors pointing in all directions + # these currently all obey minimum convention as they're much smaller than the box + vec = np.array( + list(itertools.product(range(-1, 2), range(-1, 2), range(-1, 2))), + dtype=dtype, + ) + box = np.asarray([[10, 0, 0], [0, 10, 0], [0, 0, 10]], dtype=dtype) + # box is 3x3 representation + # multiply by shift, then sum xyz components then add these to the vector + # this technically doesn't alter the vector because of periodic boundaries + shifted_vec = (vec + (box.T * shift).sum(axis=1)).astype(dtype) + box2 = np.asarray([10, 10, 10], dtype=dtype) + # empty input + res = distopia.minimize_vectors_ortho(np.array([[]], dtype=dtype), box2) + assert res.size == 0 + # pass invalid output buffer size + output = np.array([[]]).astype(dtype) + with pytest.raises(ValueError, match="must be the same length"): + res = distopia.minimize_vectors_ortho(shifted_vec, box2, output) + # pass output buffer + output = np.empty_like(shifted_vec) + res = distopia.minimize_vectors_ortho(shifted_vec, box2, output) + assert_equal(res, output) + assert_allclose(res, vec, atol=0.00001) + assert res.dtype == dtype + # don't pass output buffer + res = distopia.minimize_vectors_ortho(shifted_vec, box2) + assert_allclose(res, vec, atol=0.00001) + assert res.dtype == dtype + + # try shifts of -2 to +2 in each dimension, and all combinations of shifts + @pytest.mark.parametrize( + "shift", itertools.product(range(-2, 3), range(-2, 3), range(-2, 3)) + ) + @pytest.mark.parametrize("dtype", (np.float32, np.float64)) + def test_minimize_vectors_triclinic(self, shift, dtype): + # test vectors pointing in all directions + # these currently all obey minimum convention as they're much smaller than the box + vec = np.array( + list(itertools.product(range(-1, 2), range(-1, 2), range(-1, 2))), + dtype=dtype, + ) + box = np.asarray([[10, 0, 0], [2, 10, 0], [2, 2, 10]], dtype=dtype) + # box is 3x3 representation + # multiply by shift, then sum xyz components then add these to the vector + # this technically doesn't alter the vector because of periodic boundaries + shifted_vec = (vec + (box.T * shift).sum(axis=1)).astype(dtype) + box2 = np.asarray([[10, 0, 0, 2, 10, 0, 2, 2, 10]], dtype=dtype) + # empty input + res = distopia.minimize_vectors_triclinic(np.array([[]], dtype=dtype), box2) + assert res.size == 0 + # pass invalid output buffer size + output = np.array([[]]).astype(dtype) + with pytest.raises(ValueError, match="must be the same length"): + res = distopia.minimize_vectors_triclinic(shifted_vec, box2, output) + # pass output buffer + output = np.empty_like(shifted_vec) + res = distopia.minimize_vectors_triclinic(shifted_vec, box2, output) + assert_equal(res, output) + assert_allclose(res, vec, atol=0.00001) + assert res.dtype == dtype + # don't pass output buffer + res = distopia.minimize_vectors_triclinic(shifted_vec, box2) + assert_allclose(res, vec, atol=0.00001) + assert res.dtype == dtype diff --git a/docs/docs/api.md b/docs/docs/api.md index da74f40..7939fc3 100644 --- a/docs/docs/api.md +++ b/docs/docs/api.md @@ -46,4 +46,10 @@ :docstring: ::: distopia.self_distance_array_triclinic - :docstring: \ No newline at end of file + :docstring: + +::: distopia.minimize_vectors_ortho + :docstring: + +::: distopia.minimize_vectors_triclinic + :docstring: diff --git a/docs/docs/examples.md b/docs/docs/examples.md index d0a8479..b22d406 100644 --- a/docs/docs/examples.md +++ b/docs/docs/examples.md @@ -144,6 +144,41 @@ for i in range(N): ``` +## Minimizing Vectors + +### Orthorhombic periodic boundary conditions + +For orthorhombic boxes you need to provide the three box vector lengths `[l, l, l]` + +```python +import numpy as np +import distopia + +# make N x 3 coordinate arrays +N = 10000 +coordinates0 = np.random.rand(3 * N).reshape(N, 3).astype(np.float32) +coordinates1 = np.random.rand(3 * N).reshape(N, 3).astype(np.float32) +box = np.asarray([10, 10, 10]).astype(np.float32) +vectors = coordinates1 - coordinates0 +output = distopia.minimize_vectors_ortho(vectors, box) +``` + +### Triclinic periodic boundary conditions + +For triclinic boxes the box matrix must be provided in 3x3 matrix form. + +```python +import numpy as np +import distopia + +# make N x 3 coordinate arrays +N = 10000 +coordinates0 = np.random.rand(3 * N).reshape(N, 3).astype(np.float32) +coordinates1 = np.random.rand(3 * N).reshape(N, 3).astype(np.float32) +box = np.asarray([[10, 0, 0], [0, 10, 0], [0, 0, 10]]).astype(np.float32) +vectors = coordinates1 - coordinates0 +output = distopia.minimize_vectors_triclinic(vectors, box) +``` ## Questions diff --git a/libdistopia/include/distopia.h b/libdistopia/include/distopia.h index ee0dcde..80cec5d 100644 --- a/libdistopia/include/distopia.h +++ b/libdistopia/include/distopia.h @@ -48,6 +48,8 @@ namespace distopia { template void SelfDistanceArrayNoBoxIdx(const T *coords, const int *idx, int n, T *out); template void SelfDistanceArrayOrthoIdx(const T *coords, const int *idx, int n, const T *box, T *out); template void SelfDistanceArrayTriclinicIdx(const T *coords, const int *idx, int n, const T *box, T *out); + template void MinimizeVectorsOrtho(const T *vectors, int n, const T *box, T *output); + template void MinimizeVectorsTriclinic(const T *vectors, int n, const T *box, T *output); int GetNFloatLanes(); int GetNDoubleLanes(); std::vector DistopiaSupportedAndGeneratedTargets(); diff --git a/libdistopia/src/distopia.cpp b/libdistopia/src/distopia.cpp index 22dd0b5..8143932 100644 --- a/libdistopia/src/distopia.cpp +++ b/libdistopia/src/distopia.cpp @@ -112,6 +112,12 @@ namespace distopia { return hn::Sqrt(acc); } + + void MinimizeVectorsInPrimaryUnitCell(V &vx, + V &vy, + V &vz) const { + MinimiseVectors(vx, vy, vz); + } }; template , typename V = hn::VFromD> @@ -258,6 +264,13 @@ namespace distopia { return acc; } + + void MinimizeVectorsInPrimaryUnitCell(V &vx, + V &vy, + V &vz) const { + ShiftIntoPrimaryUnitCell(vx, vy, vz); + MinimiseVectors(vx, vy, vz); + } }; template @@ -1037,6 +1050,31 @@ namespace distopia { } } + template + void MinimizeVectors(const T *coords, int n, T *output, B &box) { + const hn::ScalableTag d; + int nlanes = hn::Lanes(d); + + auto a_x = hn::Undefined(d); + auto a_y = hn::Undefined(d); + auto a_z = hn::Undefined(d); + + const T *a_src; + T *dst; + + a_src = &coords[0]; + dst = output; + + for (int i=0; i void DistancesNoBoxHwy(const T *a, const T *b, int n, T *out) { hn::ScalableTag d; @@ -1248,6 +1286,18 @@ namespace distopia { SelfDistanceArrayIdx(coords, idx, n, out, vbox); } + template + void MinimizeVectorsOrthoHwy(const T *coords, int n, const T *box, T *output) { + hn::ScalableTag d; + const OrthogonalBox vbox(d, box); + MinimizeVectors(coords, n, output, vbox); + } + template + void MinimizeVectorsTriclinicHwy(const T *coords, int n, const T *box, T *output) { + hn::ScalableTag d; + const TriclinicBox vbox(d, box); + MinimizeVectors(coords, n, output, vbox); + } int GetNFloatLanes() { hn::ScalableTag d; @@ -1602,6 +1652,30 @@ namespace distopia { } HWY_EXPORT_AND_DYNAMIC_DISPATCH_T(SelfDistanceArrayTriclinicIdxHwy)(coords, idx, n, box, out); } + template <> void MinimizeVectorsOrtho(const float *vectors, int n, const float* box, float *output) { + if (n < GetNFloatLanes()) { + return distopia::N_SCALAR::MinimizeVectorsOrthoHwy(vectors, n, box, output); + } + HWY_EXPORT_AND_DYNAMIC_DISPATCH_T(MinimizeVectorsOrthoHwy)(vectors, n, box, output); + } + template <> void MinimizeVectorsOrtho(const double *vectors, int n, const double *box, double *output) { + if (n < GetNDoubleLanes()) { + return distopia::N_SCALAR::MinimizeVectorsOrthoHwy(vectors, n, box, output); + } + HWY_EXPORT_AND_DYNAMIC_DISPATCH_T(MinimizeVectorsOrthoHwy)(vectors, n, box, output); + } + template <> void MinimizeVectorsTriclinic(const float *vectors, int n, const float* box, float *output) { + if (n < GetNFloatLanes()) { + return distopia::N_SCALAR::MinimizeVectorsTriclinicHwy(vectors, n, box, output); + } + HWY_EXPORT_AND_DYNAMIC_DISPATCH_T(MinimizeVectorsTriclinicHwy)(vectors, n, box, output); + } + template <> void MinimizeVectorsTriclinic(const double *vectors, int n, const double *box, double *output) { + if (n < GetNDoubleLanes()) { + return distopia::N_SCALAR::MinimizeVectorsTriclinicHwy(vectors, n, box, output); + } + HWY_EXPORT_AND_DYNAMIC_DISPATCH_T(MinimizeVectorsTriclinicHwy)(vectors, n, box, output); + } std::vector DistopiaSupportedAndGeneratedTargets() { std::vector targets = hwy::SupportedAndGeneratedTargets(); From 0d5c620e39cd4d9d08aa5d3c3bd875daacffcfb4 Mon Sep 17 00:00:00 2001 From: Pardhav Maradani Date: Sun, 22 Feb 2026 19:05:36 +0530 Subject: [PATCH 2/5] * Add c++ layer tests and benchmarks comparing against mda --- libdistopia/src/distopia.cpp | 4 +- libdistopia/test/bench.cpp | 130 +++++++++++++++++ libdistopia/test/compare/calc_distances.hpp | 149 ++++++++++++++++++++ libdistopia/test/test_fixtures.h | 54 +++++++ libdistopia/test/test_mda_match.cpp | 34 +++++ 5 files changed, 369 insertions(+), 2 deletions(-) create mode 100644 libdistopia/test/compare/calc_distances.hpp diff --git a/libdistopia/src/distopia.cpp b/libdistopia/src/distopia.cpp index 8143932..2f8d4b5 100644 --- a/libdistopia/src/distopia.cpp +++ b/libdistopia/src/distopia.cpp @@ -1658,7 +1658,7 @@ namespace distopia { } HWY_EXPORT_AND_DYNAMIC_DISPATCH_T(MinimizeVectorsOrthoHwy)(vectors, n, box, output); } - template <> void MinimizeVectorsOrtho(const double *vectors, int n, const double *box, double *output) { + template <> void MinimizeVectorsOrtho(const double *vectors, int n, const double *box, double *output) { if (n < GetNDoubleLanes()) { return distopia::N_SCALAR::MinimizeVectorsOrthoHwy(vectors, n, box, output); } @@ -1670,7 +1670,7 @@ namespace distopia { } HWY_EXPORT_AND_DYNAMIC_DISPATCH_T(MinimizeVectorsTriclinicHwy)(vectors, n, box, output); } - template <> void MinimizeVectorsTriclinic(const double *vectors, int n, const double *box, double *output) { + template <> void MinimizeVectorsTriclinic(const double *vectors, int n, const double *box, double *output) { if (n < GetNDoubleLanes()) { return distopia::N_SCALAR::MinimizeVectorsTriclinicHwy(vectors, n, box, output); } diff --git a/libdistopia/test/bench.cpp b/libdistopia/test/bench.cpp index 9a7f135..dfc74f0 100644 --- a/libdistopia/test/bench.cpp +++ b/libdistopia/test/bench.cpp @@ -11,6 +11,7 @@ #include "distopia.h" #include "compare/calc_distances.h" +#include "compare/calc_distances.hpp" #define BOXSIZE 30 @@ -654,4 +655,133 @@ BENCHMARK_REGISTER_F(CoordinatesBench, DihedralsMDATriclinicOutBoxDouble)->Range +template class MinimizeVectorsBench : public benchmark::Fixture { +public: + // members + size_t ncoords; + size_t nelements; + T* coords0 = nullptr; + T* coords1 = nullptr; + T* vectors = nullptr; + T* output = nullptr; + T box[3]; + T triclinic_box[9]; + + void SetUp(benchmark::State &state) override { + ncoords = static_cast(state.range(0)); + InitCoords(state.range(0), BOXSIZE, state.range(1)); + } + + // coordinates range from 0 - delta to BOXSIZE + delta + void InitCoords(const size_t n_coords, const double boxsize, const double delta) { + ncoords = n_coords; + nelements = ncoords * 3; + coords0 = new T[nelements]; + coords1 = new T[nelements]; + vectors = new T[nelements]; + output = new T[nelements]; + + RandomFloatingPoint(coords0, nelements, 0 - delta, boxsize + delta); + RandomFloatingPoint(coords1, nelements, 0 - delta, boxsize + delta); + + for (size_t i = 0; i < nelements; ++i) { + vectors[i] = coords1[i] - coords0[i]; + } + + box[0] = boxsize; + box[1] = boxsize; + box[2] = boxsize; + + triclinic_box[0] = boxsize; + triclinic_box[1] = triclinic_box[2] = 0.0; + triclinic_box[3] = boxsize / 5.0; + triclinic_box[4] = boxsize; + triclinic_box[5] = 0.0; + triclinic_box[6] = triclinic_box[7] = boxsize / 5.0; + triclinic_box[8] = boxsize; + } + + void TearDown(const ::benchmark::State &state) override { + delete[] coords0; + delete[] coords1; + delete[] vectors; + delete[] output; + } + + void BM_minimize_vectors_ortho_MDA(benchmark::State &state) { + for (auto _ : state) { + _calc_minimize_vectors_ortho(vectors, ncoords, box, output); + } + state.SetItemsProcessed(ncoords * state.iterations()); + state.counters["Per Vector"] = benchmark::Counter( + ncoords * state.iterations(), + benchmark::Counter::kIsRate | benchmark::Counter::kInvert); + } + + void BM_minimize_vectors_ortho(benchmark::State &state) { + for (auto _ : state) { + distopia::MinimizeVectorsOrtho(vectors, ncoords, box, output); + } + state.SetItemsProcessed(ncoords * state.iterations()); + state.counters["Per Vector"] = benchmark::Counter( + ncoords * state.iterations(), + benchmark::Counter::kIsRate | benchmark::Counter::kInvert); + } + + void BM_minimize_vectors_triclinic_MDA(benchmark::State &state) { + for (auto _ : state) { + _calc_minimize_vectors_triclinic(vectors, ncoords, triclinic_box, output); + } + state.SetItemsProcessed(ncoords * state.iterations()); + state.counters["Per Vector"] = benchmark::Counter( + ncoords * state.iterations(), + benchmark::Counter::kIsRate | benchmark::Counter::kInvert); + } + + void BM_minimize_vectors_triclinic(benchmark::State &state) { + for (auto _ : state) { + distopia::MinimizeVectorsTriclinic(vectors, ncoords, triclinic_box, output); + } + state.SetItemsProcessed(ncoords * state.iterations()); + state.counters["Per Vector"] = benchmark::Counter( + ncoords * state.iterations(), + benchmark::Counter::kIsRate | benchmark::Counter::kInvert); + } +}; + + +BENCHMARK_TEMPLATE_DEFINE_F(MinimizeVectorsBench, MinimizeVectorsMDAOrthoFloat, float) +(benchmark::State &state) { BM_minimize_vectors_ortho_MDA(state); } +BENCHMARK_REGISTER_F(MinimizeVectorsBench, MinimizeVectorsMDAOrthoFloat)->RangeMultiplier(10)->Ranges({{10, 10000000}, {0, 0}, {0, 0}}); + +BENCHMARK_TEMPLATE_DEFINE_F(MinimizeVectorsBench, MinimizeVectorsMDAOrthoDouble, double) +(benchmark::State &state) { BM_minimize_vectors_ortho_MDA(state); } +BENCHMARK_REGISTER_F(MinimizeVectorsBench, MinimizeVectorsMDAOrthoDouble)->RangeMultiplier(10)->Ranges({{10, 10000000}, {0, 0}, {0, 0}}); + +BENCHMARK_TEMPLATE_DEFINE_F(MinimizeVectorsBench, MinimizeVectorsOrthoFloat, float) +(benchmark::State &state) { BM_minimize_vectors_ortho(state); } +BENCHMARK_REGISTER_F(MinimizeVectorsBench, MinimizeVectorsOrthoFloat)->RangeMultiplier(10)->Ranges({{10, 10000000}, {0, 0}, {0, 0}}); + +BENCHMARK_TEMPLATE_DEFINE_F(MinimizeVectorsBench, MinimizeVectorsOrthoDouble, double) +(benchmark::State &state) { BM_minimize_vectors_ortho(state); } +BENCHMARK_REGISTER_F(MinimizeVectorsBench, MinimizeVectorsOrthoDouble)->RangeMultiplier(10)->Ranges({{10, 10000000}, {0, 0}, {0, 0}}); + + +BENCHMARK_TEMPLATE_DEFINE_F(MinimizeVectorsBench, MinimizeVectorsMDATriclinicFloat, float) +(benchmark::State &state) { BM_minimize_vectors_triclinic_MDA(state); } +BENCHMARK_REGISTER_F(MinimizeVectorsBench, MinimizeVectorsMDATriclinicFloat)->RangeMultiplier(10)->Ranges({{10, 10000000}, {0, 0}, {0, 0}}); + +BENCHMARK_TEMPLATE_DEFINE_F(MinimizeVectorsBench, MinimizeVectorsMDATriclinicDouble, double) +(benchmark::State &state) { BM_minimize_vectors_triclinic_MDA(state); } +BENCHMARK_REGISTER_F(MinimizeVectorsBench, MinimizeVectorsMDATriclinicDouble)->RangeMultiplier(10)->Ranges({{10, 10000000}, {0, 0}, {0, 0}}); + +BENCHMARK_TEMPLATE_DEFINE_F(MinimizeVectorsBench, MinimizeVectorsTriclinicFloat, float) +(benchmark::State &state) { BM_minimize_vectors_triclinic(state); } +BENCHMARK_REGISTER_F(MinimizeVectorsBench, MinimizeVectorsTriclinicFloat)->RangeMultiplier(10)->Ranges({{10, 10000000}, {0, 0}, {0, 0}}); + +BENCHMARK_TEMPLATE_DEFINE_F(MinimizeVectorsBench, MinimizeVectorsTriclinicDouble, double) +(benchmark::State &state) { BM_minimize_vectors_triclinic(state); } +BENCHMARK_REGISTER_F(MinimizeVectorsBench, MinimizeVectorsTriclinicDouble)->RangeMultiplier(10)->Ranges({{10, 10000000}, {0, 0}, {0, 0}}); + + BENCHMARK_MAIN(); \ No newline at end of file diff --git a/libdistopia/test/compare/calc_distances.hpp b/libdistopia/test/compare/calc_distances.hpp new file mode 100644 index 0000000..bdc89a0 --- /dev/null +++ b/libdistopia/test/compare/calc_distances.hpp @@ -0,0 +1,149 @@ +/* + MDAnalysis --- https://www.mdanalysis.org + Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors + (see the file AUTHORS for the full list of names) + + Released under the Lesser GNU Public Licence, v2.1 or any higher version + + Please cite your use of MDAnalysis in published work: + + R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, + D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. + MDAnalysis: A Python package for the rapid analysis of molecular dynamics + simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th + Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. + doi: 10.25080/majora-629e541a-00e + + N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. + MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. + J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 +*/ + +#ifndef __DISTANCES_HPP__ +#define __DISTANCES_HPP__ + +#include + +template +void _minimize_vectors_ortho(T* dx, T* box, T* inverse_box) +{ + /* + Minimise dx to be the shortest vector + Operates in-place on dx! + + This version is near identical to minimum_image, with the + difference being that this version uses templatized types. + */ + int i; + T s; + T eps = std::numeric_limits::epsilon(); + for (i=0; i<3; i++) { + if (box[i] > eps) { + s = inverse_box[i] * dx[i]; + dx[i] = box[i] * (s - round(s)); + } + } +} + +template +static void _calc_minimize_vectors_ortho(T* vectors, uint64_t numvectors, T* box, T* output) +{ + T inverse_box[3]; + inverse_box[0] = 1.0 / box[0]; + inverse_box[1] = 1.0 / box[1]; + inverse_box[2] = 1.0 / box[2]; + +#ifdef PARALLEL +#pragma omp parallel for shared(vectors) +#endif + for (uint64_t i = 0; i < numvectors; i++) { + T dx[3]; + dx[0] = vectors[i * 3]; + dx[1] = vectors[i * 3 + 1]; + dx[2] = vectors[i * 3 + 2]; + _minimize_vectors_ortho(dx, box, inverse_box); + output[i * 3] = dx[0]; + output[i * 3 + 1] = dx[1]; + output[i * 3 + 2] = dx[2]; + } +} + +template +void _minimize_vectors_triclinic(T* dx, T* box, T* inverse_box) +{ + /* + Minimise dx to be the shortest vector + Operates in-place on dx! + + This version is near identical to minimum_image_triclinic, with the + difference being that this version does not assume that coordinates are at + most a single box length apart and uses templatized types. + */ + T dx_min[3] = {0.0, 0.0, 0.0}; + T s; + T dsq; + T dsq_min = DBL_MAX; + T rx; + T ry[2]; + T rz[3]; + int ix, iy, iz; + + // Shift into primary unit cell + s = round(inverse_box[2] * dx[2]); + dx[0] -= s * box[6]; + dx[1] -= s * box[7]; + dx[2] -= s * box[8]; + s = round(inverse_box[1] * dx[1]); + dx[0] -= s * box[3]; + dx[1] -= s * box[4]; + s = round(inverse_box[0] * dx[0]); + dx[0] -= s * box[0]; + + for (ix = -1; ix < 2; ++ix) { + rx = dx[0] + box[0] * ix; + for (iy = -1; iy < 2; ++iy) { + ry[0] = rx + box[3] * iy; + ry[1] = dx[1] + box[4] * iy; + for (iz = -1; iz < 2; ++iz) { + rz[0] = ry[0] + box[6] * iz; + rz[1] = ry[1] + box[7] * iz; + rz[2] = dx[2] + box[8] * iz; + dsq = rz[0] * rz[0] + rz[1] * rz[1] + rz[2] * rz[2]; + if (dsq < dsq_min) { + dsq_min = dsq; + dx_min[0] = rz[0]; + dx_min[1] = rz[1]; + dx_min[2] = rz[2]; + } + } + } + } + dx[0] = dx_min[0]; + dx[1] = dx_min[1]; + dx[2] = dx_min[2]; +} + +template +static void _calc_minimize_vectors_triclinic(T* vectors, uint64_t numvectors, T* box, T* output) +{ + T inverse_box[3]; + inverse_box[0] = 1.0 / box[0]; + inverse_box[1] = 1.0 / box[4]; + inverse_box[2] = 1.0 / box[8]; + +#ifdef PARALLEL +#pragma omp parallel for shared(vectors) +#endif + for (uint64_t i = 0; i < numvectors; i++) { + T dx[3]; + dx[0] = vectors[i * 3]; + dx[1] = vectors[i * 3 + 1]; + dx[2] = vectors[i * 3 + 2]; + _minimize_vectors_triclinic(dx, box, inverse_box); + output[i * 3] = dx[0]; + output[i * 3 + 1] = dx[1]; + output[i * 3 + 2] = dx[2]; + } +} + +#endif diff --git a/libdistopia/test/test_fixtures.h b/libdistopia/test/test_fixtures.h index ce87075..df39993 100644 --- a/libdistopia/test/test_fixtures.h +++ b/libdistopia/test/test_fixtures.h @@ -222,4 +222,58 @@ class CoordinatesIdx : public ::testing::Test { delete[] results; } }; + + +template +class MinimizeVectorsTest : public ::testing::Test { +public: + // members + size_t ncoords; + size_t nelements; + T* coords0 = nullptr; + T* coords1 = nullptr; + T* vectors = nullptr; + T* ref = nullptr; + T* output = nullptr; + T box[3]; + T triclinic_box[9]; + + void SetUp(const size_t n_coords, const double boxsize, const double delta) { + ncoords = n_coords; + nelements = ncoords * 3; + coords0 = new T[nelements]; + coords1 = new T[nelements]; + vectors = new T[nelements]; + ref = new T[nelements]; + output = new T[nelements]; + + RandomFloatingPoint(coords0, nelements, 0 - delta, boxsize + delta); + RandomFloatingPoint(coords1, nelements, 0 - delta, boxsize + delta); + + for (size_t i = 0; i < nelements; ++i) { + vectors[i] = coords1[i] - coords0[i]; + } + + box[0] = boxsize; + box[1] = boxsize; + box[2] = boxsize; + + triclinic_box[0] = boxsize; + triclinic_box[1] = triclinic_box[2] = 0.0; + triclinic_box[3] = boxsize / 5.0; + triclinic_box[4] = boxsize; + triclinic_box[5] = 0.0; + triclinic_box[6] = triclinic_box[7] = boxsize / 5.0; + triclinic_box[8] = boxsize; + } + + // Destructor with inlined cleanup + virtual void TearDown() { + delete[] coords0; + delete[] coords1; + delete[] vectors; + delete[] ref; + delete[] output; + } +}; #endif // DISTOPIA_TEST_FIXTURES_H \ No newline at end of file diff --git a/libdistopia/test/test_mda_match.cpp b/libdistopia/test/test_mda_match.cpp index e3f3668..ca7bb0d 100644 --- a/libdistopia/test/test_mda_match.cpp +++ b/libdistopia/test/test_mda_match.cpp @@ -7,6 +7,7 @@ #include "test_utils.h" #include "test_fixtures.h" #include "compare/calc_distances.h" +#include "compare/calc_distances.hpp" using testing::Types; typedef Types ScalarTypes; @@ -15,6 +16,7 @@ typedef Types ScalarTypes; // constants constexpr int BOXSIZE = 30; constexpr int NRESULTS = 1000; +constexpr int NCOORDS = 1000; constexpr int NINDICIES = 5; constexpr double abs_err = 1.0e-4; @@ -482,6 +484,38 @@ TYPED_TEST(DistanceArrayCoordinates, SelfDistanceArrayTriclinicScalar) { } +TYPED_TEST_SUITE(MinimizeVectorsTest, ScalarTypes); + +TYPED_TEST(MinimizeVectorsTest, MinimizeVectorsOrthoMatchesMDA) +{ + this->SetUp(NCOORDS, BOXSIZE, 3 * BOXSIZE); + + _calc_minimize_vectors_ortho(this->vectors, this->ncoords, this->box, this->ref); + + distopia::MinimizeVectorsOrtho(this->vectors, this->ncoords, this->box, this->output); + + for (std::size_t i = 0; i < this->nelements; i++) + { + EXPECT_NEAR(this->ref[i], this->output[i], abs_err); + } +} + + +TYPED_TEST(MinimizeVectorsTest, MinimizeVectorsTriclinicMatchesMDA) +{ + this->SetUp(NCOORDS, BOXSIZE, 3 * BOXSIZE); + + _calc_minimize_vectors_triclinic(this->vectors, this->ncoords, this->triclinic_box, this->ref); + + distopia::MinimizeVectorsTriclinic(this->vectors, this->ncoords, this->triclinic_box, this->output); + + for (std::size_t i = 0; i < this->nelements; i++) + { + EXPECT_NEAR(this->ref[i], this->output[i], abs_err); + } +} + + // TYPED_TEST_SUITE(CoordinatesIdx, ScalarTypes); From a5724b8456728e45392d47e9d17ed2faa3efa27b Mon Sep 17 00:00:00 2001 From: Pardhav Maradani Date: Tue, 3 Mar 2026 16:01:58 +0530 Subject: [PATCH 3/5] * Empty commit to run tests From 5b0fe5a5d37313bc953761abc5dbf06a6d10ded5 Mon Sep 17 00:00:00 2001 From: Pardhav Maradani Date: Tue, 3 Mar 2026 20:37:05 +0530 Subject: [PATCH 4/5] * Update google benchmark version to v1.9.5 --- libdistopia/googlebench | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libdistopia/googlebench b/libdistopia/googlebench index ca8d0f7..192ef10 160000 --- a/libdistopia/googlebench +++ b/libdistopia/googlebench @@ -1 +1 @@ -Subproject commit ca8d0f7b613ac915cd6b161ab01b7be449d1e1cd +Subproject commit 192ef10025eb2c4cdd392bc502f0c852196baa48 From 0678ae36049828f79e567aa0f0eedd6c253f557f Mon Sep 17 00:00:00 2001 From: Pardhav Maradani Date: Wed, 4 Mar 2026 10:01:44 +0530 Subject: [PATCH 5/5] * Temporarily add -v -p no:logging to pytest to debug macos,3.14 --- .github/workflows/make_and_test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/make_and_test.yml b/.github/workflows/make_and_test.yml index bb6e64b..e97e8cf 100644 --- a/.github/workflows/make_and_test.yml +++ b/.github/workflows/make_and_test.yml @@ -96,7 +96,7 @@ jobs: - name: python_test # run python API tests - run: pytest distopia/tests + run: pytest -v -p no:logging distopia/tests external_hwy_no_test: