Chun Feng Huang / MATRIX_LIB

A lite library for operations in linear algebra

MATRIX_LIB.cpp@0:33b75d52d2a7, 2017-02-10 (annotated)

Committer:
benson516
Date:
Fri Feb 10 18:31:32 2017 +0000
Revision:
0:33b75d52d2a7
Speed up the library by using constant references

Who changed what in which revision?

UserRevisionLine numberNew contents of line
benson516 0:33b75d52d2a7 1 #include "MATRIX_LIB.h"
benson516 0:33b75d52d2a7 2
benson516 0:33b75d52d2a7 3 // No pre-assignments, empty matrix
benson516 0:33b75d52d2a7 4 Mat::Mat(void):
benson516 0:33b75d52d2a7 5 _nRow(0), _nCol(0),
benson516 0:33b75d52d2a7 6 data(0, vector<float>())
benson516 0:33b75d52d2a7 7 {
benson516 0:33b75d52d2a7 8
benson516 0:33b75d52d2a7 9 }
benson516 0:33b75d52d2a7 10 // Pre-assign sizes
benson516 0:33b75d52d2a7 11 Mat::Mat(size_t nRow_in, size_t nCol_in):
benson516 0:33b75d52d2a7 12 _nRow(nRow_in), _nCol(nCol_in),
benson516 0:33b75d52d2a7 13 data(nRow_in, vector<float>(nCol_in))
benson516 0:33b75d52d2a7 14 {
benson516 0:33b75d52d2a7 15
benson516 0:33b75d52d2a7 16 }
benson516 0:33b75d52d2a7 17 // Pre-assign sizes and elements
benson516 0:33b75d52d2a7 18 Mat::Mat(size_t nRow_in, size_t nCol_in, float element_in):
benson516 0:33b75d52d2a7 19 _nRow(nRow_in), _nCol(nCol_in),
benson516 0:33b75d52d2a7 20 data(nRow_in, vector<float>(nCol_in, element_in))
benson516 0:33b75d52d2a7 21 {
benson516 0:33b75d52d2a7 22
benson516 0:33b75d52d2a7 23 }
benson516 0:33b75d52d2a7 24
benson516 0:33b75d52d2a7 25 // Public methods
benson516 0:33b75d52d2a7 26 //----------------------------------//
benson516 0:33b75d52d2a7 27 // Size ---
benson516 0:33b75d52d2a7 28 // "Get" the size of the matrix
benson516 0:33b75d52d2a7 29 size_t Mat::nRow() {
benson516 0:33b75d52d2a7 30 return this->_nRow;
benson516 0:33b75d52d2a7 31 }
benson516 0:33b75d52d2a7 32 size_t Mat::nCol() {
benson516 0:33b75d52d2a7 33 return this->_nCol;
benson516 0:33b75d52d2a7 34 }
benson516 0:33b75d52d2a7 35 size_t Mat::nElement() { // Number of elements
benson516 0:33b75d52d2a7 36 return (_nRow*_nCol);
benson516 0:33b75d52d2a7 37 }
benson516 0:33b75d52d2a7 38 // const functions
benson516 0:33b75d52d2a7 39 size_t Mat::nRow() const {
benson516 0:33b75d52d2a7 40 return this->_nRow;
benson516 0:33b75d52d2a7 41 }
benson516 0:33b75d52d2a7 42 size_t Mat::nCol() const {
benson516 0:33b75d52d2a7 43 return this->_nCol;
benson516 0:33b75d52d2a7 44 }
benson516 0:33b75d52d2a7 45 size_t Mat::nElement() const { // Number of elements
benson516 0:33b75d52d2a7 46 return (_nRow*_nCol);
benson516 0:33b75d52d2a7 47 }
benson516 0:33b75d52d2a7 48 // "Set" and "get" the new size of the matrix
benson516 0:33b75d52d2a7 49 size_t Mat::nRow(size_t nRow_in){
benson516 0:33b75d52d2a7 50 _nRow = nRow_in;
benson516 0:33b75d52d2a7 51 // Resize the matrix
benson516 0:33b75d52d2a7 52 data.resize(_nRow, vector<float>(_nCol));
benson516 0:33b75d52d2a7 53 return _nRow;
benson516 0:33b75d52d2a7 54 }
benson516 0:33b75d52d2a7 55 size_t Mat::nCol(size_t nCol_in){
benson516 0:33b75d52d2a7 56 _nCol = nCol_in;
benson516 0:33b75d52d2a7 57 // Resize the matrix
benson516 0:33b75d52d2a7 58 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 59 data[i].resize(_nCol);
benson516 0:33b75d52d2a7 60 }
benson516 0:33b75d52d2a7 61 return _nCol;
benson516 0:33b75d52d2a7 62 }
benson516 0:33b75d52d2a7 63 // "Set" and "get" the new size of the matrix with assignment to new elements
benson516 0:33b75d52d2a7 64 size_t Mat::nRow(size_t nRow_in, float element_in){
benson516 0:33b75d52d2a7 65 _nRow = nRow_in;
benson516 0:33b75d52d2a7 66 // Resize the matrix
benson516 0:33b75d52d2a7 67 data.resize(_nRow, vector<float>(_nCol, element_in));
benson516 0:33b75d52d2a7 68 return _nRow;
benson516 0:33b75d52d2a7 69 }
benson516 0:33b75d52d2a7 70 size_t Mat::nCol(size_t nCol_in, float element_in){
benson516 0:33b75d52d2a7 71 _nCol = nCol_in;
benson516 0:33b75d52d2a7 72 // Resize the matrix
benson516 0:33b75d52d2a7 73 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 74 data[i].resize(_nCol, element_in);
benson516 0:33b75d52d2a7 75 }
benson516 0:33b75d52d2a7 76 return _nCol;
benson516 0:33b75d52d2a7 77 }
benson516 0:33b75d52d2a7 78 // Resize on both row and column
benson516 0:33b75d52d2a7 79 void Mat::resize(size_t nRow_in, size_t nCol_in){
benson516 0:33b75d52d2a7 80 // Check the size
benson516 0:33b75d52d2a7 81 if (_nRow == nRow_in && _nCol == nCol_in){
benson516 0:33b75d52d2a7 82 return; // Nothing to do
benson516 0:33b75d52d2a7 83 }
benson516 0:33b75d52d2a7 84
benson516 0:33b75d52d2a7 85 // Change the nCol first to resize the old row vector
benson516 0:33b75d52d2a7 86 _nCol = nCol_in;
benson516 0:33b75d52d2a7 87 // Resize the matrix
benson516 0:33b75d52d2a7 88 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 89 data[i].resize(_nCol);
benson516 0:33b75d52d2a7 90 }
benson516 0:33b75d52d2a7 91 // Then extend the row
benson516 0:33b75d52d2a7 92 _nRow = nRow_in;
benson516 0:33b75d52d2a7 93 data.resize(_nRow, vector<float>(_nCol));
benson516 0:33b75d52d2a7 94 }
benson516 0:33b75d52d2a7 95 void Mat::resize(size_t nRow_in, size_t nCol_in, float element_in){ // Assign the new element_in to new elements
benson516 0:33b75d52d2a7 96 // Check the size
benson516 0:33b75d52d2a7 97 if (_nRow == nRow_in && _nCol == nCol_in){
benson516 0:33b75d52d2a7 98 return; // Nothing to do
benson516 0:33b75d52d2a7 99 }
benson516 0:33b75d52d2a7 100
benson516 0:33b75d52d2a7 101 // Change the nCol first to resize the old row vector
benson516 0:33b75d52d2a7 102 _nCol = nCol_in;
benson516 0:33b75d52d2a7 103 // Resize the matrix
benson516 0:33b75d52d2a7 104 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 105 data[i].resize(_nCol, element_in);
benson516 0:33b75d52d2a7 106 }
benson516 0:33b75d52d2a7 107 // Then extend the row
benson516 0:33b75d52d2a7 108 _nRow = nRow_in;
benson516 0:33b75d52d2a7 109 data.resize(_nRow, vector<float>(_nCol, element_in));
benson516 0:33b75d52d2a7 110 }
benson516 0:33b75d52d2a7 111 // Reshape the matrix, keep all the elements but change the shpae of the matrix (eg. 1 by 6 -> 2 by 3 or 3 by 2)
benson516 0:33b75d52d2a7 112 void Mat::reshape(size_t nRow_in, size_t nCol_in, bool is_Row_first){ // If is_Row_first, keep the elements' order the same in the row direction; otherwise, keep the order the same in column direction.
benson516 0:33b75d52d2a7 113 // Note: nRow_in * nCol_in should be better be the same as this->nElement(), or some elemnets will be lost or be added (add zeros, actually)
benson516 0:33b75d52d2a7 114 //
benson516 0:33b75d52d2a7 115 Mat M_new(nRow_in, nCol_in);
benson516 0:33b75d52d2a7 116 // The index for the old matrix
benson516 0:33b75d52d2a7 117 size_t io = 0;
benson516 0:33b75d52d2a7 118 size_t jo = 0;
benson516 0:33b75d52d2a7 119 //
benson516 0:33b75d52d2a7 120 if (is_Row_first){ // Row go first
benson516 0:33b75d52d2a7 121 //
benson516 0:33b75d52d2a7 122 for (size_t i = 0; i < nRow_in; ++i){
benson516 0:33b75d52d2a7 123 for (size_t j = 0; j < nCol_in; ++j){ // Row go first
benson516 0:33b75d52d2a7 124 //-------------------//
benson516 0:33b75d52d2a7 125 // In-range check:
benson516 0:33b75d52d2a7 126 if (io < this->_nRow && jo < this->_nCol){ // In range
benson516 0:33b75d52d2a7 127 M_new.data[i][j] = this->data[io][jo];
benson516 0:33b75d52d2a7 128 // Row go first
benson516 0:33b75d52d2a7 129 jo++;
benson516 0:33b75d52d2a7 130 // Change row
benson516 0:33b75d52d2a7 131 if (jo >= this->_nCol){
benson516 0:33b75d52d2a7 132 jo = 0;
benson516 0:33b75d52d2a7 133 io++;
benson516 0:33b75d52d2a7 134 }
benson516 0:33b75d52d2a7 135 //
benson516 0:33b75d52d2a7 136 }else{ // Out of range
benson516 0:33b75d52d2a7 137 M_new.data[i][j] = 0.0;
benson516 0:33b75d52d2a7 138 }
benson516 0:33b75d52d2a7 139 //-------------------//
benson516 0:33b75d52d2a7 140 }
benson516 0:33b75d52d2a7 141 }
benson516 0:33b75d52d2a7 142 }else{ // Column go first
benson516 0:33b75d52d2a7 143 //
benson516 0:33b75d52d2a7 144 for (size_t j = 0; j < nCol_in; ++j){
benson516 0:33b75d52d2a7 145 for (size_t i = 0; i < nRow_in; ++i){ // Column go first
benson516 0:33b75d52d2a7 146 //-------------------//
benson516 0:33b75d52d2a7 147 // In-range check:
benson516 0:33b75d52d2a7 148 if (io < this->_nRow && jo < this->_nCol){ // In range
benson516 0:33b75d52d2a7 149 M_new.data[i][j] = this->data[io][jo];
benson516 0:33b75d52d2a7 150 // Column go fist
benson516 0:33b75d52d2a7 151 io++;
benson516 0:33b75d52d2a7 152 // Change row
benson516 0:33b75d52d2a7 153 if (io >= this->_nRow){
benson516 0:33b75d52d2a7 154 io = 0;
benson516 0:33b75d52d2a7 155 jo++;
benson516 0:33b75d52d2a7 156 }
benson516 0:33b75d52d2a7 157 //
benson516 0:33b75d52d2a7 158 }else{ // Out of range
benson516 0:33b75d52d2a7 159 M_new.data[i][j] = 0.0;
benson516 0:33b75d52d2a7 160 }
benson516 0:33b75d52d2a7 161 //-------------------//
benson516 0:33b75d52d2a7 162 }
benson516 0:33b75d52d2a7 163 }
benson516 0:33b75d52d2a7 164 }
benson516 0:33b75d52d2a7 165 //
benson516 0:33b75d52d2a7 166 // Update this matrix
benson516 0:33b75d52d2a7 167 *this = M_new;
benson516 0:33b75d52d2a7 168 }
benson516 0:33b75d52d2a7 169 // Synchronize the _nRow and _nCol with the true size of data
benson516 0:33b75d52d2a7 170 void Mat::syncSizeInData(){
benson516 0:33b75d52d2a7 171 //
benson516 0:33b75d52d2a7 172 this->_nRow = this->data.size();
benson516 0:33b75d52d2a7 173 //
benson516 0:33b75d52d2a7 174 if (this->_nRow == 0){ // Empty
benson516 0:33b75d52d2a7 175 this->_nCol = 0;
benson516 0:33b75d52d2a7 176 }else{
benson516 0:33b75d52d2a7 177 this->_nCol = this->data[0].size();
benson516 0:33b75d52d2a7 178 }
benson516 0:33b75d52d2a7 179 }
benson516 0:33b75d52d2a7 180
benson516 0:33b75d52d2a7 181 // Assignment ---
benson516 0:33b75d52d2a7 182 void Mat::assign(std::stringstream &ss_in, size_t nRow_in, size_t nCol_in){ // The stringstream ss_in contains a nRow_in by nCol_in matrix.
benson516 0:33b75d52d2a7 183 // Resize
benson516 0:33b75d52d2a7 184 this->resize(nRow_in, nCol_in);
benson516 0:33b75d52d2a7 185 //
benson516 0:33b75d52d2a7 186 for (size_t i = 0; i < this->_nRow; ++i){
benson516 0:33b75d52d2a7 187 for (size_t j = 0; j < this->_nCol; ++j){
benson516 0:33b75d52d2a7 188 ss_in >> this->data[i][j];
benson516 0:33b75d52d2a7 189 }
benson516 0:33b75d52d2a7 190 }
benson516 0:33b75d52d2a7 191 }
benson516 0:33b75d52d2a7 192 // From primitive array in c++
benson516 0:33b75d52d2a7 193 void Mat::assign(float* Matrix_in, size_t nRow_in, size_t nCol_in){ // From primitive 1-D array in c++, Matrix_in is a nRow_in by nCol_in matrix.
benson516 0:33b75d52d2a7 194 // Resize
benson516 0:33b75d52d2a7 195 this->resize(nRow_in, nCol_in);
benson516 0:33b75d52d2a7 196 //
benson516 0:33b75d52d2a7 197 for (size_t i = 0; i < this->_nRow; ++i){
benson516 0:33b75d52d2a7 198 for (size_t j = 0; j < this->_nCol; ++j){
benson516 0:33b75d52d2a7 199 // this->data[i][j] = Matrix_in[i][j];
benson516 0:33b75d52d2a7 200 this->data[i][j] = *Matrix_in;
benson516 0:33b75d52d2a7 201 Matrix_in++;
benson516 0:33b75d52d2a7 202 }
benson516 0:33b75d52d2a7 203 }
benson516 0:33b75d52d2a7 204 }
benson516 0:33b75d52d2a7 205 void Mat::assign(const vector<float> &vec_in, bool is_RowVec){ // From 1-D vector
benson516 0:33b75d52d2a7 206 size_t n = vec_in.size();
benson516 0:33b75d52d2a7 207 //
benson516 0:33b75d52d2a7 208 if (is_RowVec){ // Row vector
benson516 0:33b75d52d2a7 209 //
benson516 0:33b75d52d2a7 210 this->resize(1,n);
benson516 0:33b75d52d2a7 211 //
benson516 0:33b75d52d2a7 212 for (size_t j = 0; j < _nCol; ++j){
benson516 0:33b75d52d2a7 213 this->data[0][j] = vec_in[j];
benson516 0:33b75d52d2a7 214 }
benson516 0:33b75d52d2a7 215
benson516 0:33b75d52d2a7 216 //
benson516 0:33b75d52d2a7 217 }else{ // Column vector
benson516 0:33b75d52d2a7 218 //
benson516 0:33b75d52d2a7 219 this->resize(n,1);
benson516 0:33b75d52d2a7 220 //
benson516 0:33b75d52d2a7 221 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 222 this->data[i][0] = vec_in[i];
benson516 0:33b75d52d2a7 223 }
benson516 0:33b75d52d2a7 224 //
benson516 0:33b75d52d2a7 225 }
benson516 0:33b75d52d2a7 226 }
benson516 0:33b75d52d2a7 227 void Mat::assign(const vector<vector<float> > &MatData_in){ // A = MatData_in, assign a primitive matrix directly
benson516 0:33b75d52d2a7 228 this->data = MatData_in;
benson516 0:33b75d52d2a7 229 this->syncSizeInData();
benson516 0:33b75d52d2a7 230 }
benson516 0:33b75d52d2a7 231 // Partial asignment
benson516 0:33b75d52d2a7 232 void Mat::setPart(const Mat &M_part, size_t m_from, size_t n_from){ // The block starts from (m_from, n_from)
benson516 0:33b75d52d2a7 233 size_t m_to, n_to;
benson516 0:33b75d52d2a7 234 // The index of the last element of the inserted block
benson516 0:33b75d52d2a7 235 m_to = m_from + M_part.nRow() - 1;
benson516 0:33b75d52d2a7 236 n_to = n_from + M_part.nCol() - 1;
benson516 0:33b75d52d2a7 237
benson516 0:33b75d52d2a7 238 //-------------------------------//
benson516 0:33b75d52d2a7 239 // Note: if the region of the insertion block exceed the region of this matrix,
benson516 0:33b75d52d2a7 240 // resize to a biger one and insert zero at the blank sides.
benson516 0:33b75d52d2a7 241
benson516 0:33b75d52d2a7 242 // The size of the new matrix
benson516 0:33b75d52d2a7 243 size_t M,N;
benson516 0:33b75d52d2a7 244 //
benson516 0:33b75d52d2a7 245 M = this->_nRow;
benson516 0:33b75d52d2a7 246 N = this->_nCol;
benson516 0:33b75d52d2a7 247 //
benson516 0:33b75d52d2a7 248 // Check if block exceeds the region
benson516 0:33b75d52d2a7 249 if (m_to >= M ){
benson516 0:33b75d52d2a7 250 M = m_to + 1;
benson516 0:33b75d52d2a7 251 }
benson516 0:33b75d52d2a7 252 if(n_to >= N){
benson516 0:33b75d52d2a7 253 N = n_to + 1;
benson516 0:33b75d52d2a7 254 }
benson516 0:33b75d52d2a7 255 // Deside if the resize action is needed
benson516 0:33b75d52d2a7 256 this->resize(M, N, 0.0);
benson516 0:33b75d52d2a7 257
benson516 0:33b75d52d2a7 258 // Start the insertion
benson516 0:33b75d52d2a7 259 //
benson516 0:33b75d52d2a7 260 // Initialization
benson516 0:33b75d52d2a7 261 size_t ip = 0;
benson516 0:33b75d52d2a7 262 size_t jp = 0;
benson516 0:33b75d52d2a7 263 for (size_t i = m_from; i <= m_to; ++i ){
benson516 0:33b75d52d2a7 264 jp = 0;
benson516 0:33b75d52d2a7 265 for (size_t j = n_from; j <= n_to; ++j){
benson516 0:33b75d52d2a7 266 // M_part.data[ip][jp] = data[i][j]; // <- get
benson516 0:33b75d52d2a7 267 this->data[i][j] = M_part.data[ip][jp]; // <- insert
benson516 0:33b75d52d2a7 268 //
benson516 0:33b75d52d2a7 269 jp++;
benson516 0:33b75d52d2a7 270 }
benson516 0:33b75d52d2a7 271 ip++;
benson516 0:33b75d52d2a7 272 }
benson516 0:33b75d52d2a7 273
benson516 0:33b75d52d2a7 274 }
benson516 0:33b75d52d2a7 275 // Spetial assignments
benson516 0:33b75d52d2a7 276 void Mat::zeros(size_t nRow_in, size_t nCol_in){ // zeros(m,n)
benson516 0:33b75d52d2a7 277 resize(0,0);
benson516 0:33b75d52d2a7 278 resize(nRow_in,nCol_in, 0.0);
benson516 0:33b75d52d2a7 279 }
benson516 0:33b75d52d2a7 280 void Mat::ones(size_t nRow_in, size_t nCol_in){ // ones(m,n)
benson516 0:33b75d52d2a7 281 resize(0,0);
benson516 0:33b75d52d2a7 282 resize(nRow_in,nCol_in, 1.0);
benson516 0:33b75d52d2a7 283 }
benson516 0:33b75d52d2a7 284 void Mat::eye(size_t nRow_in, size_t nCol_in){ // Identity matrix, eye(m,n)
benson516 0:33b75d52d2a7 285 zeros(nRow_in, nCol_in);
benson516 0:33b75d52d2a7 286 //
benson516 0:33b75d52d2a7 287 size_t K;
benson516 0:33b75d52d2a7 288 //
benson516 0:33b75d52d2a7 289 // min()
benson516 0:33b75d52d2a7 290 if (nRow_in < nCol_in){
benson516 0:33b75d52d2a7 291 K = nRow_in;
benson516 0:33b75d52d2a7 292 }else{
benson516 0:33b75d52d2a7 293 K = nCol_in;
benson516 0:33b75d52d2a7 294 }
benson516 0:33b75d52d2a7 295 //
benson516 0:33b75d52d2a7 296 for (size_t i = 0; i < K; ++i){
benson516 0:33b75d52d2a7 297 data[i][i] = 1.0;
benson516 0:33b75d52d2a7 298 }
benson516 0:33b75d52d2a7 299 //
benson516 0:33b75d52d2a7 300 }
benson516 0:33b75d52d2a7 301 void Mat::diag(const Mat &M_diag){ // Transform the row/column vector to a diagonal matrix and assign into this matrix
benson516 0:33b75d52d2a7 302 size_t M = M_diag.nRow();
benson516 0:33b75d52d2a7 303 size_t N = M_diag.nCol();
benson516 0:33b75d52d2a7 304 size_t nE = M_diag.nElement();
benson516 0:33b75d52d2a7 305 //
benson516 0:33b75d52d2a7 306 this->zeros(nE,nE); // nE by nE matrix
benson516 0:33b75d52d2a7 307 //
benson516 0:33b75d52d2a7 308 size_t id = 0;
benson516 0:33b75d52d2a7 309 size_t jd = 0;
benson516 0:33b75d52d2a7 310 for (size_t i = 0; i < nE; ++i){
benson516 0:33b75d52d2a7 311 this->data[i][i] = M_diag.data[id][jd];
benson516 0:33b75d52d2a7 312 jd++;
benson516 0:33b75d52d2a7 313 if (jd >= N){
benson516 0:33b75d52d2a7 314 jd = 0;
benson516 0:33b75d52d2a7 315 id++;
benson516 0:33b75d52d2a7 316 }
benson516 0:33b75d52d2a7 317 }
benson516 0:33b75d52d2a7 318 }
benson516 0:33b75d52d2a7 319
benson516 0:33b75d52d2a7 320
benson516 0:33b75d52d2a7 321 // Get elements ---
benson516 0:33b75d52d2a7 322 float Mat::at(size_t m, size_t n){ // Get the element at (m,n)
benson516 0:33b75d52d2a7 323 return data[m][n];
benson516 0:33b75d52d2a7 324 }
benson516 0:33b75d52d2a7 325 Mat Mat::getPart(size_t m_from, size_t m_to, size_t n_from, size_t n_to){ // Get partial, M_part = M(m_from:m_to, n_from:n_to)
benson516 0:33b75d52d2a7 326 int M = m_to - m_from + 1;
benson516 0:33b75d52d2a7 327 int N = n_to - n_from + 1;
benson516 0:33b75d52d2a7 328 //
benson516 0:33b75d52d2a7 329 if (M < 0 || N < 0){
benson516 0:33b75d52d2a7 330 Mat Empty;
benson516 0:33b75d52d2a7 331 return Empty;
benson516 0:33b75d52d2a7 332 }
benson516 0:33b75d52d2a7 333 //
benson516 0:33b75d52d2a7 334 Mat M_part(M,N);
benson516 0:33b75d52d2a7 335
benson516 0:33b75d52d2a7 336 //
benson516 0:33b75d52d2a7 337 // Initialization
benson516 0:33b75d52d2a7 338 size_t ip = 0;
benson516 0:33b75d52d2a7 339 size_t jp = 0;
benson516 0:33b75d52d2a7 340 //
benson516 0:33b75d52d2a7 341 for (size_t i = m_from; i <= m_to; ++i ){
benson516 0:33b75d52d2a7 342 jp = 0;
benson516 0:33b75d52d2a7 343 for (size_t j = n_from; j <= n_to; ++j){
benson516 0:33b75d52d2a7 344 // In-range check
benson516 0:33b75d52d2a7 345 if (i < this->_nRow && j < this->_nCol){
benson516 0:33b75d52d2a7 346 M_part.data[ip][jp] = this->data[i][j];
benson516 0:33b75d52d2a7 347 }else{ // Out of bound, set elements to zeros
benson516 0:33b75d52d2a7 348 M_part.data[ip][jp] = 0.0;
benson516 0:33b75d52d2a7 349 }
benson516 0:33b75d52d2a7 350 //
benson516 0:33b75d52d2a7 351 jp++;
benson516 0:33b75d52d2a7 352 }
benson516 0:33b75d52d2a7 353 ip++;
benson516 0:33b75d52d2a7 354 }
benson516 0:33b75d52d2a7 355 return M_part;
benson516 0:33b75d52d2a7 356 }
benson516 0:33b75d52d2a7 357 Mat Mat::getCol(size_t n){ // Get a specific column vector in 2-D matrix form
benson516 0:33b75d52d2a7 358 return getPart(0, (this->_nRow-1), n, n);
benson516 0:33b75d52d2a7 359 }
benson516 0:33b75d52d2a7 360 Mat Mat::getRow(size_t m){ // Get a specific row vector in 2-D matrix form
benson516 0:33b75d52d2a7 361 return getPart(m, m, 0, (this->_nCol-1));
benson516 0:33b75d52d2a7 362 }
benson516 0:33b75d52d2a7 363 vector<float> Mat::getColVec(size_t n){ // Return a c++ vector, M(:,n)
benson516 0:33b75d52d2a7 364 vector<float> ColVec(this->_nRow);
benson516 0:33b75d52d2a7 365 for (size_t i = 0; i < this->_nRow; ++i){
benson516 0:33b75d52d2a7 366 ColVec[i] = this->data[i][n];
benson516 0:33b75d52d2a7 367 }
benson516 0:33b75d52d2a7 368 return ColVec;
benson516 0:33b75d52d2a7 369 }
benson516 0:33b75d52d2a7 370 vector<float> Mat::getRowVec(size_t m){ // Return a c++ vector, M(m,:)
benson516 0:33b75d52d2a7 371 return this->data[m];
benson516 0:33b75d52d2a7 372 }
benson516 0:33b75d52d2a7 373
benson516 0:33b75d52d2a7 374 // Print out matrix ---
benson516 0:33b75d52d2a7 375 std::string Mat::print(void){ // Print this matrix out as string
benson516 0:33b75d52d2a7 376 /*
benson516 0:33b75d52d2a7 377 std::string str_out;
benson516 0:33b75d52d2a7 378 //
benson516 0:33b75d52d2a7 379 std::string endl("\n"); // The end line character
benson516 0:33b75d52d2a7 380
benson516 0:33b75d52d2a7 381 //
benson516 0:33b75d52d2a7 382 str_out += endl;
benson516 0:33b75d52d2a7 383 for (size_t i = 0; i < this->_nRow; ++i){
benson516 0:33b75d52d2a7 384 for (size_t j = 0; j < this->_nCol; ++j){
benson516 0:33b75d52d2a7 385 str_out += std::to_string( this->data[i][j] ) + "\t";
benson516 0:33b75d52d2a7 386 }
benson516 0:33b75d52d2a7 387 str_out += endl;
benson516 0:33b75d52d2a7 388 }
benson516 0:33b75d52d2a7 389 str_out += endl;
benson516 0:33b75d52d2a7 390 //
benson516 0:33b75d52d2a7 391 return str_out;
benson516 0:33b75d52d2a7 392 */
benson516 0:33b75d52d2a7 393
benson516 0:33b75d52d2a7 394 // Using stringstream
benson516 0:33b75d52d2a7 395 std::stringstream ss_out;
benson516 0:33b75d52d2a7 396 //
benson516 0:33b75d52d2a7 397 ss_out << std::endl;
benson516 0:33b75d52d2a7 398 //
benson516 0:33b75d52d2a7 399 for (size_t i = 0; i < this->_nRow; ++i){
benson516 0:33b75d52d2a7 400 for (size_t j = 0; j < this->_nCol; ++j){
benson516 0:33b75d52d2a7 401 ss_out << this->data[i][j] << "\t";
benson516 0:33b75d52d2a7 402 }
benson516 0:33b75d52d2a7 403 ss_out << std::endl;
benson516 0:33b75d52d2a7 404 }
benson516 0:33b75d52d2a7 405 ss_out << std::endl;
benson516 0:33b75d52d2a7 406 //
benson516 0:33b75d52d2a7 407 return ss_out.str();
benson516 0:33b75d52d2a7 408 }
benson516 0:33b75d52d2a7 409
benson516 0:33b75d52d2a7 410
benson516 0:33b75d52d2a7 411 // Operations =====
benson516 0:33b75d52d2a7 412 // Comparison
benson516 0:33b75d52d2a7 413 bool Mat::is_equal(const Mat &M_in){
benson516 0:33b75d52d2a7 414 //
benson516 0:33b75d52d2a7 415 if (M_in.nRow() != this->_nRow || M_in.nCol() != this->_nCol){
benson516 0:33b75d52d2a7 416 return false;
benson516 0:33b75d52d2a7 417 }
benson516 0:33b75d52d2a7 418
benson516 0:33b75d52d2a7 419 // They are of the same size.
benson516 0:33b75d52d2a7 420 for (size_t i = 0; i < this->_nRow; ++i){
benson516 0:33b75d52d2a7 421 for (size_t j = 0; j < this->_nCol; ++j){
benson516 0:33b75d52d2a7 422 if ( this->data[i][j] != M_in.data[i][j] ){
benson516 0:33b75d52d2a7 423 return false;
benson516 0:33b75d52d2a7 424 }
benson516 0:33b75d52d2a7 425 }
benson516 0:33b75d52d2a7 426 }
benson516 0:33b75d52d2a7 427 return true;
benson516 0:33b75d52d2a7 428 }
benson516 0:33b75d52d2a7 429
benson516 0:33b75d52d2a7 430 // Self-operation (affecting on this matrix) ---
benson516 0:33b75d52d2a7 431 void Mat::scaleUp(float scale){ // M *= scale
benson516 0:33b75d52d2a7 432 //
benson516 0:33b75d52d2a7 433 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 434 MP::Get_VectorScaleUp(data[i], scale);
benson516 0:33b75d52d2a7 435 }
benson516 0:33b75d52d2a7 436 }
benson516 0:33b75d52d2a7 437 void Mat::scaleUp_Mat(const Mat &M_right){ // M = M.times(M_right), element-wise multiplication
benson516 0:33b75d52d2a7 438 //
benson516 0:33b75d52d2a7 439 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 440 for (size_t j = 0; j < _nCol; ++j){
benson516 0:33b75d52d2a7 441 data[i][j] *= M_right.data[i][j];
benson516 0:33b75d52d2a7 442 }
benson516 0:33b75d52d2a7 443 }
benson516 0:33b75d52d2a7 444 }
benson516 0:33b75d52d2a7 445 void Mat::increase(float scale){ // M += scale, for each element in M
benson516 0:33b75d52d2a7 446 //
benson516 0:33b75d52d2a7 447 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 448 for (size_t j = 0; j < _nCol; ++j){
benson516 0:33b75d52d2a7 449 data[i][j] += scale;
benson516 0:33b75d52d2a7 450 }
benson516 0:33b75d52d2a7 451 }
benson516 0:33b75d52d2a7 452 }
benson516 0:33b75d52d2a7 453 void Mat::decrease(float scale){ // M -= scale, for each element in M
benson516 0:33b75d52d2a7 454 //
benson516 0:33b75d52d2a7 455 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 456 for (size_t j = 0; j < _nCol; ++j){
benson516 0:33b75d52d2a7 457 data[i][j] -= scale;
benson516 0:33b75d52d2a7 458 }
benson516 0:33b75d52d2a7 459 }
benson516 0:33b75d52d2a7 460 }
benson516 0:33b75d52d2a7 461 void Mat::increase(const Mat &M_right){ // M += M_right
benson516 0:33b75d52d2a7 462 //
benson516 0:33b75d52d2a7 463 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 464 MP::Get_VectorIncrement(data[i], M_right.data[i], false); // +=
benson516 0:33b75d52d2a7 465 }
benson516 0:33b75d52d2a7 466 }
benson516 0:33b75d52d2a7 467 void Mat::decrease(const Mat &M_right){ // M -= M_right
benson516 0:33b75d52d2a7 468 //
benson516 0:33b75d52d2a7 469 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 470 MP::Get_VectorIncrement(data[i], M_right.data[i], true); // -=
benson516 0:33b75d52d2a7 471 }
benson516 0:33b75d52d2a7 472 }
benson516 0:33b75d52d2a7 473 // Single-operation, output another matrix ---
benson516 0:33b75d52d2a7 474 // Transpose
benson516 0:33b75d52d2a7 475 Mat& Mat::T(void){ // Return the transposed version of this matrix
benson516 0:33b75d52d2a7 476 //
benson516 0:33b75d52d2a7 477 // Mat M_t(_nCol, _nRow);
benson516 0:33b75d52d2a7 478 static Mat M_t;
benson516 0:33b75d52d2a7 479 M_t.resize(_nCol, _nRow);
benson516 0:33b75d52d2a7 480 //
benson516 0:33b75d52d2a7 481 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 482 for (size_t j = 0; j < _nCol; ++j){
benson516 0:33b75d52d2a7 483 M_t.data[j][i] = data[i][j];
benson516 0:33b75d52d2a7 484 }
benson516 0:33b75d52d2a7 485 }
benson516 0:33b75d52d2a7 486 return M_t;
benson516 0:33b75d52d2a7 487 }
benson516 0:33b75d52d2a7 488 // Inverse
benson516 0:33b75d52d2a7 489 Mat& Mat::inverse(void){ // Return the inversion of this matrix
benson516 0:33b75d52d2a7 490 // The output matrix
benson516 0:33b75d52d2a7 491 static Mat M_inv;
benson516 0:33b75d52d2a7 492
benson516 0:33b75d52d2a7 493 if (this->_nRow != this->_nCol){
benson516 0:33b75d52d2a7 494 M_inv.resize(0,0); // Empty matrix
benson516 0:33b75d52d2a7 495 return M_inv;
benson516 0:33b75d52d2a7 496 }
benson516 0:33b75d52d2a7 497
benson516 0:33b75d52d2a7 498
benson516 0:33b75d52d2a7 499 // M_inv.resize( _nCol, _nRow); // Maybe we will impliment the psodo inverse next time
benson516 0:33b75d52d2a7 500
benson516 0:33b75d52d2a7 501 //
benson516 0:33b75d52d2a7 502 M_inv.data = MP::MatrixInversion(this->data);
benson516 0:33b75d52d2a7 503
benson516 0:33b75d52d2a7 504 // Update the _nRow and _nCol
benson516 0:33b75d52d2a7 505 M_inv.syncSizeInData();
benson516 0:33b75d52d2a7 506
benson516 0:33b75d52d2a7 507 return M_inv;
benson516 0:33b75d52d2a7 508 }
benson516 0:33b75d52d2a7 509 Mat& Mat::intPower(int power){ // M^power, M^0 -> I, M^-1 -> M_inverse
benson516 0:33b75d52d2a7 510 // The output matrix
benson516 0:33b75d52d2a7 511 static Mat M_out;
benson516 0:33b75d52d2a7 512
benson516 0:33b75d52d2a7 513 // Spetial case M^1 -> M
benson516 0:33b75d52d2a7 514 if (power == 1){
benson516 0:33b75d52d2a7 515 M_out = *this;
benson516 0:33b75d52d2a7 516 return M_out;
benson516 0:33b75d52d2a7 517 }
benson516 0:33b75d52d2a7 518 // Spetial case M^0 -> I
benson516 0:33b75d52d2a7 519 if (power == 0){
benson516 0:33b75d52d2a7 520 M_out.eye(_nRow,_nCol);
benson516 0:33b75d52d2a7 521 return M_out;
benson516 0:33b75d52d2a7 522 }
benson516 0:33b75d52d2a7 523
benson516 0:33b75d52d2a7 524 // Check for square matrix
benson516 0:33b75d52d2a7 525 if (this->_nRow != this->_nCol){
benson516 0:33b75d52d2a7 526 M_out.resize(0,0); // Error, return empty matrix
benson516 0:33b75d52d2a7 527 return M_out;
benson516 0:33b75d52d2a7 528 }
benson516 0:33b75d52d2a7 529 //-----------------------------------------------//
benson516 0:33b75d52d2a7 530
benson516 0:33b75d52d2a7 531 // Spetial case M^-1 -> this->inverse()
benson516 0:33b75d52d2a7 532 if (power == -1){
benson516 0:33b75d52d2a7 533 // M^-1, special case for speed up the frequent useage
benson516 0:33b75d52d2a7 534 return this->inverse();
benson516 0:33b75d52d2a7 535 }
benson516 0:33b75d52d2a7 536 //
benson516 0:33b75d52d2a7 537
benson516 0:33b75d52d2a7 538 // Other cases
benson516 0:33b75d52d2a7 539 //
benson516 0:33b75d52d2a7 540 if (power > 1){ // power > 1
benson516 0:33b75d52d2a7 541 M_out = *this;
benson516 0:33b75d52d2a7 542 for (size_t i = 0; i < (power-1); ++i){
benson516 0:33b75d52d2a7 543 M_out = this->dot(M_out);
benson516 0:33b75d52d2a7 544 }
benson516 0:33b75d52d2a7 545 }else if (power < -1){ // power < -1
benson516 0:33b75d52d2a7 546 Mat M_inv = this->inverse();
benson516 0:33b75d52d2a7 547 M_out = M_inv; // this->inverse()
benson516 0:33b75d52d2a7 548 for (size_t i = 0; i < (-power - 1); ++i){
benson516 0:33b75d52d2a7 549 M_out = M_inv.dot(M_out);
benson516 0:33b75d52d2a7 550 }
benson516 0:33b75d52d2a7 551 }else{
benson516 0:33b75d52d2a7 552 // Nothing, power \in {-1, 0, 1} has been processed
benson516 0:33b75d52d2a7 553 }
benson516 0:33b75d52d2a7 554 //
benson516 0:33b75d52d2a7 555 return M_out;
benson516 0:33b75d52d2a7 556
benson516 0:33b75d52d2a7 557 }
benson516 0:33b75d52d2a7 558
benson516 0:33b75d52d2a7 559 // Duo-operation, output another matrix ---
benson516 0:33b75d52d2a7 560 // Plus
benson516 0:33b75d52d2a7 561 Mat& Mat::plus(float scale){ // (M + scale), for each element in M
benson516 0:33b75d52d2a7 562 // Create the output matrix
benson516 0:33b75d52d2a7 563 // Mat M_out(_nRow,_nCol);
benson516 0:33b75d52d2a7 564 static Mat M_out;
benson516 0:33b75d52d2a7 565 M_out.resize(_nRow,_nCol);
benson516 0:33b75d52d2a7 566 //
benson516 0:33b75d52d2a7 567 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 568 for (size_t j = 0; j < _nCol; ++j){
benson516 0:33b75d52d2a7 569 M_out.data[i][j] = this->data[i][j] + scale;
benson516 0:33b75d52d2a7 570 }
benson516 0:33b75d52d2a7 571 }
benson516 0:33b75d52d2a7 572 return M_out;
benson516 0:33b75d52d2a7 573 }
benson516 0:33b75d52d2a7 574 Mat& Mat::minus(float scale){ // (M - scale), for each element in M
benson516 0:33b75d52d2a7 575 // Create the output matrix
benson516 0:33b75d52d2a7 576 // Mat M_out(_nRow,_nCol);
benson516 0:33b75d52d2a7 577 static Mat M_out;
benson516 0:33b75d52d2a7 578 M_out.resize(_nRow,_nCol);
benson516 0:33b75d52d2a7 579 //
benson516 0:33b75d52d2a7 580 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 581 for (size_t j = 0; j < _nCol; ++j){
benson516 0:33b75d52d2a7 582 M_out.data[i][j] = this->data[i][j] - scale;
benson516 0:33b75d52d2a7 583 }
benson516 0:33b75d52d2a7 584 }
benson516 0:33b75d52d2a7 585 return M_out;
benson516 0:33b75d52d2a7 586 }
benson516 0:33b75d52d2a7 587 Mat& Mat::minus(float scale, bool is_reversed){ // is_reversed -> (scale - M), for each element in M
benson516 0:33b75d52d2a7 588 // Create the output matrix
benson516 0:33b75d52d2a7 589 // Mat M_out(_nRow,_nCol);
benson516 0:33b75d52d2a7 590 static Mat M_out;
benson516 0:33b75d52d2a7 591 M_out.resize(_nRow,_nCol);
benson516 0:33b75d52d2a7 592 //
benson516 0:33b75d52d2a7 593 if (is_reversed){ // Reversed
benson516 0:33b75d52d2a7 594 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 595 for (size_t j = 0; j < _nCol; ++j){
benson516 0:33b75d52d2a7 596 M_out.data[i][j] = scale - this->data[i][j]; // (scale - M)
benson516 0:33b75d52d2a7 597 }
benson516 0:33b75d52d2a7 598 }
benson516 0:33b75d52d2a7 599 }else{
benson516 0:33b75d52d2a7 600 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 601 for (size_t j = 0; j < _nCol; ++j){
benson516 0:33b75d52d2a7 602 M_out.data[i][j] = this->data[i][j] - scale; // (M - scale)
benson516 0:33b75d52d2a7 603 }
benson516 0:33b75d52d2a7 604 }
benson516 0:33b75d52d2a7 605 }
benson516 0:33b75d52d2a7 606 return M_out;
benson516 0:33b75d52d2a7 607 }
benson516 0:33b75d52d2a7 608 Mat& Mat::plus(const Mat &M_right){
benson516 0:33b75d52d2a7 609 // Create the output matrix
benson516 0:33b75d52d2a7 610 // Mat M_out(_nRow,_nCol);
benson516 0:33b75d52d2a7 611 static Mat M_out;
benson516 0:33b75d52d2a7 612 M_out.resize(_nRow,_nCol);
benson516 0:33b75d52d2a7 613 //
benson516 0:33b75d52d2a7 614 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 615 for (size_t j = 0; j < _nCol; ++j){
benson516 0:33b75d52d2a7 616 M_out.data[i][j] = this->data[i][j] + M_right.data[i][j];
benson516 0:33b75d52d2a7 617 }
benson516 0:33b75d52d2a7 618 }
benson516 0:33b75d52d2a7 619 return M_out;
benson516 0:33b75d52d2a7 620 }
benson516 0:33b75d52d2a7 621 Mat& Mat::minus(const Mat &M_right){
benson516 0:33b75d52d2a7 622 // Create the output matrix
benson516 0:33b75d52d2a7 623 // Mat M_out(_nRow,_nCol);
benson516 0:33b75d52d2a7 624 static Mat M_out;
benson516 0:33b75d52d2a7 625 M_out.resize(_nRow,_nCol);
benson516 0:33b75d52d2a7 626 //
benson516 0:33b75d52d2a7 627 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 628 for (size_t j = 0; j < _nCol; ++j){
benson516 0:33b75d52d2a7 629 M_out.data[i][j] = this->data[i][j] - M_right.data[i][j];
benson516 0:33b75d52d2a7 630 }
benson516 0:33b75d52d2a7 631 }
benson516 0:33b75d52d2a7 632 return M_out;
benson516 0:33b75d52d2a7 633 }
benson516 0:33b75d52d2a7 634 // Concatenation
benson516 0:33b75d52d2a7 635 Mat Mat::cat_below(const Mat &M_in){ // Below this matrix, [A; B]
benson516 0:33b75d52d2a7 636 Mat M_cat;
benson516 0:33b75d52d2a7 637 // Using setPart(), the size will be automatically changed.
benson516 0:33b75d52d2a7 638 //------------------------------//
benson516 0:33b75d52d2a7 639 // Put this matrix at (0, 0)
benson516 0:33b75d52d2a7 640 M_cat.setPart(*this, 0, 0);
benson516 0:33b75d52d2a7 641 // Put M_in at (_nRow, 0)
benson516 0:33b75d52d2a7 642 M_cat.setPart(M_in, this->_nRow, 0);
benson516 0:33b75d52d2a7 643 //------------------------------//
benson516 0:33b75d52d2a7 644 return M_cat;
benson516 0:33b75d52d2a7 645 }
benson516 0:33b75d52d2a7 646 Mat Mat::cat_right(const Mat &M_in){ // Right-side of this matrix, [A, B]
benson516 0:33b75d52d2a7 647 Mat M_cat;
benson516 0:33b75d52d2a7 648 // Using setPart(), the size will be automatically changed.
benson516 0:33b75d52d2a7 649 //------------------------------//
benson516 0:33b75d52d2a7 650 // Put this matrix at (0, 0)
benson516 0:33b75d52d2a7 651 M_cat.setPart(*this, 0, 0);
benson516 0:33b75d52d2a7 652 // Put M_in at (0, _nCol)
benson516 0:33b75d52d2a7 653 M_cat.setPart(M_in, 0, this->_nCol);
benson516 0:33b75d52d2a7 654 //------------------------------//
benson516 0:33b75d52d2a7 655 return M_cat;
benson516 0:33b75d52d2a7 656 }
benson516 0:33b75d52d2a7 657 Mat Mat::cat(const Mat &M_in, bool is_horizontal){ // is_horizontal --> cat_Right(); otherwise --> cat_Below()
benson516 0:33b75d52d2a7 658 //
benson516 0:33b75d52d2a7 659 if(is_horizontal){ // [A, B]
benson516 0:33b75d52d2a7 660 return this->cat_right(M_in);
benson516 0:33b75d52d2a7 661 }else{ // [A; B]
benson516 0:33b75d52d2a7 662 return this->cat_below(M_in);
benson516 0:33b75d52d2a7 663 }
benson516 0:33b75d52d2a7 664 }
benson516 0:33b75d52d2a7 665
benson516 0:33b75d52d2a7 666 // Scalar/Element-wise multiplication
benson516 0:33b75d52d2a7 667 Mat& Mat::times(float scale){ // Scalar multiplication
benson516 0:33b75d52d2a7 668 // Create the output matrix
benson516 0:33b75d52d2a7 669 // Mat M_out(_nRow,_nCol);
benson516 0:33b75d52d2a7 670 //
benson516 0:33b75d52d2a7 671 static Mat M_out;
benson516 0:33b75d52d2a7 672 M_out.resize(_nRow, _nCol);
benson516 0:33b75d52d2a7 673
benson516 0:33b75d52d2a7 674
benson516 0:33b75d52d2a7 675 if(scale == -1.0){
benson516 0:33b75d52d2a7 676 //
benson516 0:33b75d52d2a7 677 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 678 for (size_t j = 0; j < _nCol; ++j){
benson516 0:33b75d52d2a7 679 M_out.data[i][j] = -(this->data[i][j]);
benson516 0:33b75d52d2a7 680 }
benson516 0:33b75d52d2a7 681 }
benson516 0:33b75d52d2a7 682 }else{
benson516 0:33b75d52d2a7 683 //
benson516 0:33b75d52d2a7 684 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 685 for (size_t j = 0; j < _nCol; ++j){
benson516 0:33b75d52d2a7 686 M_out.data[i][j] = this->data[i][j] * scale;
benson516 0:33b75d52d2a7 687 }
benson516 0:33b75d52d2a7 688 }
benson516 0:33b75d52d2a7 689 }
benson516 0:33b75d52d2a7 690 return M_out;
benson516 0:33b75d52d2a7 691 }
benson516 0:33b75d52d2a7 692 Mat& Mat::times(const Mat &M_right){ // Element-wise multiplication
benson516 0:33b75d52d2a7 693 // Create the output matrix
benson516 0:33b75d52d2a7 694 // Mat M_out(_nRow,_nCol);
benson516 0:33b75d52d2a7 695 static Mat M_out;
benson516 0:33b75d52d2a7 696 M_out.resize(_nRow,_nCol);
benson516 0:33b75d52d2a7 697 //
benson516 0:33b75d52d2a7 698 for (size_t i = 0; i < _nRow; ++i){
benson516 0:33b75d52d2a7 699 for (size_t j = 0; j < _nCol; ++j){
benson516 0:33b75d52d2a7 700 M_out.data[i][j] = this->data[i][j] * M_right.data[i][j];
benson516 0:33b75d52d2a7 701 }
benson516 0:33b75d52d2a7 702 }
benson516 0:33b75d52d2a7 703 return M_out;
benson516 0:33b75d52d2a7 704 }
benson516 0:33b75d52d2a7 705
benson516 0:33b75d52d2a7 706
benson516 0:33b75d52d2a7 707 // Matrix multiplication
benson516 0:33b75d52d2a7 708 // Note: this matrix is the "left" one
benson516 0:33b75d52d2a7 709 Mat& Mat::dot(const Mat &M_right){ // Similar to the nomenclature of numpy in Python
benson516 0:33b75d52d2a7 710
benson516 0:33b75d52d2a7 711 // Size check
benson516 0:33b75d52d2a7 712 // mxn = (mxk)*(kxn)
benson516 0:33b75d52d2a7 713 size_t M, K, N;
benson516 0:33b75d52d2a7 714 //
benson516 0:33b75d52d2a7 715 M = _nRow;
benson516 0:33b75d52d2a7 716 K = _nCol;
benson516 0:33b75d52d2a7 717 N = M_right.nCol();
benson516 0:33b75d52d2a7 718
benson516 0:33b75d52d2a7 719 // The output matrix
benson516 0:33b75d52d2a7 720 // Mat M_out(M, N);
benson516 0:33b75d52d2a7 721 static Mat M_out;
benson516 0:33b75d52d2a7 722 M_out.resize(M, N);
benson516 0:33b75d52d2a7 723
benson516 0:33b75d52d2a7 724
benson516 0:33b75d52d2a7 725 // cout << "here in\n";
benson516 0:33b75d52d2a7 726 float* ptr_data = NULL;
benson516 0:33b75d52d2a7 727 // Using indexing
benson516 0:33b75d52d2a7 728 for (size_t i = 0; i < M; ++i){ // row in m_out
benson516 0:33b75d52d2a7 729 //
benson516 0:33b75d52d2a7 730 // cout << "#i: " << i << "\n";
benson516 0:33b75d52d2a7 731 //
benson516 0:33b75d52d2a7 732 for (size_t j = 0; j < N; ++j){ // column in m_out
benson516 0:33b75d52d2a7 733 //
benson516 0:33b75d52d2a7 734 // cout << "#j: " << j << "\n";
benson516 0:33b75d52d2a7 735 //
benson516 0:33b75d52d2a7 736 ptr_data = &(M_out.data[i][j]);
benson516 0:33b75d52d2a7 737 //
benson516 0:33b75d52d2a7 738 *ptr_data = 0.0;
benson516 0:33b75d52d2a7 739 for (size_t k = 0; k < K; ++k){
benson516 0:33b75d52d2a7 740 if (data[i][k] != 0.0 && M_right.data[k][j] != 0.0)
benson516 0:33b75d52d2a7 741 *ptr_data += data[i][k]*M_right.data[k][j];
benson516 0:33b75d52d2a7 742 }
benson516 0:33b75d52d2a7 743 }
benson516 0:33b75d52d2a7 744 }
benson516 0:33b75d52d2a7 745
benson516 0:33b75d52d2a7 746 return M_out;
benson516 0:33b75d52d2a7 747 }
benson516 0:33b75d52d2a7 748 Mat& Mat::dot(bool Transpose_left, const Mat &M_right, bool Transpose_right){ // Extended version for conbining the ability of transpose of both mtrices
benson516 0:33b75d52d2a7 749
benson516 0:33b75d52d2a7 750 // Size check
benson516 0:33b75d52d2a7 751 // mxn = (mxk)*(kxn)
benson516 0:33b75d52d2a7 752 size_t M,K,N;
benson516 0:33b75d52d2a7 753 // Left transpose
benson516 0:33b75d52d2a7 754 if (Transpose_left){
benson516 0:33b75d52d2a7 755 M = _nCol;
benson516 0:33b75d52d2a7 756 K = _nRow;
benson516 0:33b75d52d2a7 757 }else{
benson516 0:33b75d52d2a7 758 M = _nRow;
benson516 0:33b75d52d2a7 759 K = _nCol;
benson516 0:33b75d52d2a7 760 }
benson516 0:33b75d52d2a7 761 // Right transpose
benson516 0:33b75d52d2a7 762 if (Transpose_right){
benson516 0:33b75d52d2a7 763 N = M_right.nRow();
benson516 0:33b75d52d2a7 764 }else{
benson516 0:33b75d52d2a7 765 N = M_right.nCol();
benson516 0:33b75d52d2a7 766 }
benson516 0:33b75d52d2a7 767
benson516 0:33b75d52d2a7 768 // The output matrix
benson516 0:33b75d52d2a7 769 // Mat M_out(M, N);
benson516 0:33b75d52d2a7 770 static Mat M_out;
benson516 0:33b75d52d2a7 771 M_out.resize(M, N);
benson516 0:33b75d52d2a7 772
benson516 0:33b75d52d2a7 773 //
benson516 0:33b75d52d2a7 774 float* ptr_data = NULL;
benson516 0:33b75d52d2a7 775 //
benson516 0:33b75d52d2a7 776
benson516 0:33b75d52d2a7 777 // Check the conditions of transpotations
benson516 0:33b75d52d2a7 778 if(Transpose_left){
benson516 0:33b75d52d2a7 779 if(Transpose_right){ // Both transposed
benson516 0:33b75d52d2a7 780 //
benson516 0:33b75d52d2a7 781 // Using indexing
benson516 0:33b75d52d2a7 782 for (size_t i = 0; i < M; ++i){ // row in m_out
benson516 0:33b75d52d2a7 783 for (size_t j = 0; j < N; ++j){ // column in m_out
benson516 0:33b75d52d2a7 784 // M_out.data[i][j] = 0.0;
benson516 0:33b75d52d2a7 785 ptr_data = &(M_out.data[i][j]);
benson516 0:33b75d52d2a7 786 *ptr_data = 0.0;
benson516 0:33b75d52d2a7 787 //
benson516 0:33b75d52d2a7 788 for (size_t k = 0; k < K; ++k){
benson516 0:33b75d52d2a7 789 if (data[k][i] != 0.0 && M_right.data[j][k] != 0.0) // (i,k) -> (k,i), (k,j) -> (j,k)
benson516 0:33b75d52d2a7 790 *ptr_data += data[k][i]*M_right.data[j][k]; // (i,k) -> (k,i), (k,j) -> (j,k)
benson516 0:33b75d52d2a7 791 }
benson516 0:33b75d52d2a7 792 }
benson516 0:33b75d52d2a7 793 }
benson516 0:33b75d52d2a7 794 //
benson516 0:33b75d52d2a7 795 }else{ // Only left transpose
benson516 0:33b75d52d2a7 796 //
benson516 0:33b75d52d2a7 797 // Using indexing
benson516 0:33b75d52d2a7 798 for (size_t i = 0; i < M; ++i){ // row in m_out
benson516 0:33b75d52d2a7 799 for (size_t j = 0; j < N; ++j){ // column in m_out
benson516 0:33b75d52d2a7 800 // M_out.data[i][j] = 0.0;
benson516 0:33b75d52d2a7 801 ptr_data = &(M_out.data[i][j]);
benson516 0:33b75d52d2a7 802 *ptr_data = 0.0;
benson516 0:33b75d52d2a7 803 //
benson516 0:33b75d52d2a7 804 for (size_t k = 0; k < K; ++k){
benson516 0:33b75d52d2a7 805 if (data[k][i] != 0.0 && M_right.data[k][j] != 0.0) // (i,k) -> (k,i)
benson516 0:33b75d52d2a7 806 *ptr_data += data[k][i]*M_right.data[k][j]; // (i,k) -> (k,i)
benson516 0:33b75d52d2a7 807 }
benson516 0:33b75d52d2a7 808 }
benson516 0:33b75d52d2a7 809 }
benson516 0:33b75d52d2a7 810 //
benson516 0:33b75d52d2a7 811 }
benson516 0:33b75d52d2a7 812 }else{
benson516 0:33b75d52d2a7 813 if(Transpose_right){ // Only right transpose
benson516 0:33b75d52d2a7 814 //
benson516 0:33b75d52d2a7 815 // Using indexing
benson516 0:33b75d52d2a7 816 for (size_t i = 0; i < M; ++i){ // row in m_out
benson516 0:33b75d52d2a7 817 for (size_t j = 0; j < N; ++j){ // column in m_out
benson516 0:33b75d52d2a7 818 // M_out.data[i][j] = 0.0;
benson516 0:33b75d52d2a7 819 ptr_data = &(M_out.data[i][j]);
benson516 0:33b75d52d2a7 820 *ptr_data = 0.0;
benson516 0:33b75d52d2a7 821 //
benson516 0:33b75d52d2a7 822 for (size_t k = 0; k < K; ++k){
benson516 0:33b75d52d2a7 823 if (data[i][k] != 0.0 && M_right.data[j][k] != 0.0) // (k,j) -> (j,k)
benson516 0:33b75d52d2a7 824 *ptr_data += data[i][k]*M_right.data[j][k]; // (k,j) -> (j,k)
benson516 0:33b75d52d2a7 825 }
benson516 0:33b75d52d2a7 826 }
benson516 0:33b75d52d2a7 827 }
benson516 0:33b75d52d2a7 828 //
benson516 0:33b75d52d2a7 829 }else{ // Normal
benson516 0:33b75d52d2a7 830 //
benson516 0:33b75d52d2a7 831 // Using indexing
benson516 0:33b75d52d2a7 832 for (size_t i = 0; i < M; ++i){ // row in m_out
benson516 0:33b75d52d2a7 833 for (size_t j = 0; j < N; ++j){ // column in m_out
benson516 0:33b75d52d2a7 834 // M_out.data[i][j] = 0.0;
benson516 0:33b75d52d2a7 835 ptr_data = &(M_out.data[i][j]);
benson516 0:33b75d52d2a7 836 *ptr_data = 0.0;
benson516 0:33b75d52d2a7 837 //
benson516 0:33b75d52d2a7 838 for (size_t k = 0; k < K; ++k){
benson516 0:33b75d52d2a7 839 if (data[i][k] != 0.0 && M_right.data[k][j] != 0.0)
benson516 0:33b75d52d2a7 840 *ptr_data += data[i][k]*M_right.data[k][j];
benson516 0:33b75d52d2a7 841 }
benson516 0:33b75d52d2a7 842 }
benson516 0:33b75d52d2a7 843 }
benson516 0:33b75d52d2a7 844 //
benson516 0:33b75d52d2a7 845 }
benson516 0:33b75d52d2a7 846 }
benson516 0:33b75d52d2a7 847
benson516 0:33b75d52d2a7 848 return M_out;
benson516 0:33b75d52d2a7 849 }
benson516 0:33b75d52d2a7 850 // Operator overloading
benson516 0:33b75d52d2a7 851 //---------------------------------------//
benson516 0:33b75d52d2a7 852 // Member
benson516 0:33b75d52d2a7 853 // A <- b, assignment
benson516 0:33b75d52d2a7 854 Mat& Mat::operator = (float scale){ // Assign a real number as 1 by 1 matrix
benson516 0:33b75d52d2a7 855 // this->resize(1,1, scale);
benson516 0:33b75d52d2a7 856 this->resize(1,1);
benson516 0:33b75d52d2a7 857 this->data[0][0] = scale;
benson516 0:33b75d52d2a7 858 return *this;
benson516 0:33b75d52d2a7 859 }
benson516 0:33b75d52d2a7 860 Mat& Mat::operator = (const vector<float> &colVec_in){ // A = vec_in, assign as a column vector
benson516 0:33b75d52d2a7 861 this->assign(colVec_in, false); // A column vector
benson516 0:33b75d52d2a7 862 return *this;
benson516 0:33b75d52d2a7 863 }
benson516 0:33b75d52d2a7 864 Mat& Mat::operator = (const vector<vector<float> > &MatData_in){ // A = MatData_in, assign a primitive matrix directly
benson516 0:33b75d52d2a7 865 this->assign(MatData_in);
benson516 0:33b75d52d2a7 866 return *this;
benson516 0:33b75d52d2a7 867 }
benson516 0:33b75d52d2a7 868 // b <- A, get value, using implicit type conversion
benson516 0:33b75d52d2a7 869 // Note: no implicit conversion to float, since it would be ambiguous in other operator overriding
benson516 0:33b75d52d2a7 870 Mat::operator float (){ // Get the first element as float, good for a 1 by 1 matrix
benson516 0:33b75d52d2a7 871 return this->data[0][0];
benson516 0:33b75d52d2a7 872 }
benson516 0:33b75d52d2a7 873 Mat::operator std::vector<float> (){ // Get the column vector as std::vector<float> in c++
benson516 0:33b75d52d2a7 874 return this->getColVec(0);
benson516 0:33b75d52d2a7 875 }
benson516 0:33b75d52d2a7 876 // A[]
benson516 0:33b75d52d2a7 877 vector<float>& Mat::operator [] (size_t m){ // Indexing the m-th row, equivalent to data[m]
benson516 0:33b75d52d2a7 878 return this->data[m];
benson516 0:33b75d52d2a7 879 }
benson516 0:33b75d52d2a7 880 // A == B
benson516 0:33b75d52d2a7 881 bool Mat::operator == (Mat const& B){ // is_equal()
benson516 0:33b75d52d2a7 882 return this->is_equal(B);
benson516 0:33b75d52d2a7 883 }
benson516 0:33b75d52d2a7 884 // A^z, z \in Z
benson516 0:33b75d52d2a7 885 Mat& Mat::operator ^ (int power){ // A^power, this->intPower()
benson516 0:33b75d52d2a7 886 return this->intPower(power);
benson516 0:33b75d52d2a7 887 }
benson516 0:33b75d52d2a7 888
benson516 0:33b75d52d2a7 889
benson516 0:33b75d52d2a7 890 // Non-member
benson516 0:33b75d52d2a7 891 //
benson516 0:33b75d52d2a7 892 // -A
benson516 0:33b75d52d2a7 893 Mat& operator - (const Mat &A){
benson516 0:33b75d52d2a7 894 return Mat(A).times(-1.0);
benson516 0:33b75d52d2a7 895 }
benson516 0:33b75d52d2a7 896 // A + B
benson516 0:33b75d52d2a7 897 Mat& operator + (float a, const Mat &B){
benson516 0:33b75d52d2a7 898 return Mat(B).plus(a);
benson516 0:33b75d52d2a7 899 }
benson516 0:33b75d52d2a7 900 Mat& operator + (const Mat &A, float b){
benson516 0:33b75d52d2a7 901 return Mat(A).plus(b);
benson516 0:33b75d52d2a7 902 }
benson516 0:33b75d52d2a7 903 Mat& operator + (const Mat &A, const Mat &B){
benson516 0:33b75d52d2a7 904 return Mat(A).plus(B);
benson516 0:33b75d52d2a7 905 }
benson516 0:33b75d52d2a7 906 // A - B
benson516 0:33b75d52d2a7 907 Mat& operator - (float a, const Mat &B){
benson516 0:33b75d52d2a7 908 return Mat(B).minus(a, true); // (a - B)
benson516 0:33b75d52d2a7 909 }
benson516 0:33b75d52d2a7 910 Mat& operator - (const Mat &A, float b){
benson516 0:33b75d52d2a7 911 return Mat(A).minus(b);
benson516 0:33b75d52d2a7 912 }
benson516 0:33b75d52d2a7 913 Mat& operator - (const Mat &A, const Mat &B){
benson516 0:33b75d52d2a7 914 return Mat(A).minus(B);
benson516 0:33b75d52d2a7 915 }
benson516 0:33b75d52d2a7 916 // A * B
benson516 0:33b75d52d2a7 917 Mat& operator * (float a, const Mat &B){
benson516 0:33b75d52d2a7 918 return Mat(B).times(a);
benson516 0:33b75d52d2a7 919 }
benson516 0:33b75d52d2a7 920 Mat& operator * (const Mat &A, float b){
benson516 0:33b75d52d2a7 921 return Mat(A).times(b);
benson516 0:33b75d52d2a7 922 }
benson516 0:33b75d52d2a7 923 Mat& operator * (const Mat &A, const Mat &B){
benson516 0:33b75d52d2a7 924 // Matrix multiplication
benson516 0:33b75d52d2a7 925 return Mat(A).dot(B);
benson516 0:33b75d52d2a7 926 }
benson516 0:33b75d52d2a7 927 // A/B, including matrix inversion (eg, (1.0/B) <--> B^-1)
benson516 0:33b75d52d2a7 928 Mat& operator / (float a, const Mat &B){ // a*(B^-1), for that B to be inversed
benson516 0:33b75d52d2a7 929 if (a == 1.0){
benson516 0:33b75d52d2a7 930 return Mat(B).inverse();
benson516 0:33b75d52d2a7 931 }
benson516 0:33b75d52d2a7 932 return ( ( Mat(B).inverse()).times(a) );
benson516 0:33b75d52d2a7 933 }
benson516 0:33b75d52d2a7 934 Mat& operator / (const Mat &A, float b){ // A*(1/b), scalar multiplication of 1/b
benson516 0:33b75d52d2a7 935 return Mat(A).times(1.0/b);
benson516 0:33b75d52d2a7 936 }
benson516 0:33b75d52d2a7 937 Mat& operator / (const Mat &A, const Mat &B){ // A*(B^-1), for that B to be inversed
benson516 0:33b75d52d2a7 938 return ( Mat(A).dot(Mat(B).inverse()) );
benson516 0:33b75d52d2a7 939 }
benson516 0:33b75d52d2a7 940 // A += B
benson516 0:33b75d52d2a7 941 Mat& operator += (Mat &A, float b){
benson516 0:33b75d52d2a7 942 A.increase(b);
benson516 0:33b75d52d2a7 943 return A;
benson516 0:33b75d52d2a7 944 }
benson516 0:33b75d52d2a7 945 Mat& operator += (Mat &A, const Mat &B){
benson516 0:33b75d52d2a7 946 A.increase(B);
benson516 0:33b75d52d2a7 947 return A;
benson516 0:33b75d52d2a7 948 }
benson516 0:33b75d52d2a7 949 // A -= B
benson516 0:33b75d52d2a7 950 Mat& operator -= (Mat &A, float b){
benson516 0:33b75d52d2a7 951 A.decrease(b);
benson516 0:33b75d52d2a7 952 return A;
benson516 0:33b75d52d2a7 953 }
benson516 0:33b75d52d2a7 954 Mat& operator -= (Mat &A, const Mat &B){
benson516 0:33b75d52d2a7 955 A.decrease(B);
benson516 0:33b75d52d2a7 956 return A;
benson516 0:33b75d52d2a7 957 }
benson516 0:33b75d52d2a7 958 // A *= B
benson516 0:33b75d52d2a7 959 Mat& operator *= (Mat &A, float b){
benson516 0:33b75d52d2a7 960 A.scaleUp(b);
benson516 0:33b75d52d2a7 961 return A;
benson516 0:33b75d52d2a7 962 }
benson516 0:33b75d52d2a7 963 Mat& operator *= (Mat &A, const Mat &B){ // Matrix multiplication
benson516 0:33b75d52d2a7 964 A = A.dot(B);
benson516 0:33b75d52d2a7 965 return A;
benson516 0:33b75d52d2a7 966 }
benson516 0:33b75d52d2a7 967
benson516 0:33b75d52d2a7 968 //---------------------------------------//
benson516 0:33b75d52d2a7 969 // end Operator overloading