diff --git a/src/ivoc/ocmatrix.cpp b/src/ivoc/ocmatrix.cpp index a735199fa2..d70efd74fa 100644 --- a/src/ivoc/ocmatrix.cpp +++ b/src/ivoc/ocmatrix.cpp @@ -41,7 +41,7 @@ static void Vect2VEC(Vect* v1, VEC& v2) { } OcMatrix::OcMatrix(int type) { - obj_ = NULL; + obj_ = nullptr; type_ = type; } OcMatrix::~OcMatrix() {} @@ -97,8 +97,8 @@ OcFullMatrix* OcMatrix::full() { OcFullMatrix::OcFullMatrix(int nrow, int ncol) : OcMatrix(MFULL) { - lu_factor_ = NULL; - lu_pivot_ = NULL; + lu_factor_ = nullptr; + lu_pivot_ = nullptr; m_ = m_get(nrow, ncol); } OcFullMatrix::~OcFullMatrix() { @@ -165,7 +165,7 @@ void OcFullMatrix::symmeigen(Matrix* mout, Vect* vout) { void OcFullMatrix::svd1(Matrix* u, Matrix* v, Vect* d) { VEC v1; Vect2VEC(d, v1); - svd(m_, u ? u->full()->m_ : NULL, v ? v->full()->m_ : NULL, &v1); + svd(m_, u ? u->full()->m_ : nullptr, v ? v->full()->m_ : nullptr, &v1); } void OcFullMatrix::getrow(int k, Vect* out) { @@ -193,6 +193,8 @@ void OcFullMatrix::getdiag(int k, Vect* out) { #endif } } else { + // Yes for negative diagonal we set the vector from the middle + // The output vector should ALWAYS be the size of biggest diagonal for (i = -k, j = 0; i < row && j < col; ++i, ++j) { #ifdef WIN32 v_elem(out, i) = m_entry(m_, i, j); @@ -228,6 +230,8 @@ void OcFullMatrix::setdiag(int k, Vect* in) { #endif } } else { + // Yes for negative diagonal we set the vector from the middle + // The input vector should ALWAYS have `nrows` elements. for (i = -k, j = 0; i < row && j < col; ++i, ++j) { #ifdef WIN32 m_set_val(m_, i, j, v_elem(in, i)); @@ -313,15 +317,6 @@ double OcFullMatrix::det(int* e) { PERM* piv = px_get(n); m_copy(m_, lu); LUfactor(lu, piv); -#if 0 -printf("LU\n"); -for (int i = 0; i < n; ++i) { - for (int j = 0; j < n; ++j) { - printf(" %g", lu->me[i][j]); - } - printf("\t%d\n", piv->pe[i]); -} -#endif double m = 1.0; *e = 0; for (int i = 0; i < n; ++i) { @@ -364,8 +359,8 @@ OcSparseMatrix::OcSparseMatrix(int nrow, int ncol) int len = 4; m_ = sp_get(nrow, ncol, len); - lu_factor_ = NULL; - lu_pivot_ = NULL; + lu_factor_ = nullptr; + lu_pivot_ = nullptr; } OcSparseMatrix::~OcSparseMatrix() { if (lu_factor_) { @@ -375,14 +370,14 @@ OcSparseMatrix::~OcSparseMatrix() { SP_FREE(m_); } -// returns pointer to sparse element. NULL if it does not exist. +// returns pointer to sparse element. nullptr if it does not exist. double* OcSparseMatrix::pelm(int i, int j) { SPROW* r = m_->row + i; int idx = sprow_idx(r, j); if (idx >= 0) { return &r->elt[idx].val; } else { - return NULL; + return nullptr; } } @@ -445,7 +440,7 @@ void OcSparseMatrix::setrow(int k, Vect* in) { int i, n = ncol(); double* p; for (i = 0; i < n; ++i) { - if ((p = pelm(k, i)) != NULL) { + if ((p = pelm(k, i)) != nullptr) { #ifdef WIN32 *p = v_elem(in, i); } else if (v_elem(in, i)) { @@ -465,7 +460,7 @@ void OcSparseMatrix::setcol(int k, Vect* in) { int i, n = nrow(); double* p; for (i = 0; i < n; ++i) { - if ((p = pelm(i, k)) != NULL) { + if ((p = pelm(i, k)) != nullptr) { #ifdef WIN32 *p = v_elem(in, i); } else if (v_elem(in, i)) { @@ -486,7 +481,7 @@ void OcSparseMatrix::setdiag(int k, Vect* in) { double* p; if (k >= 0) { for (i = 0, j = k; i < row && j < col; ++i, ++j) { - if ((p = pelm(i, j)) != NULL) { + if ((p = pelm(i, j)) != nullptr) { #ifdef WIN32 *p = v_elem(in, i); } else if (v_elem(in, i)) { @@ -500,7 +495,9 @@ void OcSparseMatrix::setdiag(int k, Vect* in) { } } else { for (i = -k, j = 0; i < row && j < col; ++i, ++j) { - if ((p = pelm(i, j)) != NULL) { + // Yes for negative diagonal we set the vector from the middle + // The input vector should ALWAYS have `nrows` elements. + if ((p = pelm(i, j)) != nullptr) { #ifdef WIN32 *p = v_elem(in, i); } else if (v_elem(in, i)) { diff --git a/src/ivoc/ocmatrix.h b/src/ivoc/ocmatrix.h index 994c691176..b1545aa04e 100644 --- a/src/ivoc/ocmatrix.h +++ b/src/ivoc/ocmatrix.h @@ -24,7 +24,7 @@ class OcMatrix { virtual double* mep(int i, int j) { unimp(); - return NULL; + return nullptr; } // matrix element pointer inline double& operator()(int i, int j) { return *mep(i, j); @@ -202,7 +202,7 @@ class OcSparseMatrix: public OcMatrix { // type 2 virtual ~OcSparseMatrix(); virtual double* mep(int, int); - virtual double* pelm(int, int); // NULL if element does not exist + virtual double* pelm(int, int); // nullptr if element does not exist virtual int nrow(); virtual int ncol(); virtual double getval(int, int); diff --git a/src/nrniv/CMakeLists.txt b/src/nrniv/CMakeLists.txt index 0b4bd3c283..c296eaedae 100644 --- a/src/nrniv/CMakeLists.txt +++ b/src/nrniv/CMakeLists.txt @@ -502,7 +502,7 @@ if(NRN_ENABLE_INTERVIEWS) include_directories(${IV_INCLUDE_DIR}) target_link_libraries(nrniv_lib interviews) else() - include_directories(nrniv_lib ${NRN_IVOS_SRC_DIR} ${PROJECT_BINARY_DIR}/src/ivos) + target_include_directories(nrniv_lib PUBLIC ${NRN_IVOS_SRC_DIR} ${PROJECT_BINARY_DIR}/src/ivos) endif() if(IV_ENABLE_X11_DYNAMIC) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 2892ab648e..0fb11833a2 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -17,6 +17,7 @@ add_executable( testneuron common/catch2_main.cpp unit_tests/basic.cpp + unit_tests/matrix.cpp unit_tests/container/container.cpp unit_tests/container/generic_data_handle.cpp unit_tests/container/mechanism.cpp diff --git a/test/unit_tests/matrix.cpp b/test/unit_tests/matrix.cpp new file mode 100644 index 0000000000..fc91b654ca --- /dev/null +++ b/test/unit_tests/matrix.cpp @@ -0,0 +1,404 @@ +#include +#include +#include +#include "ivocvect.h" + +#include +using namespace Catch::literals; + +template +class ApproxOrOpposite: public Catch::MatcherBase> { + std::vector vec; + + public: + ApproxOrOpposite(std::vector vec) + : vec(vec) {} + + bool match(std::vector const& in) const override { + if (in.size() != vec.size()) { + return false; + } + bool matched = true; + for (int i = 0; i < in.size(); ++i) { + if (in[i] != approx(vec[i])) { + matched = false; + break; + } + } + if (matched) { + return true; + } + matched = true; + for (int i = 0; i < in.size(); ++i) { + if (in[i] != approx(-1 * vec[i])) { + matched = false; + break; + } + } + + return matched; + } + + std::string describe() const override { + std::ostringstream ss; + ss << "is not approx or opposite approx of " << Catch::Detail::stringify(vec); + return ss.str(); + } + template ::value>::type> + ApproxOrOpposite& epsilon(T const& newEpsilon) { + approx.epsilon(newEpsilon); + return *this; + } + template ::value>::type> + ApproxOrOpposite& margin(T const& newMargin) { + approx.margin(newMargin); + return *this; + } + template ::value>::type> + ApproxOrOpposite& scale(T const& newScale) { + approx.scale(newScale); + return *this; + } + + mutable Catch::Detail::Approx approx = Catch::Detail::Approx::custom(); +}; + +bool compareMatrix(OcMatrix& m, const std::vector>& ref) { + REQUIRE(m.nrow() == ref.size()); + for (int i = 0; i < m.nrow(); ++i) { + REQUIRE(m.ncol() == ref[i].size()); + for (int j = 0; j < m.ncol(); ++j) { + REQUIRE(m.getval(i, j) == Catch::Detail::Approx(ref[i][j]).margin(1e-10)); + } + } + return true; +} + +SCENARIO("A Matrix", "[neuron_ivoc][OcMatrix]") { + GIVEN("A 3x3 Full matrix") { + OcFullMatrix m{3, 3}; + REQUIRE(m.nrow() == 3); + REQUIRE(m.ncol() == 3); + { + m.ident(); + REQUIRE(compareMatrix(m, {{1., 0., 0.}, {0., 1., 0}, {0., 0., 1.}})); + } + { + double* value = m.mep(0, 0); + REQUIRE(*value == 1); + *value = 3; + REQUIRE(m.getval(0, 0) == 3); + } + m.resize(4, 3); + { + m.setrow(3, 2.0); + REQUIRE(compareMatrix(m, {{3., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}, {2., 2., 2.}})); + } + { + std::vector x, y; + m.nonzeros(x, y); + std::vector res_x = {0, 1, 2, 3, 3, 3}; + std::vector res_y = {0, 1, 2, 0, 1, 2}; + REQUIRE(x == res_x); + REQUIRE(y == res_y); + } + { + m.setcol(1, 4.0); + REQUIRE(compareMatrix(m, {{3., 4., 0.}, {0., 4., 0.}, {0., 4., 1.}, {2., 4., 2.}})); + } + { + m.setdiag(0, 5.0); + REQUIRE(compareMatrix(m, {{5., 4., 0.}, {0., 5., 0.}, {0., 4., 5.}, {2., 4., 2.}})); + } + { + m.setdiag(1, 6.0); + REQUIRE(compareMatrix(m, {{5., 6., 0.}, {0., 5., 6.}, {0., 4., 5.}, {2., 4., 2.}})); + } + { + m.setdiag(-1, 7.0); + REQUIRE(compareMatrix(m, {{5., 6., 0.}, {7., 5., 6.}, {0., 7., 5.}, {2., 4., 7.}})); + } + + { + OcFullMatrix n(4, 3); + n.ident(); + m.add(&n, &m); + REQUIRE(compareMatrix(m, {{6., 6., 0.}, {7., 6., 6.}, {0., 7., 6.}, {2., 4., 7.}})); + } + { + OcFullMatrix n(4, 3); + m.bcopy(&n, 1, 1, 3, 2, 0, 0); + REQUIRE(compareMatrix(n, {{6., 6., 0.}, {7., 6., 0.}, {4., 7., 0.}, {0., 0., 0.}})); + } + { + OcFullMatrix n(4, 3); + m.transpose(&n); + REQUIRE(compareMatrix(n, {{6., 7., 0., 2.}, {6., 6., 7., 4.}, {0., 6., 6., 7.}})); + } + { + IvocVect v(3); + m.getrow(2, &v); + REQUIRE_THAT(v.vec(), Catch::Matchers::Approx(std::vector({0., 7., 6.}))); + m.setrow(0, &v); + REQUIRE(compareMatrix(m, {{0., 7., 6.}, {7., 6., 6.}, {0., 7., 6.}, {2., 4., 7.}})); + } + { + IvocVect v(4); + m.getcol(2, &v); + REQUIRE_THAT(v.vec(), Catch::Matchers::Approx(std::vector({6., 6., 6., 7.}))); + m.setcol(1, &v); + REQUIRE(compareMatrix(m, {{0., 6., 6.}, {7., 6., 6.}, {0., 6., 6.}, {2., 7., 7.}})); + } + { + m.resize(3, 3); + REQUIRE(compareMatrix(m, {{0., 6., 6.}, {7., 6., 6.}, {0., 6., 6.}})); + } + { + OcFullMatrix n(4, 3); + m.exp(&n); + REQUIRE(n(0, 0) == Catch::Detail::Approx(442925.)); + REQUIRE(compareMatrix(n, + {{442925., 938481., 938481.}, + {651970., 1381407., 1381407.}, + {442926., 938481., 938482.}})); + } + { + m.pow(2, &m); + REQUIRE(compareMatrix(m, {{42., 72., 72.}, {42., 114., 114.}, {42., 72., 72.}})); + } + { // Mescach computing of det has an error + // int e{}; + // double det = m.det(&e); + // REQUIRE(det == 0.); + // REQUIRE(e == 0); + } + *m.mep(2, 0) = 1; + *m.mep(2, 2) = 2; + { + // Mescach computing of det has an error + // int e{}; + // double det = m.det(&e); + // REQUIRE(det == -1.2348_a); + // REQUIRE(e == 5); + } { + OcFullMatrix n(4, 3); + m.inverse(&n); + n.resize(3, 3); // ??? + REQUIRE(compareMatrix(n, + {{0.064625850, -0.040816326, 0.}, + {-0.00024295432, -0.0000971817, 0.01428571}, + {-0.023566569, 0.0239067055, -0.014285714}})); + n.zero(); + REQUIRE(compareMatrix(n, {{0., 0., 0.}, {0., 0., 0.}, {0., 0., 0.}})); + } + { + IvocVect v(3); + m.getdiag(1, &v); + REQUIRE_THAT(v.vec(), Catch::Matchers::Approx(std::vector({72., 114., 0.}))); + v.vec() = {0., 72., 114.}; + m.setdiag(-1, &v); + REQUIRE(compareMatrix(m, {{42., 72., 72.}, {72., 114., 114.}, {1., 114., 2.}})); + } + { + IvocVect v(3); + m.getdiag(-2, &v); + REQUIRE(v.vec()[2] == Catch::Detail::Approx({1.})); + v.vec() = {1., 0., 0.}; + m.setdiag(2, &v); + REQUIRE(compareMatrix(m, {{42., 72., 1.}, {72., 114., 114.}, {1., 114., 2.}})); + } + { + IvocVect v(3); + v.vec() = {1, 1, 1}; + IvocVect vout(3); + m.mulv(&v, &vout); + REQUIRE_THAT(vout.vec(), + Catch::Matchers::Approx(std::vector({115., 300., 117.}))); + } + { + OcFullMatrix n(3, 3); + m.copy(&n); + REQUIRE(compareMatrix(n, {{42., 72., 1.}, {72., 114., 114.}, {1., 114., 2.}})); + OcFullMatrix o(3, 3); + m.mulm(&n, &o); + REQUIRE(compareMatrix( + o, {{6949., 11346., 8252.}, {11346., 31176., 13296.}, {8252., 13296., 13001.}})); + } + { + OcFullMatrix n(3, 3); + m.muls(2, &n); + REQUIRE(compareMatrix(n, {{84., 144., 2.}, {144., 228., 228.}, {2., 228., 4.}})); + } + { + IvocVect v(3); + v.vec() = {1, 1, 1}; + IvocVect vout(3); + m.solv(&v, &vout, false); + REQUIRE_THAT(vout.vec(), + Catch::Matchers::Approx( + std::vector({0.0088700, 0.0087927, -0.00562299}))); + m.solv(&v, &vout, true); + REQUIRE_THAT(vout.vec(), + Catch::Matchers::Approx( + std::vector({0.0088700, 0.0087927, -0.00562299}))); + } + { + IvocVect v(3); + OcFullMatrix n(3, 3); + v.vec() = {1, 2, 3}; + m.setrow(0, &v); + v.vec() = {2, 1, 4}; + m.setrow(1, &v); + v.vec() = {3, 4, 1}; + m.setrow(2, &v); + m.symmeigen(&n, &v); + REQUIRE_THAT(v.vec(), + Catch::Matchers::Approx( + std::vector({7.074673, -0.88679, -3.18788}))); + n.getcol(0, &v); + REQUIRE_THAT(v.vec(), ApproxOrOpposite({0.50578, 0.5843738, 0.634577})); + n.getcol(1, &v); + REQUIRE_THAT(v.vec(), ApproxOrOpposite({-0.8240377, 0.544925, 0.154978})); + n.getcol(2, &v); + REQUIRE_THAT(v.vec(), ApproxOrOpposite({-0.255231, -0.601301, 0.7571611})); + } + { + m.resize(2, 2); + OcFullMatrix u(2, 2); + OcFullMatrix v(2, 2); + IvocVect d(2); + m.svd1(&u, &v, &d); + REQUIRE_THAT(d.vec(), Catch::Matchers::Approx(std::vector({3., 1.}))); + // For comparison of u and v and problems with signs, see: + // https://www.educative.io/blog/sign-ambiguity-in-singular-value-decomposition + IvocVect c(4); + c.vec() = {u(0, 0), u(0, 1), v(0, 0), v(0, 1)}; + CHECK_THAT(c.vec(), ApproxOrOpposite({0.70710, 0.70710, 0.70710, 0.70710})); + c.vec() = {u(1, 0), u(1, 1), v(1, 0), v(1, 1)}; + CHECK_THAT(c.vec(), ApproxOrOpposite({0.70710, -0.70710, -0.70710, 0.70710})); + } + { + m.resize(2, 3); + { + IvocVect s(3); + s.vec() = {3., 2., 2.}; + m.setrow(0, &s); + s.vec() = {2., 3., -2.}; + m.setrow(1, &s); + } + OcFullMatrix u(2, 2); + OcFullMatrix v(3, 3); + IvocVect d(2); + m.svd1(&u, &v, &d); + REQUIRE_THAT(d.vec(), Catch::Matchers::Approx(std::vector({5., 3.}))); + // For comparison of u and v and problems with signs, see: + // https://www.educative.io/blog/sign-ambiguity-in-singular-value-decomposition + IvocVect c(5); + c.vec() = {u(0, 0), u(0, 1), v(0, 0), v(0, 1), v(0, 2)}; + CHECK_THAT(c.vec(), + ApproxOrOpposite({0.70710, 0.70710, 0.70710, 0.70710, 0.}).margin(1e-10)); + c.vec() = {u(1, 0), u(1, 1), v(1, 0), v(1, 1), v(1, 2)}; + CHECK_THAT(c.vec(), + ApproxOrOpposite({0.70710, -0.70710, 0.235702, -0.235702, 0.942809})); + c.vec() = {0., 0., v(2, 0), v(2, 1), v(2, 2)}; + CHECK_THAT(c.vec(), ApproxOrOpposite({0., 0., 0.66666, -0.66666, -0.3333333})); + } + } + GIVEN("A 3x3 Sparse matrix") { + OcSparseMatrix m{3, 3}; + REQUIRE(m.nrow() == 3); + REQUIRE(m.ncol() == 3); + { + m.ident(); + REQUIRE(compareMatrix(m, {{1., 0., 0.}, {0., 1., 0}, {0., 0., 1.}})); + REQUIRE(m.sprowlen(1) == 1); + } + { + std::vector x, y, result = {0, 1, 2}; + m.nonzeros(x, y); + REQUIRE(x == result); + REQUIRE(y == result); + } + { + double* pmep = m.mep(1, 1); + REQUIRE(*pmep == 1); + pmep = m.mep(1, 0); + REQUIRE(*pmep == 0); + pmep = m.pelm(0, 1); + REQUIRE(pmep == nullptr); + } + { + int col{}; + double value = m.spgetrowval(2, 0, &col); + REQUIRE(col == 2); + REQUIRE(value == Catch::Detail::Approx(1.0)); + } + { // m.zero() don't erase the matrix but only replace existing values by 0. + m.zero(); + REQUIRE(m.sprowlen(2) == 1); + REQUIRE(compareMatrix(m, {{0., 0., 0.}, {0., 0., 0}, {0., 0., 0.}})); + } + { + m.setrow(1, 2); + REQUIRE(compareMatrix(m, {{0., 0., 0.}, {2., 2., 2.}, {0., 0., 0.}})); + } + { + m.setcol(0, 3); + REQUIRE(compareMatrix(m, {{3., 0., 0.}, {3., 2., 2.}, {3., 0., 0.}})); + } + { + m.setdiag(0, 1); + REQUIRE(compareMatrix(m, {{1., 0., 0.}, {3., 1., 2.}, {3., 0., 1.}})); + } + { + m.setdiag(-1, 4); + REQUIRE(compareMatrix(m, {{1., 0., 0.}, {4., 1., 2.}, {3., 4., 1.}})); + } + { + m.setdiag(2, 5); + REQUIRE(compareMatrix(m, {{1., 0., 5.}, {4., 1., 2.}, {3., 4., 1.}})); + } + REQUIRE(m.sprowlen(1) == 3); + + { + IvocVect v(3); + v.vec() = {1, 2, 3}; + m.setrow(0, &v); + REQUIRE(compareMatrix(m, {{1., 2., 3.}, {4., 1., 2.}, {3., 4., 1.}})); + } + { + IvocVect v(3); + v.vec() = {1, 2, 3}; + m.setcol(0, &v); + REQUIRE(compareMatrix(m, {{1., 2., 3.}, {2., 1., 2.}, {3., 4., 1.}})); + } + { + IvocVect v(3); + v.vec() = {1, 2, 3}; + m.setdiag(0, &v); + REQUIRE(compareMatrix(m, {{1., 2., 3.}, {2., 2., 2.}, {3., 4., 3.}})); + } + { + IvocVect v(3); + v.vec() = {0., 1., 2.}; + m.setdiag(-1, &v); + REQUIRE(compareMatrix(m, {{1., 2., 3.}, {1., 2., 2.}, {3., 2., 3.}})); + } + { + IvocVect v(3); + v.vec() = {1, 2, 3}; + IvocVect out(3); + m.mulv(&v, &out); + REQUIRE_THAT(out.vec(), Catch::Matchers::Approx(std::vector({14., 11., 16.}))); + } + { + IvocVect v(3); + v.vec() = {1, 1, 1}; + IvocVect vout(3); + m.solv(&v, &vout, false); + REQUIRE_THAT(vout.vec(), Catch::Matchers::Approx(std::vector({0., 0.5, 0.}))); + m.solv(&v, &vout, true); + REQUIRE_THAT(vout.vec(), Catch::Matchers::Approx(std::vector({0., 0.5, 0.}))); + } + } +}