diff --git a/changelog b/changelog index 736adc6..87d3923 100644 --- a/changelog +++ b/changelog @@ -1,5 +1,6 @@ Version 9.0.5 - Added exists to compound_factory +- Added sub_matrix, fix and extend determinant calculation Version 9.0.4 - Fix various stopping and reconstruction errors diff --git a/include/cif++/matrix.hpp b/include/cif++/matrix.hpp index 909aec4..192f24c 100644 --- a/include/cif++/matrix.hpp +++ b/include/cif++/matrix.hpp @@ -124,6 +124,23 @@ class matrix_expression return os; } + + template + constexpr bool operator==(const matrix_expression &m) const + { + bool same = false; + if (dim_m() == m.dim_m() and dim_n() == m.dim_n()) + { + same = true; + for (std::size_t i = 0; same and i < m.dim_m(); ++i) + { + for (std::size_t j = 0; same and j < m.dim_n(); ++j) + same = operator()(i, j) == m(i, j); + } + } + + return same; + } }; // -------------------------------------------------------------------- @@ -594,6 +611,35 @@ auto operator*(const matrix_expression &m1, const matrix_expression &m2) // -------------------------------------------------------------------- +template +class sub_matrix : public matrix_expression> +{ + public: + sub_matrix(const M2 &m, int i, int j) + : m_m(m) + , m_i(i) + , m_j(j) + { + } + + constexpr std::size_t dim_m() const { return m_m.dim_m() - 1; } ///< Return dimension m + constexpr std::size_t dim_n() const { return m_m.dim_n() - 1; } ///< Return dimension n + + /** Access to the value of element [ @a i, @a j ] */ + constexpr auto operator()(std::size_t i, std::size_t j) const + { + return m_m( + i >= m_i ? i + 1 : i, + j >= m_j ? j + 1 : j); + } + + private: + const M2 &m_m; + std::size_t m_i, m_j; +}; + +// -------------------------------------------------------------------- + /** Generic routine to calculate the determinant of a matrix * * @note This is currently only implemented for fixed matrices of size 3x3 @@ -605,11 +651,23 @@ auto determinant(const M &m); template auto determinant(const matrix3x3 &m) { - return (m(0, 0) * (m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1)) + - m(0, 1) * (m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2)) + - m(0, 2) * (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0))); + return (m(0, 0) * ((m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1))) + + m(0, 1) * ((m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2))) + + m(0, 2) * ((m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0)))); } +/** Implementation of the determinant function for fixed size matrices of size 4x4 */ +template +F determinant(const matrix4x4 &m) +{ + return m(0, 0) * determinant(matrix3x3(sub_matrix(m, 0, 0))) - + m(0, 1) * determinant(matrix3x3(sub_matrix(m, 0, 1))) + + m(0, 2) * determinant(matrix3x3(sub_matrix(m, 0, 2))) - + m(0, 3) * determinant(matrix3x3(sub_matrix(m, 0, 3))); +} + +// -------------------------------------------------------------------- + /** Generic routine to calculate the inverse of a matrix * * @note This is currently only implemented for fixed matrices of size 3x3 diff --git a/include/cif++/point.hpp b/include/cif++/point.hpp index b4e8a6e..7e5d612 100644 --- a/include/cif++/point.hpp +++ b/include/cif++/point.hpp @@ -30,6 +30,7 @@ #include #include #include +#include #include #include @@ -365,11 +366,18 @@ class quaternion_type } /// \brief test for all zero values - constexpr operator bool() const + constexpr explicit operator bool() const { return a != 0 or b != 0 or c != 0 or d != 0; } + /// \brief for debugging e.g. + friend std::ostream &operator<<(std::ostream &os, const quaternion_type &rhs) + { + os << std::format("{{ a: {}, b: {}, c: {}, d: {} }}", rhs.a, rhs.b, rhs.c, rhs.d); + return os; + } + private: value_type a, b, c, d; }; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index e75820a..5fadf31 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -49,6 +49,7 @@ list( spinner reconstruction validate-pdbx + matrix ) add_library(test-main OBJECT "${CMAKE_CURRENT_SOURCE_DIR}/test-main.cpp") diff --git a/test/matrix-test.cpp b/test/matrix-test.cpp new file mode 100644 index 0000000..4320aac --- /dev/null +++ b/test/matrix-test.cpp @@ -0,0 +1,101 @@ +/*- + * SPDX-License-Identifier: BSD-2-Clause + * + * Copyright (c) 2025 NKI/AVL, Netherlands Cancer Institute + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#include "cif++/matrix.hpp" +#include "test-main.hpp" + +#include +#include + +TEST_CASE("m1") +{ + cif::matrix3x3 m = cif::identity_matrix(3); + + CHECK(cif::determinant(m) == 1); +} + +TEST_CASE("m2") +{ + cif::matrix4x4 m = cif::identity_matrix(4); + + cif::sub_matrix> ms(m, 1, 1); + CHECK(ms == cif::identity_matrix(3)); +} + +TEST_CASE("m3") +{ + cif::matrix4x4 m{ + { 1, 2, 3, 4, // + 5, 6, 7, 8, // + 9, 10, 11, 12, // + 13, 14, 15, 16 } + }; + cif::sub_matrix> ms(m, 1, 1); + + cif::matrix3x3 t{ + { 1, 3, 4, 9, 11, 12, 13, 15, 16 } + }; + + CHECK(ms == t); +} + +TEST_CASE("m4") +{ + cif::matrix4x4 m{ + { + -2, + 3, + 1, + 0, + 4, + 1, + -3, + 2, + 0, + -1, + 2, + 5, + 3, + 2, + 0, + -4, + } + }; + + std::cout << m << "\n\n"; + + std::cout << cif::matrix3x3(cif::sub_matrix(m, 0, 0)) << "\n\n"; + std::cout << cif::matrix3x3(cif::sub_matrix(m, 0, 1)) << "\n\n"; + std::cout << cif::matrix3x3(cif::sub_matrix(m, 0, 2)) << "\n\n"; + std::cout << cif::matrix3x3(cif::sub_matrix(m, 0, 3)) << "\n\n"; + + std::cout << cif::determinant(cif::matrix3x3(cif::sub_matrix(m, 0, 0))) << "\n\n"; + std::cout << cif::determinant(cif::matrix3x3(cif::sub_matrix(m, 0, 1))) << "\n\n"; + std::cout << cif::determinant(cif::matrix3x3(cif::sub_matrix(m, 0, 2))) << "\n\n"; + std::cout << cif::determinant(cif::matrix3x3(cif::sub_matrix(m, 0, 3))) << "\n\n"; + + CHECK(cif::determinant(m) == 332); +}