Skip to content
Permalink

Comparing changes

Choose two branches to see what’s changed or to start a new pull request. If you need to, you can also or learn more about diff comparisons.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also . Learn more about diff comparisons here.
base repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected base ref is valid, then try again.
Loading
base: ac7a114a0752
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 599f202ef721
Choose a head ref
  • 6 commits
  • 3 files changed
  • 1 contributor

Commits on Oct 4, 2020

  1. Upper triangular matrix.

    pleroy committed Oct 4, 2020
    Copy the full SHA
    84dac9a View commit details
  2. Tests and a fancy indexing.

    pleroy committed Oct 4, 2020
    Copy the full SHA
    c81e27d View commit details
  3. Cleanup.

    pleroy committed Oct 4, 2020
    Copy the full SHA
    6808e64 View commit details
  4. After egg's review.

    pleroy committed Oct 4, 2020
    Copy the full SHA
    392c10d View commit details
  5. After egg's review.

    pleroy committed Oct 4, 2020
    Copy the full SHA
    5389192 View commit details
  6. Merge pull request #2748 from pleroy/NewProjection

    Add support for upper-triangular matrices
    pleroy authored Oct 4, 2020
    Copy the full SHA
    599f202 View commit details
Showing with 373 additions and 53 deletions.
  1. +74 −2 numerics/unbounded_arrays.hpp
  2. +174 −0 numerics/unbounded_arrays_body.hpp
  3. +125 −51 numerics/unbounded_arrays_test.cpp
76 changes: 74 additions & 2 deletions numerics/unbounded_arrays.hpp
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@


#pragma once

#include <initializer_list>
@@ -74,7 +74,7 @@ class UnboundedLowerTriangularMatrix final {

bool operator==(UnboundedLowerTriangularMatrix const& right) const;

// For 0 < j <= i < rows, the entry a_ij is accessed as |a[i][j]|.
// For 0 j i < rows, the entry a_ij is accessed as |a[i][j]|.
// if i and j do not satisfy these conditions, the expression |a[i][j]| is
// erroneous.
Scalar* operator[](int index);
@@ -90,6 +90,73 @@ class UnboundedLowerTriangularMatrix final {
UnboundedLowerTriangularMatrix<S> const& matrix);
};

template<typename Scalar>
class UnboundedUpperTriangularMatrix final {
public:
explicit UnboundedUpperTriangularMatrix(int columns);
UnboundedUpperTriangularMatrix(int columns, uninitialized_t);

// The |data| must be in row-major format.
UnboundedUpperTriangularMatrix(std::initializer_list<Scalar> const& data);

void Extend(int extra_columns);
void Extend(int extra_columns, uninitialized_t);

// The |data| must be in row-major format.
void Extend(std::initializer_list<Scalar> const& data);

void EraseToEnd(int begin_column_index);

int columns() const;
int dimension() const;

bool operator==(UnboundedUpperTriangularMatrix const& right) const;

// A helper class for indexing column-major data in a human-friendly manner.
template<typename Matrix>
class Row {
public:
Scalar& operator[](int column);
Scalar const& operator[](int column) const;

private:
explicit Row(Matrix& matrix, int row);

Matrix& matrix_;
int row_;

template<typename S>
friend class UnboundedUpperTriangularMatrix;
};

// For 0 ≤ i ≤ j < columns, the entry a_ij is accessed as |a[i][j]|.
// if i and j do not satisfy these conditions, the expression |a[i][j]| is
// erroneous.
Row<UnboundedUpperTriangularMatrix> operator[](int row);
Row<UnboundedUpperTriangularMatrix const> operator[](int row) const;

private:
// For ease of writing matrices in tests, the input data is received in row-
// major format. This translates a trapezoidal slice to make it column-major.
static std::vector<Scalar, uninitialized_allocator<Scalar>> Transpose(
std::initializer_list<Scalar> const& data,
int current_columns,
int extra_columns);

int columns_;
// Stored in column-major format, so the data passed the public API must be
// transposed.
std::vector<Scalar, uninitialized_allocator<Scalar>> data_;

template<typename S>
friend std::ostream& operator<<(
std::ostream& out,
UnboundedUpperTriangularMatrix<S> const& matrix);

template<typename R>
friend class Row;
};

