F5/unison
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

core: Additional MatrixArray utilities

Browse Source 
New features include determinant, Frobenius norm, identity matrix, matrix page copies and matrix page joining
This commit is contained in:
Gabriel Ferreira
2024-02-29 15:52:21 +01:00
parent 49a6e2ef7b
commit 95f959932b
3 changed files with 355 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

142
src/core/model/matrix-array.cc
View File

@@ -156,6 +156,70 @@ MatrixArray<T>::Transpose() const
return res;
}
template <class T>
MatrixArray<T>
MatrixArray<T>::Determinant() const
{
MatrixArray<T> res{1, 1, m_numPages};
NS_ASSERT_MSG(m_numRows == m_numCols, "Matrix is not square");
// In case of small matrices, we use a fast path
if (m_numRows == 1)
{
return *this;
}
// Calculate determinant for each matrix
for (size_t page = 0; page < m_numPages; ++page)
{
res(0, 0, page) = 0;
auto pageValues = GetPagePtr(page);
// Fast path for 2x2 matrices
if (m_numRows == 2)
{
res(0, 0, page) = pageValues[0] * pageValues[3] - pageValues[1] * pageValues[2];
continue;
}
for (size_t detN = 0; detN < m_numRows; detN++)
{
auto partDetP = T{0} + 1.0;
auto partDetN = T{0} + 1.0;
for (size_t row = 0; row < m_numRows; row++)
{
// Wraparound not to have to extend the matrix
// Positive determinant
size_t col = (row + detN) % m_numCols;
partDetP *= pageValues[row * m_numCols + col];
// Negative determinant
col = m_numCols - 1 - (row + detN) % m_numCols;
partDetN *= pageValues[row * m_numCols + col];
}
res(0, 0, page) += partDetP - partDetN;
}
}
return res;
}
template <class T>
MatrixArray<T>
MatrixArray<T>::FrobeniusNorm() const
{
MatrixArray<T> res{1, 1, m_numPages};
for (size_t page = 0; page < m_numPages; ++page)
{
// Calculate the sum of squared absolute values of each matrix page
res[page] = 0;
auto pagePtr = this->GetPagePtr(page);
for (size_t i = 0; i < m_numRows * m_numCols; i++)
{
auto absVal = std::abs(pagePtr[i]);
res[page] += absVal * absVal;
}
res[page] = sqrt(res[page]);
}
return res;
}
template <class T>
MatrixArray<T>
MatrixArray<T>::MultiplyByLeftAndRightMatrix(const MatrixArray<T>& lMatrix,
@@ -229,6 +293,84 @@ MatrixArray<T>::HermitianTranspose() const
return retMatrix;
}
template <class T>
MatrixArray<T>
MatrixArray<T>::MakeNCopies(size_t nCopies) const
{
NS_ASSERT_MSG(m_numPages == 1, "The MatrixArray should have only one page to be copied.");
auto copiedMatrix = MatrixArray<T>{m_numRows, m_numCols, nCopies};
for (size_t copy = 0; copy < nCopies; copy++)
{
for (size_t i = 0; i < m_numRows * m_numCols; i++)
{
copiedMatrix.GetPagePtr(copy)[i] = m_values[i];
}
}
return copiedMatrix;
}
template <class T>
MatrixArray<T>
MatrixArray<T>::ExtractPage(size_t page) const
{
NS_ASSERT_MSG(page < m_numPages, "The page to extract from the MatrixArray is out of bounds.");
auto extractedPage = MatrixArray<T>{m_numRows, m_numCols, 1};
for (size_t i = 0; i < m_numRows * m_numCols; ++i)
{
extractedPage.m_values[i] = GetPagePtr(page)[i];
}
return extractedPage;
}
template <class T>
MatrixArray<T>
MatrixArray<T>::JoinPages(const std::vector<MatrixArray<T>>& pages)
{
auto jointMatrix =
MatrixArray<T>{pages.front().GetNumRows(), pages.front().GetNumCols(), pages.size()};
for (size_t page = 0; page < jointMatrix.GetNumPages(); page++)
{
NS_ASSERT_MSG(pages[page].GetNumRows() == jointMatrix.GetNumRows(),
"All page matrices should have the same number of rows");
NS_ASSERT_MSG(pages[page].GetNumCols() == jointMatrix.GetNumCols(),
"All page matrices should have the same number of columns");
NS_ASSERT_MSG(pages[page].GetNumPages() == 1,
"All page matrices should have a single page");
size_t i = 0;
for (auto a : pages[page].GetValues())
{
jointMatrix.GetPagePtr(page)[i] = a;
i++;
}
}
return jointMatrix;
}
template <class T>
MatrixArray<T>
MatrixArray<T>::IdentityMatrix(const size_t size, const size_t pages)
{
auto identityMatrix = MatrixArray<T>{size, size, pages};
for (std::size_t page = 0; page < pages; page++)
{
for (std::size_t i = 0; i < size; i++)
{
identityMatrix(i, i, page) = 1.0;
}
}
return identityMatrix;
}
template <class T>
MatrixArray<T>
MatrixArray<T>::IdentityMatrix(const MatrixArray& likeme)
{
NS_ASSERT_MSG(likeme.GetNumRows() == likeme.GetNumCols(), "Template array is not square.");
return IdentityMatrix(likeme.GetNumRows(), likeme.GetNumPages());
}
template MatrixArray<std::complex<double>> MatrixArray<std::complex<double>>::HermitianTranspose()
const;
template class MatrixArray<std::complex<double>>;

