diff --git a/src/cppmat/fix_diagonal_matrix.cpp b/src/cppmat/fix_diagonal_matrix.cpp index 322f1c2..ac67e71 100644 --- a/src/cppmat/fix_diagonal_matrix.cpp +++ b/src/cppmat/fix_diagonal_matrix.cpp @@ -221,7 +221,7 @@ matrix matrix::CopyDense(Iterator first, Iterator last) template inline -std::vector matrix::asVector() const +matrix::operator std::vector () const { std::vector out(mSize); @@ -230,6 +230,17 @@ std::vector matrix::asVector() const return out; } +// ================================================================================================= +// modify bounds check +// ================================================================================================= + +template +inline +void matrix::setPeriodic(bool periodic) +{ + mPeriodic = periodic; +} + // ================================================================================================= // get dimensions // ================================================================================================= @@ -320,7 +331,44 @@ const X& matrix::operator[](size_t i) const template inline -X& matrix::operator()(size_t a, size_t b) +X& matrix::operator()(int a, int b) +{ + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + assert( ( b < n && b >= -n ) or mPeriodic ); + + if ( a != b ) return mZero[0]; + + size_t A = static_cast( (n+(a%n)) % n ); + + return mData[A]; +} + +// ------------------------------------------------------------------------------------------------- + +template +inline +const X& matrix::operator()(int a, int b) const +{ + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + assert( ( b < n && b >= -n ) or mPeriodic ); + + if ( a != b ) return mZero[0]; + + size_t A = static_cast( (n+(a%n)) % n ); + + return mData[A]; +} + +// ------------------------------------------------------------------------------------------------- + +template +template +inline +X& matrix::operator()(T a, T b) { assert( a < N ); assert( b < N ); @@ -332,8 +380,9 @@ X& matrix::operator()(size_t a, size_t b) // ------------------------------------------------------------------------------------------------- template +template inline -const X& matrix::operator()(size_t a, size_t b) const +const X& matrix::operator()(T a, T b) const { assert( a < N ); assert( b < N ); @@ -348,7 +397,25 @@ const X& matrix::operator()(size_t a, size_t b) const template inline -size_t matrix::compress(size_t a, size_t b) const +size_t matrix::compress(int a, int b) const +{ + assert( a == b ); + + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + + return A; +} + +// ------------------------------------------------------------------------------------------------- + +template +template +inline +size_t matrix::compress(T a, T b) const { assert( a < N ); assert( b < N ); @@ -465,7 +532,42 @@ auto matrix::index(size_t i) const template inline -auto matrix::item(size_t a, size_t b) +auto matrix::item(int a, int b) +{ + assert( a == b ); + + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + + return this->begin() + A; +} + +// ------------------------------------------------------------------------------------------------- + +template +inline +auto matrix::item(int a, int b) const +{ + assert( a == b ); + + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + + return this->begin() + A; +} + +// ------------------------------------------------------------------------------------------------- + +template +template +inline +auto matrix::item(T a, T b) { assert( a < N ); assert( b < N ); @@ -477,8 +579,9 @@ auto matrix::item(size_t a, size_t b) // ------------------------------------------------------------------------------------------------- template +template inline -auto matrix::item(size_t a, size_t b) const +auto matrix::item(T a, T b) const { assert( a < N ); assert( b < N ); diff --git a/src/cppmat/fix_diagonal_matrix.h b/src/cppmat/fix_diagonal_matrix.h index c82b68d..180af85 100644 --- a/src/cppmat/fix_diagonal_matrix.h +++ b/src/cppmat/fix_diagonal_matrix.h @@ -28,10 +28,11 @@ class matrix protected: - static const size_t mSize=N; // total size == data.size() - static const size_t mRank=2; // rank (number of axes) - X mData[N]; // data container - X mZero[1]; // pointer to a zero entry + static const size_t mSize=N; // total size == data.size() + static const size_t mRank=2; // rank (number of axes) + X mData[N]; // data container + X mZero[1]; // pointer to a zero entry + bool mPeriodic=false; // if true: disable bounds-check where possible public: @@ -65,7 +66,10 @@ class matrix template static matrix CopyDense(Itr first, Itr last); // return plain storage as vector - std::vector asVector() const; + operator std::vector () const; + + // modify bounds-checks + void setPeriodic(bool periodic); // get dimensions size_t size() const; @@ -79,11 +83,22 @@ class matrix const X& operator[](size_t i) const; // index operators: access using matrix-indices - X& operator()(size_t a, size_t b); - const X& operator()(size_t a, size_t b) const; + X& operator()(int a, int b); + const X& operator()(int a, int b) const; + + // index operators: access using matrix-indices + template::value,void>::type> + X& operator()(T a, T b); + + template::value,void>::type> + const X& operator()(T a, T b) const; // index operators: matrix-indices -> plain storage (a,b -> i) - size_t compress(size_t a, size_t b) const; + size_t compress(int a, int b) const; + + // index operators: matrix-indices -> plain storage (a,b -> i) + template::value,void>::type> + size_t compress(T a, T b) const; // index operators: plain storage -> matrix-indices (i -> a,b) std::vector decompress(size_t i) const; @@ -103,8 +118,15 @@ class matrix auto index(size_t i) const; // iterator to specific entry: access using matrix-indices - auto item(size_t a, size_t b); - auto item(size_t a, size_t b) const; + auto item(int a, int b); + auto item(int a, int b) const; + + // iterator to specific entry: access using matrix-indices + template::value,void>::type> + auto item(T a, T b); + + template::value,void>::type> + auto item(T a, T b) const; // initialization void setRandom(X lower=(X)0, X upper=(X)1); diff --git a/src/cppmat/fix_regular_array.cpp b/src/cppmat/fix_regular_array.cpp index 682f6df..a6950e2 100644 --- a/src/cppmat/fix_regular_array.cpp +++ b/src/cppmat/fix_regular_array.cpp @@ -224,6 +224,17 @@ array::operator std::vector () const return out; } +// ================================================================================================= +// modify bounds check +// ================================================================================================= + +template +inline +void array::setPeriodic(bool periodic) +{ + mPeriodic = periodic; +} + // ================================================================================================= // get dimensions // ================================================================================================= @@ -341,7 +352,7 @@ X& array::operator()(int a) { int na = static_cast(mShape[0]); - assert( a < na && a >= -na ); + assert( ( a < na && a >= -na ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); @@ -357,7 +368,7 @@ const X& array::operator()(int a) const { int na = static_cast(mShape[0]); - assert( a < na && a >= -na ); + assert( ( a < na && a >= -na ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); @@ -374,8 +385,8 @@ X& array::operator()(int a, int b) int na = static_cast(mShape[0]); int nb = static_cast(mShape[1]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -394,8 +405,8 @@ const X& array::operator()(int a, int b) const int na = static_cast(mShape[0]); int nb = static_cast(mShape[1]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -415,9 +426,9 @@ X& array::operator()(int a, int b, int c) int nb = static_cast(mShape[1]); int nc = static_cast(mShape[2]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -439,9 +450,9 @@ const X& array::operator()(int a, int b, int c) const int nb = static_cast(mShape[1]); int nc = static_cast(mShape[2]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -464,10 +475,10 @@ X& array::operator()(int a, int b, int c, int d) int nc = static_cast(mShape[2]); int nd = static_cast(mShape[3]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -492,10 +503,10 @@ const X& array::operator()(int a, int b, int c, int d) const int nc = static_cast(mShape[2]); int nd = static_cast(mShape[3]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -521,11 +532,11 @@ X& array::operator()(int a, int b, int c, int d, int e) int nd = static_cast(mShape[3]); int ne = static_cast(mShape[4]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -553,11 +564,11 @@ const X& array::operator()(int a, int b, int c, int d, int e int nd = static_cast(mShape[3]); int ne = static_cast(mShape[4]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -586,12 +597,12 @@ X& array::operator()(int a, int b, int c, int d, int e, int int ne = static_cast(mShape[4]); int nf = static_cast(mShape[5]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); - assert( f < nf && f >= -nf ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); + assert( ( f < nf && f >= -nf ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -622,12 +633,12 @@ const X& array::operator()(int a, int b, int c, int d, int e int ne = static_cast(mShape[4]); int nf = static_cast(mShape[5]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); - assert( f < nf && f >= -nf ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); + assert( ( f < nf && f >= -nf ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -939,7 +950,7 @@ size_t array::compress(int a) const { int na = static_cast(mShape[0]); - assert( a < na && a >= -na ); + assert( ( a < na && a >= -na ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); @@ -955,8 +966,8 @@ size_t array::compress(int a, int b) const int na = static_cast(mShape[0]); int nb = static_cast(mShape[1]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -975,9 +986,9 @@ size_t array::compress(int a, int b, int c) const int nb = static_cast(mShape[1]); int nc = static_cast(mShape[2]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -999,10 +1010,10 @@ size_t array::compress(int a, int b, int c, int d) const int nc = static_cast(mShape[2]); int nd = static_cast(mShape[3]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1027,11 +1038,11 @@ size_t array::compress(int a, int b, int c, int d, int e) co int nd = static_cast(mShape[3]); int ne = static_cast(mShape[4]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1059,12 +1070,12 @@ size_t array::compress(int a, int b, int c, int d, int e, in int ne = static_cast(mShape[4]); int nf = static_cast(mShape[5]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); - assert( f < nf && f >= -nf ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); + assert( ( f < nf && f >= -nf ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1298,7 +1309,7 @@ auto array::item(int a) { int na = static_cast(mShape[0]); - assert( a < na && a >= -na ); + assert( ( a < na && a >= -na ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); @@ -1314,7 +1325,7 @@ auto array::item(int a) const { int na = static_cast(mShape[0]); - assert( a < na && a >= -na ); + assert( ( a < na && a >= -na ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); @@ -1331,8 +1342,8 @@ auto array::item(int a, int b) int na = static_cast(mShape[0]); int nb = static_cast(mShape[1]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1351,8 +1362,8 @@ auto array::item(int a, int b) const int na = static_cast(mShape[0]); int nb = static_cast(mShape[1]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1372,9 +1383,9 @@ auto array::item(int a, int b, int c) int nb = static_cast(mShape[1]); int nc = static_cast(mShape[2]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1396,9 +1407,9 @@ auto array::item(int a, int b, int c) const int nb = static_cast(mShape[1]); int nc = static_cast(mShape[2]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1421,10 +1432,10 @@ auto array::item(int a, int b, int c, int d) int nc = static_cast(mShape[2]); int nd = static_cast(mShape[3]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1449,10 +1460,10 @@ auto array::item(int a, int b, int c, int d) const int nc = static_cast(mShape[2]); int nd = static_cast(mShape[3]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1478,11 +1489,11 @@ auto array::item(int a, int b, int c, int d, int e) int nd = static_cast(mShape[3]); int ne = static_cast(mShape[4]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1510,11 +1521,11 @@ auto array::item(int a, int b, int c, int d, int e) const int nd = static_cast(mShape[3]); int ne = static_cast(mShape[4]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1543,12 +1554,12 @@ auto array::item(int a, int b, int c, int d, int e, int f) int ne = static_cast(mShape[4]); int nf = static_cast(mShape[5]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); - assert( f < nf && f >= -nf ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); + assert( ( f < nf && f >= -nf ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1579,12 +1590,12 @@ auto array::item(int a, int b, int c, int d, int e, int f) c int ne = static_cast(mShape[4]); int nf = static_cast(mShape[5]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); - assert( f < nf && f >= -nf ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); + assert( ( f < nf && f >= -nf ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); diff --git a/src/cppmat/fix_regular_array.h b/src/cppmat/fix_regular_array.h index 29baac7..700b773 100644 --- a/src/cppmat/fix_regular_array.h +++ b/src/cppmat/fix_regular_array.h @@ -32,12 +32,13 @@ class array protected: - static const size_t MAX_DIM=6; // maximum number of dimensions - static const size_t mSize=I*J*K*L*M*N; // total size == data.size() == prod(shape) - static const size_t mRank=RANK; // rank (number of axes) - size_t mShape [MAX_DIM]; // number of entries along each axis - size_t mStrides[MAX_DIM]; // stride length for each index - X mData[I*J*K*L*M*N]; // data container + static const size_t MAX_DIM=6; // maximum number of dimensions + static const size_t mSize=I*J*K*L*M*N; // total size == data.size() == prod(shape) + static const size_t mRank=RANK; // rank (number of axes) + size_t mShape [MAX_DIM]; // number of entries along each axis + size_t mStrides[MAX_DIM]; // stride length for each index + X mData[I*J*K*L*M*N]; // data container + bool mPeriodic=false; // if true: disable bounds-check where possible public: @@ -71,6 +72,9 @@ class array // return plain storage as vector operator std::vector () const; + // modify bounds-checks + void setPeriodic(bool periodic); + // get dimensions size_t size() const; size_t rank() const; diff --git a/src/cppmat/fix_symmetric_matrix.cpp b/src/cppmat/fix_symmetric_matrix.cpp index 89a8d2c..c746a50 100644 --- a/src/cppmat/fix_symmetric_matrix.cpp +++ b/src/cppmat/fix_symmetric_matrix.cpp @@ -231,7 +231,7 @@ matrix matrix::CopyDense(Iterator first, Iterator last) template inline -std::vector matrix::asVector() const +matrix::operator std::vector () const { std::vector out(mSize); @@ -240,6 +240,17 @@ std::vector matrix::asVector() const return out; } +// ================================================================================================= +// modify bounds check +// ================================================================================================= + +template +inline +void matrix::setPeriodic(bool periodic) +{ + mPeriodic = periodic; +} + // ================================================================================================= // get dimensions // ================================================================================================= @@ -330,7 +341,44 @@ const X& matrix::operator[](size_t i) const template inline -X& matrix::operator()(size_t a, size_t b) +X& matrix::operator()(int a, int b) +{ + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + assert( ( b < n && b >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + size_t B = static_cast( (n+(b%n)) % n ); + + if (A <= B) return mData[ A*N - (A-1)*A/2 + B - A ]; + else return mData[ B*N - (B-1)*B/2 + A - B ]; +} + +// ------------------------------------------------------------------------------------------------- + +template +inline +const X& matrix::operator()(int a, int b) const +{ + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + assert( ( b < n && b >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + size_t B = static_cast( (n+(b%n)) % n ); + + if (A <= B) return mData[ A*N - (A-1)*A/2 + B - A ]; + else return mData[ B*N - (B-1)*B/2 + A - B ]; +} + +// ------------------------------------------------------------------------------------------------- + +template +template +inline +X& matrix::operator()(T a, T b) { assert( a < N ); assert( b < N ); @@ -342,8 +390,9 @@ X& matrix::operator()(size_t a, size_t b) // ------------------------------------------------------------------------------------------------- template +template inline -const X& matrix::operator()(size_t a, size_t b) const +const X& matrix::operator()(T a, T b) const { assert( a < N ); assert( b < N ); @@ -358,7 +407,26 @@ const X& matrix::operator()(size_t a, size_t b) const template inline -size_t matrix::compress(size_t a, size_t b) const +size_t matrix::compress(int a, int b) const +{ + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + assert( ( b < n && b >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + size_t B = static_cast( (n+(b%n)) % n ); + + if (A <= B) return A*N - (A-1)*A/2 + B - A; + else return B*N - (B-1)*B/2 + A - B; +} + +// ------------------------------------------------------------------------------------------------- + +template +template +inline +size_t matrix::compress(T a, T b) const { assert( a < N ); assert( b < N ); @@ -483,7 +551,44 @@ auto matrix::index(size_t i) const template inline -auto matrix::item(size_t a, size_t b) +auto matrix::item(int a, int b) +{ + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + assert( ( b < n && b >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + size_t B = static_cast( (n+(b%n)) % n ); + + if (A <= B) return begin() + ( A*N - (A-1)*A/2 + B - A ); + else return begin() + ( B*N - (B-1)*B/2 + A - B ); +} + +// ------------------------------------------------------------------------------------------------- + +template +inline +auto matrix::item(int a, int b) const +{ + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + assert( ( b < n && b >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + size_t B = static_cast( (n+(b%n)) % n ); + + if (A <= B) return begin() + ( A*N - (A-1)*A/2 + B - A ); + else return begin() + ( B*N - (B-1)*B/2 + A - B ); +} + +// ------------------------------------------------------------------------------------------------- + +template +template +inline +auto matrix::item(T a, T b) { assert( a < N ); assert( b < N ); @@ -495,8 +600,9 @@ auto matrix::item(size_t a, size_t b) // ------------------------------------------------------------------------------------------------- template +template inline -auto matrix::item(size_t a, size_t b) const +auto matrix::item(T a, T b) const { assert( a < N ); assert( b < N ); diff --git a/src/cppmat/fix_symmetric_matrix.h b/src/cppmat/fix_symmetric_matrix.h index baa2487..628c6ae 100644 --- a/src/cppmat/fix_symmetric_matrix.h +++ b/src/cppmat/fix_symmetric_matrix.h @@ -31,6 +31,7 @@ class matrix static const size_t mSize=(N+1)*N/2; // total size == data.size() static const size_t mRank=2; // rank (number of axes) X mData[(N+1)*N/2]; // data container + bool mPeriodic=false; // if true: disable bounds-check where possible public: @@ -67,7 +68,10 @@ class matrix template static matrix CopyDense(Itr first, Itr last); // return plain storage as vector - std::vector asVector() const; + operator std::vector () const; + + // modify bounds-checks + void setPeriodic(bool periodic); // get dimensions size_t size() const; @@ -81,11 +85,22 @@ class matrix const X& operator[](size_t i) const; // index operators: access using matrix-indices - X& operator()(size_t a, size_t b); - const X& operator()(size_t a, size_t b) const; + X& operator()(int a, int b); + const X& operator()(int a, int b) const; + + // index operators: access using matrix-indices + template::value,void>::type> + X& operator()(T a, T b); + + template::value,void>::type> + const X& operator()(T a, T b) const; // index operators: matrix-indices -> plain storage (a,b -> i) - size_t compress(size_t a, size_t b) const; + size_t compress(int a, int b) const; + + // index operators: matrix-indices -> plain storage (a,b -> i) + template::value,void>::type> + size_t compress(T a, T b) const; // index operators: plain storage -> matrix-indices (i -> a,b) std::vector decompress(size_t i) const; @@ -105,8 +120,15 @@ class matrix auto index(size_t i) const; // iterator to specific entry: access using matrix-indices - auto item(size_t a, size_t b); - auto item(size_t a, size_t b) const; + auto item(int a, int b); + auto item(int a, int b) const; + + // iterator to specific entry: access using matrix-indices + template::value,void>::type> + auto item(T a, T b); + + template::value,void>::type> + auto item(T a, T b) const; // initialization void setRandom(X lower=(X)0, X upper=(X)1); diff --git a/src/cppmat/map_diagonal_matrix.cpp b/src/cppmat/map_diagonal_matrix.cpp index f0b1bc8..5a2da6c 100644 --- a/src/cppmat/map_diagonal_matrix.cpp +++ b/src/cppmat/map_diagonal_matrix.cpp @@ -73,7 +73,7 @@ matrix matrix::Map(const X *D) template inline -std::vector matrix::asVector() const +matrix::operator std::vector () const { std::vector out(mSize); @@ -82,6 +82,17 @@ std::vector matrix::asVector() const return out; } +// ================================================================================================= +// modify bounds check +// ================================================================================================= + +template +inline +void matrix::setPeriodic(bool periodic) +{ + mPeriodic = periodic; +} + // ================================================================================================= // get dimensions // ================================================================================================= @@ -161,7 +172,26 @@ const X& matrix::operator[](size_t i) const template inline -const X& matrix::operator()(size_t a, size_t b) const +const X& matrix::operator()(int a, int b) const +{ + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + assert( ( b < n && b >= -n ) or mPeriodic ); + + if ( a != b ) return mZero[0]; + + size_t A = static_cast( (n+(a%n)) % n ); + + return mData[A]; +} + +// ------------------------------------------------------------------------------------------------- + +template +template +inline +const X& matrix::operator()(T a, T b) const { assert( a < N ); assert( b < N ); @@ -176,7 +206,25 @@ const X& matrix::operator()(size_t a, size_t b) const template inline -size_t matrix::compress(size_t a, size_t b) const +size_t matrix::compress(int a, int b) const +{ + assert( a == b ); + + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + + return A; +} + +// ------------------------------------------------------------------------------------------------- + +template +template +inline +size_t matrix::compress(T a, T b) const { assert( a < N ); assert( b < N ); @@ -255,7 +303,25 @@ auto matrix::index(size_t i) const template inline -auto matrix::item(size_t a, size_t b) const +auto matrix::item(int a, int b) const +{ + assert( a == b ); + + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + + return this->begin() + A; +} + +// ------------------------------------------------------------------------------------------------- + +template +template +inline +auto matrix::item(T a, T b) const { assert( a < N ); assert( b < N ); diff --git a/src/cppmat/map_diagonal_matrix.h b/src/cppmat/map_diagonal_matrix.h index f6a68ea..b561192 100644 --- a/src/cppmat/map_diagonal_matrix.h +++ b/src/cppmat/map_diagonal_matrix.h @@ -28,10 +28,11 @@ class matrix protected: - static const size_t mSize=N; // total size == data.size() - static const size_t mRank=2; // rank (number of axes) - const X *mData; // data container - X mZero[1]; // pointer to a zero entry + static const size_t mSize=N; // total size == data.size() + static const size_t mRank=2; // rank (number of axes) + const X *mData; // data container + X mZero[1]; // pointer to a zero entry + bool mPeriodic=false; // if true: disable bounds-check where possible public: @@ -48,7 +49,10 @@ class matrix static matrix Map(const X *D); // return plain storage as vector - std::vector asVector() const; + operator std::vector () const; + + // modify bounds-checks + void setPeriodic(bool periodic); // get dimensions size_t size() const; @@ -61,10 +65,18 @@ class matrix const X& operator[](size_t i) const; // index operators: access using matrix-indices - const X& operator()(size_t a, size_t b) const; + const X& operator()(int a, int b) const; + + // index operators: access using matrix-indices + template::value,void>::type> + const X& operator()(T a, T b) const; // index operators: matrix-indices -> plain storage (a,b -> i) - size_t compress(size_t a, size_t b) const; + size_t compress(int a, int b) const; + + // index operators: matrix-indices -> plain storage (a,b -> i) + template::value,void>::type> + size_t compress(T a, T b) const; // index operators: plain storage -> matrix-indices (i -> a,b) std::vector decompress(size_t i) const; @@ -80,7 +92,11 @@ class matrix auto index(size_t i) const; // iterator to specific entry: access using matrix-indices - auto item(size_t a, size_t b) const; + auto item(int a, int b) const; + + // iterator to specific entry: access using matrix-indices + template::value,void>::type> + auto item(T a, T b) const; // initialization void setMap(const X *D); diff --git a/src/cppmat/map_regular_array.cpp b/src/cppmat/map_regular_array.cpp index 34da176..fdccc8d 100644 --- a/src/cppmat/map_regular_array.cpp +++ b/src/cppmat/map_regular_array.cpp @@ -91,6 +91,17 @@ array::operator std::vector () const return out; } +// ================================================================================================= +// modify bounds check +// ================================================================================================= + +template +inline +void array::setPeriodic(bool periodic) +{ + mPeriodic = periodic; +} + // ================================================================================================= // get dimensions // ================================================================================================= @@ -197,7 +208,7 @@ const X& array::operator()(int a) const { int na = static_cast(mShape[0]); - assert( a < na && a >= -na ); + assert( ( a < na && a >= -na ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); @@ -214,8 +225,8 @@ const X& array::operator()(int a, int b) const int na = static_cast(mShape[0]); int nb = static_cast(mShape[1]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -235,9 +246,9 @@ const X& array::operator()(int a, int b, int c) const int nb = static_cast(mShape[1]); int nc = static_cast(mShape[2]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -260,10 +271,10 @@ const X& array::operator()(int a, int b, int c, int d) const int nc = static_cast(mShape[2]); int nd = static_cast(mShape[3]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -289,11 +300,11 @@ const X& array::operator()(int a, int b, int c, int d, int e int nd = static_cast(mShape[3]); int ne = static_cast(mShape[4]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -322,12 +333,12 @@ const X& array::operator()(int a, int b, int c, int d, int e int ne = static_cast(mShape[4]); int nf = static_cast(mShape[5]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); - assert( f < nf && f >= -nf ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); + assert( ( f < nf && f >= -nf ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -498,7 +509,7 @@ size_t array::compress(int a) const { int na = static_cast(mShape[0]); - assert( a < na && a >= -na ); + assert( ( a < na && a >= -na ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); @@ -514,8 +525,8 @@ size_t array::compress(int a, int b) const int na = static_cast(mShape[0]); int nb = static_cast(mShape[1]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -534,9 +545,9 @@ size_t array::compress(int a, int b, int c) const int nb = static_cast(mShape[1]); int nc = static_cast(mShape[2]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -558,10 +569,10 @@ size_t array::compress(int a, int b, int c, int d) const int nc = static_cast(mShape[2]); int nd = static_cast(mShape[3]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -586,11 +597,11 @@ size_t array::compress(int a, int b, int c, int d, int e) co int nd = static_cast(mShape[3]); int ne = static_cast(mShape[4]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -618,12 +629,12 @@ size_t array::compress(int a, int b, int c, int d, int e, in int ne = static_cast(mShape[4]); int nf = static_cast(mShape[5]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); - assert( f < nf && f >= -nf ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); + assert( ( f < nf && f >= -nf ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -819,7 +830,7 @@ auto array::item(int a) const { int na = static_cast(mShape[0]); - assert( a < na && a >= -na ); + assert( ( a < na && a >= -na ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); @@ -836,8 +847,8 @@ auto array::item(int a, int b) const int na = static_cast(mShape[0]); int nb = static_cast(mShape[1]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -857,9 +868,9 @@ auto array::item(int a, int b, int c) const int nb = static_cast(mShape[1]); int nc = static_cast(mShape[2]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -882,10 +893,10 @@ auto array::item(int a, int b, int c, int d) const int nc = static_cast(mShape[2]); int nd = static_cast(mShape[3]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -911,11 +922,11 @@ auto array::item(int a, int b, int c, int d, int e) const int nd = static_cast(mShape[3]); int ne = static_cast(mShape[4]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -944,12 +955,12 @@ auto array::item(int a, int b, int c, int d, int e, int f) c int ne = static_cast(mShape[4]); int nf = static_cast(mShape[5]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); - assert( f < nf && f >= -nf ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); + assert( ( f < nf && f >= -nf ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); diff --git a/src/cppmat/map_regular_array.h b/src/cppmat/map_regular_array.h index 188dd9a..0b4a662 100644 --- a/src/cppmat/map_regular_array.h +++ b/src/cppmat/map_regular_array.h @@ -25,12 +25,13 @@ class array { protected: - static const size_t MAX_DIM=6; // maximum number of dimensions - static const size_t mSize=I*J*K*L*M*N; // total size == data.size() == prod(shape) - static const size_t mRank=RANK; // rank (number of axes) - size_t mShape [MAX_DIM]; // number of entries along each axis - size_t mStrides[MAX_DIM]; // stride length for each index - const X *mData; // data container + static const size_t MAX_DIM=6; // maximum number of dimensions + static const size_t mSize=I*J*K*L*M*N; // total size == data.size() == prod(shape) + static const size_t mRank=RANK; // rank (number of axes) + size_t mShape [MAX_DIM]; // number of entries along each axis + size_t mStrides[MAX_DIM]; // stride length for each index + const X *mData; // data container + bool mPeriodic=false; // if true: disable bounds-check where possible public: @@ -49,6 +50,9 @@ class array // return plain storage as vector operator std::vector () const; + // modify bounds-checks + void setPeriodic(bool periodic); + // get dimensions size_t size() const; size_t rank() const; diff --git a/src/cppmat/map_symmetric_matrix.cpp b/src/cppmat/map_symmetric_matrix.cpp index 114ac27..e43bfcc 100644 --- a/src/cppmat/map_symmetric_matrix.cpp +++ b/src/cppmat/map_symmetric_matrix.cpp @@ -70,7 +70,7 @@ matrix matrix::Map(const X *D) template inline -std::vector matrix::asVector() const +matrix::operator std::vector () const { std::vector out(mSize); @@ -79,6 +79,17 @@ std::vector matrix::asVector() const return out; } +// ================================================================================================= +// modify bounds check +// ================================================================================================= + +template +inline +void matrix::setPeriodic(bool periodic) +{ + mPeriodic = periodic; +} + // ================================================================================================= // get dimensions // ================================================================================================= @@ -158,7 +169,26 @@ const X& matrix::operator[](size_t i) const template inline -const X& matrix::operator()(size_t a, size_t b) const +const X& matrix::operator()(int a, int b) const +{ + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + assert( ( b < n && b >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + size_t B = static_cast( (n+(b%n)) % n ); + + if (A <= B) return mData[ A*N - (A-1)*A/2 + B - A ]; + else return mData[ B*N - (B-1)*B/2 + A - B ]; +} + +// ------------------------------------------------------------------------------------------------- + +template +template +inline +const X& matrix::operator()(T a, T b) const { assert( a < N ); assert( b < N ); @@ -173,7 +203,26 @@ const X& matrix::operator()(size_t a, size_t b) const template inline -size_t matrix::compress(size_t a, size_t b) const +size_t matrix::compress(int a, int b) const +{ + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + assert( ( b < n && b >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + size_t B = static_cast( (n+(b%n)) % n ); + + if (A <= B) return A*N - (A-1)*A/2 + B - A; + else return B*N - (B-1)*B/2 + A - B; +} + +// ------------------------------------------------------------------------------------------------- + +template +template +inline +size_t matrix::compress(T a, T b) const { assert( a < N ); assert( b < N ); @@ -260,7 +309,26 @@ auto matrix::index(size_t i) const template inline -auto matrix::item(size_t a, size_t b) const +auto matrix::item(int a, int b) const +{ + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + assert( ( b < n && b >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + size_t B = static_cast( (n+(b%n)) % n ); + + if (A <= B) return begin() + ( A*N - (A-1)*A/2 + B - A ); + else return begin() + ( B*N - (B-1)*B/2 + A - B ); +} + +// ------------------------------------------------------------------------------------------------- + +template +template +inline +auto matrix::item(T a, T b) const { assert( a < N ); assert( b < N ); diff --git a/src/cppmat/map_symmetric_matrix.h b/src/cppmat/map_symmetric_matrix.h index 14e842c..4a2e49e 100644 --- a/src/cppmat/map_symmetric_matrix.h +++ b/src/cppmat/map_symmetric_matrix.h @@ -31,6 +31,7 @@ class matrix static const size_t mSize=(N+1)*N/2; // total size == data.size() static const size_t mRank=2; // rank (number of axes) const X *mData; // data container + bool mPeriodic=false; // if true: disable bounds-check where possible public: @@ -47,7 +48,10 @@ class matrix static matrix Map(const X *D); // return plain storage as vector - std::vector asVector() const; + operator std::vector () const; + + // modify bounds-checks + void setPeriodic(bool periodic); // get dimensions size_t size() const; @@ -60,10 +64,18 @@ class matrix const X& operator[](size_t i) const; // index operators: access using matrix-indices - const X& operator()(size_t a, size_t b) const; + const X& operator()(int a, int b) const; + + // index operators: access using matrix-indices + template::value,void>::type> + const X& operator()(T a, T b) const; // index operators: matrix-indices -> plain storage (a,b -> i) - size_t compress(size_t a, size_t b) const; + size_t compress(int a, int b) const; + + // index operators: matrix-indices -> plain storage (a,b -> i) + template::value,void>::type> + size_t compress(T a, T b) const; // index operators: plain storage -> matrix-indices (i -> a,b) std::vector decompress(size_t i) const; @@ -79,7 +91,11 @@ class matrix auto index(size_t i) const; // iterator to specific entry: access using matrix-indices - auto item(size_t a, size_t b) const; + auto item(int a, int b) const; + + // iterator to specific entry: access using matrix-indices + template::value,void>::type> + auto item(T a, T b) const; // initialization void setMap(const X *D); diff --git a/src/cppmat/stl.cpp b/src/cppmat/stl.cpp index 00b0271..bd47463 100644 --- a/src/cppmat/stl.cpp +++ b/src/cppmat/stl.cpp @@ -15,6 +15,8 @@ namespace cppmat { +// ================================================================================================= +// delete a specific item from a vector // ================================================================================================= template @@ -35,6 +37,8 @@ std::vector del(const std::vector &A, int idx) return B; } +// ================================================================================================= +// delete a specific item from a vector // ================================================================================================= template @@ -50,6 +54,8 @@ std::vector del(const std::vector &A, size_t idx) return B; } +// ================================================================================================= +// return the indices that would sort the vector // ================================================================================================= template @@ -77,6 +83,8 @@ std::vector argsort(const std::vector &v, bool ascending) return jdx; } +// ================================================================================================= +// join items to string // ================================================================================================= template @@ -97,7 +105,28 @@ std::string to_string(const std::vector &A, std::string join) } // namespace ... -// ------------------------------------------------------------------------------------------------- +// ================================================================================================= +// print operator +// ================================================================================================= + +#ifndef CPPMAT_NOSTD +template +inline +std::ostream& operator<<(std::ostream& out, const std::vector& src) +{ + auto w = out.width(); + auto p = out.precision(); + + for ( size_t j = 0 ; j < src.size() ; ++j ) { + out << std::setw(w) << std::setprecision(p) << src[j]; + if ( j != src.size()-1 ) out << ", "; + } + + return out; +} +#endif + +// ================================================================================================= #endif diff --git a/src/cppmat/stl.h b/src/cppmat/stl.h index d98a4d5..96ba268 100644 --- a/src/cppmat/stl.h +++ b/src/cppmat/stl.h @@ -21,13 +21,9 @@ namespace cppmat { template std::vector del(const std::vector &A, int idx); template std::vector del(const std::vector &A, size_t idx); -// ------------------------------------------------------------------------------------------------- - // return the indices that would sort the vector template std::vector argsort(const std::vector &v, bool ascending=true); -// ------------------------------------------------------------------------------------------------- - // convert vector items to string, and join these string together using the "join" string template std::string to_string(const std::vector &A, std::string join=", "); @@ -35,7 +31,14 @@ template std::string to_string(const std::vector &A, std::string joi } // namespace ... -// ------------------------------------------------------------------------------------------------- +// ================================================================================================= + +// print operator +#ifndef CPPMAT_NOSTD +template std::ostream& operator<<(std::ostream& out, const std::vector& src); +#endif + +// ================================================================================================= #endif diff --git a/src/cppmat/var_diagonal_matrix.cpp b/src/cppmat/var_diagonal_matrix.cpp index 1a2d0dd..9daff95 100644 --- a/src/cppmat/var_diagonal_matrix.cpp +++ b/src/cppmat/var_diagonal_matrix.cpp @@ -227,7 +227,7 @@ matrix matrix::CopyDense(size_t m, size_t n, Iterator first, Iterator last template inline -std::vector matrix::asVector() const +matrix::operator std::vector () const { return mData; } @@ -255,6 +255,17 @@ void matrix::resize(size_t m, size_t n) if ( mSize != size ) mData.resize(mSize); } +// ================================================================================================= +// modify bounds check +// ================================================================================================= + +template +inline +void matrix::setPeriodic(bool periodic) +{ + mPeriodic = periodic; +} + // ================================================================================================= // get dimensions // ================================================================================================= @@ -345,7 +356,44 @@ const X& matrix::operator[](size_t i) const template inline -X& matrix::operator()(size_t a, size_t b) +X& matrix::operator()(int a, int b) +{ + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + assert( ( b < n && b >= -n ) or mPeriodic ); + + if ( a != b ) return mZero[0]; + + size_t A = static_cast( (n+(a%n)) % n ); + + return mData[A]; +} + +// ------------------------------------------------------------------------------------------------- + +template +inline +const X& matrix::operator()(int a, int b) const +{ + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + assert( ( b < n && b >= -n ) or mPeriodic ); + + if ( a != b ) return mZero[0]; + + size_t A = static_cast( (n+(a%n)) % n ); + + return mData[A]; +} + +// ------------------------------------------------------------------------------------------------- + +template +template +inline +X& matrix::operator()(T a, T b) { assert( a < N ); assert( b < N ); @@ -357,8 +405,9 @@ X& matrix::operator()(size_t a, size_t b) // ------------------------------------------------------------------------------------------------- template +template inline -const X& matrix::operator()(size_t a, size_t b) const +const X& matrix::operator()(T a, T b) const { assert( a < N ); assert( b < N ); @@ -373,7 +422,25 @@ const X& matrix::operator()(size_t a, size_t b) const template inline -size_t matrix::compress(size_t a, size_t b) const +size_t matrix::compress(int a, int b) const +{ + assert( a == b ); + + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + + return A; +} + +// ------------------------------------------------------------------------------------------------- + +template +template +inline +size_t matrix::compress(T a, T b) const { assert( a < N ); assert( b < N ); @@ -490,7 +557,42 @@ auto matrix::index(size_t i) const template inline -auto matrix::item(size_t a, size_t b) +auto matrix::item(int a, int b) +{ + assert( a == b ); + + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + + return this->begin() + A; +} + +// ------------------------------------------------------------------------------------------------- + +template +inline +auto matrix::item(int a, int b) const +{ + assert( a == b ); + + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + + return this->begin() + A; +} + +// ------------------------------------------------------------------------------------------------- + +template +template +inline +auto matrix::item(T a, T b) { assert( a < N ); assert( b < N ); @@ -502,8 +604,9 @@ auto matrix::item(size_t a, size_t b) // ------------------------------------------------------------------------------------------------- template +template inline -auto matrix::item(size_t a, size_t b) const +auto matrix::item(T a, T b) const { assert( a < N ); assert( b < N ); diff --git a/src/cppmat/var_diagonal_matrix.h b/src/cppmat/var_diagonal_matrix.h index ee7005c..c585669 100644 --- a/src/cppmat/var_diagonal_matrix.h +++ b/src/cppmat/var_diagonal_matrix.h @@ -25,11 +25,12 @@ class matrix { protected: - size_t mSize=0; // total size == data.size() - static const size_t mRank=2; // rank (number of axes) - size_t N=0; // number of rows/columns - std::vector mData; // data container - X mZero[1]; // pointer to a zero entry + size_t mSize=0; // total size == data.size() + static const size_t mRank=2; // rank (number of axes) + size_t N=0; // number of rows/columns + std::vector mData; // data container + X mZero[1]; // pointer to a zero entry + bool mPeriodic=false; // if true: disable bounds-check where possible public: @@ -63,11 +64,14 @@ class matrix template static matrix CopyDense(size_t m, size_t n, Itr first, Itr last); // return plain storage as vector - std::vector asVector() const; + operator std::vector () const; // resize void resize(size_t m, size_t n); + // modify bounds-checks + void setPeriodic(bool periodic); + // get dimensions size_t size() const; size_t rank() const; @@ -80,11 +84,22 @@ class matrix const X& operator[](size_t i) const; // index operators: access using matrix-indices - X& operator()(size_t a, size_t b); - const X& operator()(size_t a, size_t b) const; + X& operator()(int a, int b); + const X& operator()(int a, int b) const; + + // index operators: access using matrix-indices + template::value,void>::type> + X& operator()(T a, T b); + + template::value,void>::type> + const X& operator()(T a, T b) const; // index operators: matrix-indices -> plain storage (a,b -> i) - size_t compress(size_t a, size_t b) const; + size_t compress(int a, int b) const; + + // index operators: matrix-indices -> plain storage (a,b -> i) + template::value,void>::type> + size_t compress(T a, T b) const; // index operators: plain storage -> matrix-indices (i -> a,b) std::vector decompress(size_t i) const; @@ -104,8 +119,15 @@ class matrix auto index(size_t i) const; // iterator to specific entry: access using matrix-indices - auto item(size_t a, size_t b); - auto item(size_t a, size_t b) const; + auto item(int a, int b); + auto item(int a, int b) const; + + // iterator to specific entry: access using matrix-indices + template::value,void>::type> + auto item(T a, T b); + + template::value,void>::type> + auto item(T a, T b) const; // initialization void setRandom(X lower=(X)0, X upper=(X)1); diff --git a/src/cppmat/var_regular_array.cpp b/src/cppmat/var_regular_array.cpp index 55c9699..73ef295 100644 --- a/src/cppmat/var_regular_array.cpp +++ b/src/cppmat/var_regular_array.cpp @@ -259,6 +259,17 @@ void array::chrank(size_t rank) mRank = rank; } +// ================================================================================================= +// modify bounds check +// ================================================================================================= + +template +inline +void array::setPeriodic(bool periodic) +{ + mPeriodic = periodic; +} + // ================================================================================================= // get dimensions // ================================================================================================= @@ -376,7 +387,7 @@ X& array::operator()(int a) { int na = static_cast(mShape[0]); - assert( a < na && a >= -na ); + assert( ( a < na && a >= -na ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); @@ -392,7 +403,7 @@ const X& array::operator()(int a) const { int na = static_cast(mShape[0]); - assert( a < na && a >= -na ); + assert( ( a < na && a >= -na ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); @@ -409,8 +420,8 @@ X& array::operator()(int a, int b) int na = static_cast(mShape[0]); int nb = static_cast(mShape[1]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -429,8 +440,8 @@ const X& array::operator()(int a, int b) const int na = static_cast(mShape[0]); int nb = static_cast(mShape[1]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -450,9 +461,9 @@ X& array::operator()(int a, int b, int c) int nb = static_cast(mShape[1]); int nc = static_cast(mShape[2]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -474,9 +485,9 @@ const X& array::operator()(int a, int b, int c) const int nb = static_cast(mShape[1]); int nc = static_cast(mShape[2]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -499,10 +510,10 @@ X& array::operator()(int a, int b, int c, int d) int nc = static_cast(mShape[2]); int nd = static_cast(mShape[3]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -527,10 +538,10 @@ const X& array::operator()(int a, int b, int c, int d) const int nc = static_cast(mShape[2]); int nd = static_cast(mShape[3]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -556,11 +567,11 @@ X& array::operator()(int a, int b, int c, int d, int e) int nd = static_cast(mShape[3]); int ne = static_cast(mShape[4]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -588,11 +599,11 @@ const X& array::operator()(int a, int b, int c, int d, int e) const int nd = static_cast(mShape[3]); int ne = static_cast(mShape[4]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -621,12 +632,12 @@ X& array::operator()(int a, int b, int c, int d, int e, int f) int ne = static_cast(mShape[4]); int nf = static_cast(mShape[5]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); - assert( f < nf && f >= -nf ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); + assert( ( f < nf && f >= -nf ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -657,12 +668,12 @@ const X& array::operator()(int a, int b, int c, int d, int e, int f) const int ne = static_cast(mShape[4]); int nf = static_cast(mShape[5]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); - assert( f < nf && f >= -nf ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); + assert( ( f < nf && f >= -nf ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -974,7 +985,7 @@ size_t array::compress(int a) const { int na = static_cast(mShape[0]); - assert( a < na && a >= -na ); + assert( ( a < na && a >= -na ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); @@ -990,8 +1001,8 @@ size_t array::compress(int a, int b) const int na = static_cast(mShape[0]); int nb = static_cast(mShape[1]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1010,9 +1021,9 @@ size_t array::compress(int a, int b, int c) const int nb = static_cast(mShape[1]); int nc = static_cast(mShape[2]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1034,10 +1045,10 @@ size_t array::compress(int a, int b, int c, int d) const int nc = static_cast(mShape[2]); int nd = static_cast(mShape[3]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1062,11 +1073,11 @@ size_t array::compress(int a, int b, int c, int d, int e) const int nd = static_cast(mShape[3]); int ne = static_cast(mShape[4]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1094,12 +1105,12 @@ size_t array::compress(int a, int b, int c, int d, int e, int f) const int ne = static_cast(mShape[4]); int nf = static_cast(mShape[5]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); - assert( f < nf && f >= -nf ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); + assert( ( f < nf && f >= -nf ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1333,7 +1344,7 @@ auto array::item(int a) { int na = static_cast(mShape[0]); - assert( a < na && a >= -na ); + assert( ( a < na && a >= -na ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); @@ -1349,7 +1360,7 @@ auto array::item(int a) const { int na = static_cast(mShape[0]); - assert( a < na && a >= -na ); + assert( ( a < na && a >= -na ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); @@ -1366,8 +1377,8 @@ auto array::item(int a, int b) int na = static_cast(mShape[0]); int nb = static_cast(mShape[1]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1386,8 +1397,8 @@ auto array::item(int a, int b) const int na = static_cast(mShape[0]); int nb = static_cast(mShape[1]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1407,9 +1418,9 @@ auto array::item(int a, int b, int c) int nb = static_cast(mShape[1]); int nc = static_cast(mShape[2]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1431,9 +1442,9 @@ auto array::item(int a, int b, int c) const int nb = static_cast(mShape[1]); int nc = static_cast(mShape[2]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1456,10 +1467,10 @@ auto array::item(int a, int b, int c, int d) int nc = static_cast(mShape[2]); int nd = static_cast(mShape[3]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1484,10 +1495,10 @@ auto array::item(int a, int b, int c, int d) const int nc = static_cast(mShape[2]); int nd = static_cast(mShape[3]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1513,11 +1524,11 @@ auto array::item(int a, int b, int c, int d, int e) int nd = static_cast(mShape[3]); int ne = static_cast(mShape[4]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1545,11 +1556,11 @@ auto array::item(int a, int b, int c, int d, int e) const int nd = static_cast(mShape[3]); int ne = static_cast(mShape[4]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1578,12 +1589,12 @@ auto array::item(int a, int b, int c, int d, int e, int f) int ne = static_cast(mShape[4]); int nf = static_cast(mShape[5]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); - assert( f < nf && f >= -nf ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); + assert( ( f < nf && f >= -nf ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); @@ -1614,12 +1625,12 @@ auto array::item(int a, int b, int c, int d, int e, int f) const int ne = static_cast(mShape[4]); int nf = static_cast(mShape[5]); - assert( a < na && a >= -na ); - assert( b < nb && b >= -nb ); - assert( c < nc && c >= -nc ); - assert( d < nd && d >= -nd ); - assert( e < ne && e >= -ne ); - assert( f < nf && f >= -nf ); + assert( ( a < na && a >= -na ) or mPeriodic ); + assert( ( b < nb && b >= -nb ) or mPeriodic ); + assert( ( c < nc && c >= -nc ) or mPeriodic ); + assert( ( d < nd && d >= -nd ) or mPeriodic ); + assert( ( e < ne && e >= -ne ) or mPeriodic ); + assert( ( f < nf && f >= -nf ) or mPeriodic ); size_t A = static_cast( (na+(a%na)) % na ); size_t B = static_cast( (nb+(b%nb)) % nb ); diff --git a/src/cppmat/var_regular_array.h b/src/cppmat/var_regular_array.h index 62cd876..a45b685 100644 --- a/src/cppmat/var_regular_array.h +++ b/src/cppmat/var_regular_array.h @@ -30,6 +30,7 @@ class array size_t mShape [MAX_DIM]; // number of entries along each axis size_t mStrides[MAX_DIM]; // stride length for each index std::vector mData; // data container + bool mPeriodic=false; // if true: disable bounds-check where possible public: @@ -70,6 +71,9 @@ class array void reshape(const std::vector &shape); void chrank (size_t rank); + // modify bounds-checks + void setPeriodic(bool periodic); + // get dimensions size_t size() const; size_t rank() const; @@ -166,7 +170,6 @@ class array template::value,void>::type> size_t compress(T a, T b, T c, T d, T e, T f) const; - // index operators: plain storage -> array-indices (i -> a,b,c,...) std::vector decompress(size_t i) const; diff --git a/src/cppmat/var_symmetric_matrix.cpp b/src/cppmat/var_symmetric_matrix.cpp index f66f44b..de0e157 100644 --- a/src/cppmat/var_symmetric_matrix.cpp +++ b/src/cppmat/var_symmetric_matrix.cpp @@ -229,7 +229,7 @@ matrix matrix::CopyDense(size_t m, size_t n, Iterator first, Iterator last template inline -std::vector matrix::asVector() const +matrix::operator std::vector () const { return mData; } @@ -257,6 +257,17 @@ void matrix::resize(size_t m, size_t n) if ( mSize != size ) mData.resize(mSize); } +// ================================================================================================= +// modify bounds check +// ================================================================================================= + +template +inline +void matrix::setPeriodic(bool periodic) +{ + mPeriodic = periodic; +} + // ================================================================================================= // get dimensions // ================================================================================================= @@ -347,7 +358,44 @@ const X& matrix::operator[](size_t i) const template inline -X& matrix::operator()(size_t a, size_t b) +X& matrix::operator()(int a, int b) +{ + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + assert( ( b < n && b >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + size_t B = static_cast( (n+(b%n)) % n ); + + if (A <= B) return mData[ A*N - (A-1)*A/2 + B - A ]; + else return mData[ B*N - (B-1)*B/2 + A - B ]; +} + +// ------------------------------------------------------------------------------------------------- + +template +inline +const X& matrix::operator()(int a, int b) const +{ + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + assert( ( b < n && b >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + size_t B = static_cast( (n+(b%n)) % n ); + + if (A <= B) return mData[ A*N - (A-1)*A/2 + B - A ]; + else return mData[ B*N - (B-1)*B/2 + A - B ]; +} + +// ------------------------------------------------------------------------------------------------- + +template +template +inline +X& matrix::operator()(T a, T b) { assert( a < N ); assert( b < N ); @@ -359,8 +407,9 @@ X& matrix::operator()(size_t a, size_t b) // ------------------------------------------------------------------------------------------------- template +template inline -const X& matrix::operator()(size_t a, size_t b) const +const X& matrix::operator()(T a, T b) const { assert( a < N ); assert( b < N ); @@ -375,7 +424,26 @@ const X& matrix::operator()(size_t a, size_t b) const template inline -size_t matrix::compress(size_t a, size_t b) const +size_t matrix::compress(int a, int b) const +{ + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + assert( ( b < n && b >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + size_t B = static_cast( (n+(b%n)) % n ); + + if (A <= B) return A*N - (A-1)*A/2 + B - A; + else return B*N - (B-1)*B/2 + A - B; +} + +// ------------------------------------------------------------------------------------------------- + +template +template +inline +size_t matrix::compress(T a, T b) const { assert( a < N ); assert( b < N ); @@ -500,7 +568,44 @@ auto matrix::index(size_t i) const template inline -auto matrix::item(size_t a, size_t b) +auto matrix::item(int a, int b) +{ + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + assert( ( b < n && b >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + size_t B = static_cast( (n+(b%n)) % n ); + + if (A <= B) return begin() + ( A*N - (A-1)*A/2 + B - A ); + else return begin() + ( B*N - (B-1)*B/2 + A - B ); +} + +// ------------------------------------------------------------------------------------------------- + +template +inline +auto matrix::item(int a, int b) const +{ + int n = static_cast(N); + + assert( ( a < n && a >= -n ) or mPeriodic ); + assert( ( b < n && b >= -n ) or mPeriodic ); + + size_t A = static_cast( (n+(a%n)) % n ); + size_t B = static_cast( (n+(b%n)) % n ); + + if (A <= B) return begin() + ( A*N - (A-1)*A/2 + B - A ); + else return begin() + ( B*N - (B-1)*B/2 + A - B ); +} + +// ------------------------------------------------------------------------------------------------- + +template +template +inline +auto matrix::item(T a, T b) { assert( a < N ); assert( b < N ); @@ -512,8 +617,9 @@ auto matrix::item(size_t a, size_t b) // ------------------------------------------------------------------------------------------------- template +template inline -auto matrix::item(size_t a, size_t b) const +auto matrix::item(T a, T b) const { assert( a < N ); assert( b < N ); diff --git a/src/cppmat/var_symmetric_matrix.h b/src/cppmat/var_symmetric_matrix.h index bc4b063..099770b 100644 --- a/src/cppmat/var_symmetric_matrix.h +++ b/src/cppmat/var_symmetric_matrix.h @@ -25,10 +25,11 @@ class matrix { protected: - size_t mSize=0; // total size == data.size() - static const size_t mRank=2; // rank (number of axes) - size_t N=0; // number of rows/columns - std::vector mData; // data container + size_t mSize=0; // total size == data.size() + static const size_t mRank=2; // rank (number of axes) + size_t N=0; // number of rows/columns + std::vector mData; // data container + bool mPeriodic=false; // if true: disable bounds-check where possible public: @@ -65,11 +66,14 @@ class matrix template static matrix CopyDense(size_t m, size_t n, Itr first, Itr last); // return plain storage as vector - std::vector asVector() const; + operator std::vector () const; // resize void resize(size_t m, size_t n); + // modify bounds-checks + void setPeriodic(bool periodic); + // get dimensions size_t size() const; size_t rank() const; @@ -82,11 +86,22 @@ class matrix const X& operator[](size_t i) const; // index operators: access using matrix-indices - X& operator()(size_t a, size_t b); - const X& operator()(size_t a, size_t b) const; + X& operator()(int a, int b); + const X& operator()(int a, int b) const; + + // index operators: access using matrix-indices + template::value,void>::type> + X& operator()(T a, T b); + + template::value,void>::type> + const X& operator()(T a, T b) const; // index operators: matrix-indices -> plain storage (a,b -> i) - size_t compress(size_t a, size_t b) const; + size_t compress(int a, int b) const; + + // index operators: matrix-indices -> plain storage (a,b -> i) + template::value,void>::type> + size_t compress(T a, T b) const; // index operators: plain storage -> matrix-indices (i -> a,b) std::vector decompress(size_t i) const; @@ -106,8 +121,15 @@ class matrix auto index(size_t i) const; // iterator to specific entry: access using matrix-indices - auto item(size_t a, size_t b); - auto item(size_t a, size_t b) const; + auto item(int a, int b); + auto item(int a, int b) const; + + // iterator to specific entry: access using matrix-indices + template::value,void>::type> + auto item(T a, T b); + + template::value,void>::type> + auto item(T a, T b) const; // initialization void setRandom(X lower=(X)0, X upper=(X)1);