template<typename Scalar>
std::ostream& operator<<(std::ostream& out,
UnboundedVector<Scalar> const& vector);
@@ -98,9 +165,14 @@ template<typename Scalar>
std::ostream& operator<<(std::ostream& out,
UnboundedLowerTriangularMatrix<Scalar> const& matrix);

template<typename Scalar>
std::ostream& operator<<(std::ostream& out,
UnboundedUpperTriangularMatrix<Scalar> const& matrix);

} // namespace internal_unbounded_arrays

using internal_unbounded_arrays::UnboundedLowerTriangularMatrix;
using internal_unbounded_arrays::UnboundedUpperTriangularMatrix;
using internal_unbounded_arrays::UnboundedVector;

} // namespace numerics
174 changes: 174 additions & 0 deletions numerics/unbounded_arrays_body.hpp
Original file line number Diff line number Diff line change
@@ -4,6 +4,7 @@
#include "numerics/unbounded_arrays.hpp"

#include <cmath>
#include <vector>

#include "base/macros.hpp"
#include "quantities/elementary_functions.hpp"
@@ -157,6 +158,165 @@ Scalar const* UnboundedLowerTriangularMatrix<Scalar>::operator[](
return &data_[index * (index + 1) / 2];
}

template<typename Scalar>
UnboundedUpperTriangularMatrix<Scalar>::UnboundedUpperTriangularMatrix(
int const columns)
: columns_(columns),
data_(columns_ * (columns_ + 1) / 2, Scalar{}) {}

template<typename Scalar>
UnboundedUpperTriangularMatrix<Scalar>::UnboundedUpperTriangularMatrix(
int const columns,
uninitialized_t)
: columns_(columns),
data_(columns_ * (columns_ + 1) / 2) {}

template<typename Scalar>
UnboundedUpperTriangularMatrix<Scalar>::UnboundedUpperTriangularMatrix(
std::initializer_list<Scalar> const& data)
: columns_(
static_cast<int>(std::lround((-1 + Sqrt(8 * data.size())) * 0.5))),
data_(Transpose(data,
/*current_columns=*/0,
/*extra_columns=*/columns_)) {
DCHECK_EQ(data_.size(), columns_ * (columns_ + 1) / 2);
}

template<typename Scalar>
void UnboundedUpperTriangularMatrix<Scalar>::Extend(int const extra_columns) {
columns_ += extra_columns;
data_.resize(columns_ * (columns_ + 1) / 2, Scalar{});
}

template<typename Scalar>
void UnboundedUpperTriangularMatrix<Scalar>::Extend(int const extra_columns,
uninitialized_t) {
columns_ += extra_columns;
data_.resize(columns_ * (columns_ + 1) / 2);
}

template<typename Scalar>
void UnboundedUpperTriangularMatrix<Scalar>::Extend(
std::initializer_list<Scalar> const& data) {
int const new_columns = static_cast<int>(
std::lround((-1 + Sqrt(8 * (data_.size() + data.size()))) * 0.5));
auto transposed_data = Transpose(data,
/*current_columns=*/columns_,
/*extra_columns=*/new_columns - columns_);
columns_ = new_columns;
std::move(transposed_data.begin(),
transposed_data.end(),
std::back_inserter(data_));
DCHECK_EQ(data_.size(), columns_ * (columns_ + 1) / 2);
}

template<typename Scalar>
void UnboundedUpperTriangularMatrix<Scalar>::EraseToEnd(
int const begin_column_index) {
columns_ = begin_column_index;
data_.erase(data_.begin() + begin_column_index * (begin_column_index + 1) / 2,
data_.end());
}