44
src/core/model/matrix-array.h
View File

@@ -189,6 +189,16 @@ class MatrixArray : public ValArray<T>
* \return The resulting MatrixArray composed of the array of transposed matrices.
*/
MatrixArray Transpose() const;
/**
* \brief This operator calculates a vector o determinants, one for each page
* \return The resulting MatrixArray with the determinants for each page.
*/
MatrixArray Determinant() const;
/**
* \brief This operator calculates a vector of Frobenius norm, one for each page
* \return The resulting MatrixArray with the Frobenius norm for each page.
*/
MatrixArray FrobeniusNorm() const;
/**
*\brief Multiply each matrix in the array by the left and the right matrix.
* For each page of this MatrixArray the operation performed is
@@ -229,6 +239,40 @@ class MatrixArray : public ValArray<T>
typename = std::enable_if_t<(std::is_same_v<T, std::complex<double>> && EnableBool)>>
MatrixArray<T> HermitianTranspose() const;
/**
* \brief Function that copies the current 1-page matrix into a new matrix with n copies of the
* original matrix
* \param nCopies Number of copies to make of the input 1-page MatrixArray
* \return Returns a new matrix with n page copies
*/
MatrixArray<T> MakeNCopies(size_t nCopies) const;
/**
* \brief Function extracts a page from a MatrixArray
* \param page Index of the page to be extracted from the n-page MatrixArray
* \return Returns an extracted page of the original MatrixArray
*/
MatrixArray<T> ExtractPage(size_t page) const;
/**
* \brief Function joins multiple pages into a single MatrixArray
* \param pages Vector containing n 1-page MatrixArrays to be merged into a n-page MatrixArray
* \return Returns a joint MatrixArray
*/
static MatrixArray<T> JoinPages(const std::vector<MatrixArray<T>>& pages);
/**
* \brief Function produces an identity MatrixArray with the specified size
* \param size Size of the output identity matrix
* \param pages Number of copies of the identity matrix
* \return Returns an identity MatrixArray
*/
static MatrixArray<T> IdentityMatrix(const size_t size, const size_t pages = 1);
/**
* \brief Function produces an identity MatrixArray with the same size
* as the input MatrixArray
* \param likeme Input matrix used to determine the size of the identity matrix
* \return Returns an identity MatrixArray
*/
static MatrixArray<T> IdentityMatrix(const MatrixArray& likeme);
protected:
// To simplify functions in MatrixArray that are using members from the template base class
using ValArray<T>::m_numRows;

169
src/core/test/matrix-array-test-suite.cc
View File