template<typename Scalar>
int UnboundedUpperTriangularMatrix<Scalar>::columns() const {
return columns_;
}

template<typename Scalar>
int UnboundedUpperTriangularMatrix<Scalar>::dimension() const {
return data_.size();
}

template<typename Scalar>
bool UnboundedUpperTriangularMatrix<Scalar>::operator==(
UnboundedUpperTriangularMatrix const& right) const {
return columns_ == right.columns_ && data_ == right.data_;
}

template<typename Scalar>
template<typename Matrix>
Scalar& UnboundedUpperTriangularMatrix<Scalar>::Row<Matrix>::operator[](
int const column) {
DCHECK_LT(column, matrix_.columns_);
return matrix_.data_[column * (column + 1) / 2 + row_];
}

template<typename Scalar>
template<typename Matrix>
Scalar const&
UnboundedUpperTriangularMatrix<Scalar>::Row<Matrix>::operator[](
int const column) const {
DCHECK_LT(column, matrix_.columns_);
return matrix_.data_[column * (column + 1) / 2 + row_];
}

template<typename Scalar>
template<typename Matrix>
UnboundedUpperTriangularMatrix<Scalar>::Row<Matrix>::Row(Matrix& matrix,
int const row)
: matrix_(matrix),
row_(row) {}

template<typename Scalar>
auto UnboundedUpperTriangularMatrix<Scalar>::operator[](int const row)
-> Row<UnboundedUpperTriangularMatrix> {
return Row<UnboundedUpperTriangularMatrix>{*this, row};
}

template<typename Scalar>
auto UnboundedUpperTriangularMatrix<Scalar>::operator[](int const row) const
-> Row<UnboundedUpperTriangularMatrix const> {
return Row<UnboundedUpperTriangularMatrix const>{*this, row};
}

template<typename Scalar>
auto
UnboundedUpperTriangularMatrix<Scalar>::Transpose(
std::initializer_list<Scalar> const& data,
int const current_columns,
int const extra_columns) ->
std::vector<Scalar, uninitialized_allocator<Scalar>> {
// |data| is a trapezoidal slice at the end of the matrix. This is
// inconvenient to index, so we start by constructing a rectangular array with
// |extra_columns| columns and |current_columns + extra_columns| rows padded
// with junk.
std::vector<Scalar, uninitialized_allocator<Scalar>> padded;
{
padded.reserve(2 * data.size()); // An overestimate.
int row = 0;
int column = 0;
for (auto it = data.begin(); it != data.end();) {
if (row <= current_columns + column) {
padded.push_back(*it);
++it;
} else {
padded.emplace_back();
}
++column;
if (column == extra_columns) {
column = 0;
++row;
}
}
}

// Scan the padded array by column and append the part above the diagonal to
// the result.
std::vector<Scalar, uninitialized_allocator<Scalar>> result;
result.reserve(data.size());
int const number_of_rows = current_columns + extra_columns;
for (int column = 0; column < extra_columns; ++column) {
for (int row = 0; row < number_of_rows; ++row) {
if (row <= current_columns + column) {
result.push_back(padded[row * extra_columns + column]);
}
}
}

return result;
}

template<typename Scalar>
std::ostream& operator<<(std::ostream& out,
UnboundedVector<Scalar> const& vector) {
@@ -183,6 +343,20 @@ std::ostream& operator<<(std::ostream& out,
return out;
}

template<typename Scalar>
std::ostream& operator<<(std::ostream& out,
UnboundedUpperTriangularMatrix<Scalar> const& matrix) {
std::stringstream s;
// TODO(phl): Triangular printout.
s << "columns: " << matrix.columns_ << " ";
for (int i = 0; i < matrix.data_.size(); ++i) {
s << (i == 0 ? "{" : "") << matrix.data_[i]
<< (i == matrix.data_.size() - 1 ? "}" : ", ");
}
out << s.str();
return out;
}

} // namespace internal_unbounded_arrays
} // namespace numerics
} // namespace principia
Loading