@@ -40,6 +40,27 @@ namespace tests
NS_LOG_COMPONENT_DEFINE("MatrixArrayTest");
/**
* \brief Function casts an input valArray "in" (type IN) to an output valArray "out" (type T)
* \param in Input valarray to be casted
* \param out Output valarray to receive casted values
*/
template <typename IN, typename T>
void
CastStdValarray(const std::valarray<IN>& in, std::valarray<T>& out)
{
// Ensure output valarray is the right size
if (out.size() != in.size())
{
out.resize(in.size());
}
// Perform the cast operation
std::transform(std::begin(in), std::end(in), std::begin(out), [](IN i) {
return static_cast<T>(i);
});
}
/**
* \ingroup matrixArray-tests
* MatrixArray test case for testing constructors, operators and other functions
@@ -443,6 +464,154 @@ MatrixArrayTestCase<T>::DoRun()
MatrixArray<T> m27 = MatrixArray<T>(std::move(jCasted));
NS_LOG_INFO("m27.GetSize ()" << m27.GetSize());
NS_TEST_ASSERT_MSG_EQ(jCastedSize, m27.GetSize(), "The number of elements are not equal.");
// test determinant
{
std::vector<std::pair<std::valarray<int>, T>> detTestCases{
// right-wraparound
{{1, 0, 7, 4, 2, 0, 6, 5, 3}, 62},
// left-wraparound
{{1, 4, 6, 0, 2, 5, 7, 0, 3}, 62},
// identity rank 3
{{1, 0, 0, 0, 1, 0, 0, 0, 1}, 1},
// identity rank 4
{{1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1}, 1},
// permutation matrix rank 4
{{0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0}, -1},
// positive det matrix rank 2
{{36, -5, -5, 43}, 1523},
// single value matrix rank 1
{{1}, 1},
};
for (const auto& [detVal, detRef] : detTestCases)
{
std::valarray<T> detCast(detVal.size());
CastStdValarray(detVal, detCast);
auto side = sqrt(detVal.size());
MatrixArray<T> detMat = MatrixArray<T>(side, side, std::move(detCast));
NS_TEST_ASSERT_MSG_EQ(detMat.Determinant()(0),
static_cast<T>(detRef),
"The determinants are not equal.");
}
}
// clang-format off
std::valarray<int> multiPageMatrixValues{
// page 0: identity matrix
1, 0, 0,
0, 1, 0,
0, 0, 1,
// page 1: permutation matrix
0, 0, 1,
0, 1, 0,
1, 0, 0,
// page 2: random matrix
1, 4, 6,
0, 2, 5,
7, 0, 3,
// page 3: upper triangular
5, -9, 2,
0, -4, -10,
0, 0, -3,
// page 4: lower triangular
-7, 0, 0,
-1, -9, 0,
-5, 6, -5,
// page 5: zero diagonal
0, 1, 2,
1, 0, 1,
2, 1, 0,
// page 6: zero bidiagonal
0, 1, 0,
1, 0, 1,
0, 1, 0
};
// clang-format on
std::valarray<T> castMultiPageMatrixValues(multiPageMatrixValues.size());
CastStdValarray(multiPageMatrixValues, castMultiPageMatrixValues);
MatrixArray<T> multiPageMatrix = MatrixArray<T>(3,
3,
multiPageMatrixValues.size() / 9,
std::move(castMultiPageMatrixValues));
// test the determinant in a multi-page MatrixArray
std::vector<int> determinants = {1, -1, 62, 60, -315, 4, 0};
auto det = multiPageMatrix.Determinant();
for (size_t page = 0; page < multiPageMatrix.GetNumPages(); page++)
{
NS_TEST_ASSERT_MSG_EQ(det(page),
static_cast<T>(determinants[page]),
"The determinants from the page " << std::to_string(page)
<< " are not equal.");
}
// test Frobenius norm in a multi-page MatrixArray
std::vector<double> fnorms = {sqrt(3), sqrt(3), 11.8322, 15.3297, 14.7309, 3.4641, 2};
auto frob = multiPageMatrix.FrobeniusNorm();
for (size_t page = 0; page < multiPageMatrix.GetNumPages(); page++)
{
NS_TEST_ASSERT_MSG_EQ_TOL(std::abs(frob(page)),
std::abs(static_cast<T>(fnorms[page])),
0.0001,
"The Frobenius norm from the page " << std::to_string(page)
<< " are not equal.");
}
// test page copying
for (size_t noOfCopies = 1; noOfCopies < 4; noOfCopies++)
{
auto copies = m27.MakeNCopies(noOfCopies);
NS_TEST_ASSERT_MSG_EQ(copies.GetNumPages(),
noOfCopies,
"Creating " << std::to_string(noOfCopies) << " copies failed.");
NS_TEST_ASSERT_MSG_EQ(copies.GetNumRows(),
m27.GetNumRows(),
"The copy doesn't have the same number of rows as the original.");
NS_TEST_ASSERT_MSG_EQ(copies.GetNumCols(),
m27.GetNumCols(),
"The copy doesn't have the same number of columns as the original.");
for (size_t page = 0; page < copies.GetNumPages(); page++)
{
T diff{};
for (size_t row = 0; row < copies.GetNumRows(); row++)
{
for (size_t col = 0; col < copies.GetNumCols(); col++)
{
diff += m27(row, col, 0) - copies(row, col, page);
}
}
NS_TEST_ASSERT_MSG_EQ(diff, T{}, "Mismatch in copied values.");
}
}
// test page 1 and 0 extraction
std::vector<MatrixArray<T>> pages{multiPageMatrix.ExtractPage(1),
multiPageMatrix.ExtractPage(0)};
// test page 1 and 0 joining
auto jointPagesMatrix = MatrixArray<T>::JoinPages(pages);
NS_TEST_ASSERT_MSG_EQ(jointPagesMatrix.GetNumPages(), 2, "Mismatch in number of join pages.");
for (size_t page = 0; page < jointPagesMatrix.GetNumPages(); page++)
{
T diff{};
for (size_t row = 0; row < jointPagesMatrix.GetNumRows(); row++)
{
for (size_t col = 0; col < jointPagesMatrix.GetNumCols(); col++)
{
diff += multiPageMatrix(row, col, 1 - page) - jointPagesMatrix(row, col, page);
}
}
NS_TEST_ASSERT_MSG_EQ(diff, T{}, "Mismatching pages.");
}
// test identity matrix
auto identityRank3Reference = multiPageMatrix.ExtractPage(0);
auto identityRank3 = MatrixArray<T>::IdentityMatrix(3);
NS_TEST_ASSERT_MSG_EQ(identityRank3, identityRank3Reference, "Mismatch in identity matrices.");
identityRank3 = MatrixArray<T>::IdentityMatrix(3, 10).ExtractPage(9);
NS_TEST_ASSERT_MSG_EQ(identityRank3, identityRank3Reference, "Mismatch in identity matrices.");
identityRank3 = MatrixArray<T>::IdentityMatrix(identityRank3Reference);
NS_TEST_ASSERT_MSG_EQ(identityRank3, identityRank3Reference, "Mismatch in identity matrices.");
}
/**