1 #ifndef MATRIX_H_INCLUDED 2 #define MATRIX_H_INCLUDED 3 4 #include <vector> 5 #include <list> 6 #include <assert.h> 7 #include <algorithm> 8 #include "vektor.h" 9 #include "printer.h" 10 11 using namespace std; 12 13 template <class typ> class Matrix{ 14 //public: 15 int width,height; 16 // vector<Vektor<typ> > rows; 17 vector<typ> data; 18 public: getHeight()19 inline int getHeight()const{return height;}; getWidth()20 inline int getWidth()const{return width;}; Matrix(const Matrix & a)21 Matrix(const Matrix &a):data(a.data),width(a.getWidth()),height(a.getHeight()){ 22 } Matrix(int height_,int width_)23 Matrix(int height_, int width_):data(height_*width_),height(height_),width(width_){ 24 assert(height>=0); 25 assert(width>=0); 26 for(int i=0;i<height_*width_;i++)data[i]=0; 27 }; ~Matrix()28 ~Matrix(){ 29 }; Matrix()30 Matrix():width(0),height(0){ 31 }; column(int i)32 Vektor<typ> column(int i)const 33 { 34 assert(i>=0); 35 assert(i<getWidth()); 36 Vektor<typ> ret(getHeight()); 37 for(int j=0;j<getHeight();j++)ret[j]=data[j*getWidth()+i]; 38 return ret; 39 } transposed()40 Matrix transposed()const 41 { 42 Matrix ret(getWidth(),getHeight()); 43 for(int i=0;i<getWidth();i++) 44 for(int j=0;j<getHeight();j++) 45 ret.data[i*getHeight()+j]=data[j*getWidth()+i]; 46 return ret; 47 } identity(int n)48 static Matrix identity(int n) 49 { 50 Matrix m(n,n); 51 for(int i=0;i<n;i++)m.data[i*(n+1)]=typ(1); 52 return m; 53 } append(Matrix const & m)54 void append(Matrix const &m) 55 { 56 assert(getWidth()==m.getWidth()); 57 assert(&m!=this);//If the two matrices are the same object the following will not work: 58 data.insert(data.end(),m.data.begin(),m.data.end()); 59 height+=m.height; 60 } appendRow(Vektor<typ> const & r)61 void appendRow(Vektor<typ> const &r) 62 { 63 assert(r.size()==width); 64 int oldsize=data.size(); 65 data.resize(data.size()+width); 66 for(int i=0;i<width;i++)data[oldsize+i]=r.v[i]; 67 height++; 68 } 69 /* void prependRow(Vektor<typ> const &r) 70 { 71 assert(r.size()==width); 72 rows.front_back(r); 73 height++; 74 }*/ setRow(int i,IntegerVector const & v)75 void setRow(int i,IntegerVector const &v) 76 { 77 assert(v.size()==getWidth()); 78 assert(i>=0 && i<getHeight()); 79 for(int j=0;j<getWidth();j++)data[i*getWidth()+j]=v.v[j]; 80 } getRow(int i)81 IntegerVector getRow(int i)const 82 { 83 IntegerVector ret(getWidth()); 84 for(int j=0;j<getWidth();j++)ret.v[j]=data[i*getWidth()+j]; 85 return ret; 86 } getRows()87 IntegerVectorList getRows()const 88 { 89 IntegerVectorList ret; 90 for(int i=0;i<height;i++)ret.push_back(getRow(i)); 91 return ret; 92 } dotRow(IntegerVector const & v,int i)93 typ dotRow(IntegerVector const &v, int i)const 94 { 95 assert(v.size()==getWidth()); 96 typ s=0; 97 for(int j=0;j<getWidth();j++) 98 s+=data[j+i*getWidth()]*v.v[j]; 99 return s; 100 } dotRowLong(IntegerVector const & v,int i)101 typ dotRowLong(IntegerVector const &v, int i)const 102 { 103 assert(v.size()==getWidth()); 104 int64 s=0; 105 for(int j=0;j<getWidth();j++) 106 s+=((int64)(data[j+i*getWidth()]))*v.v[j]; 107 return s; 108 } vectormultiply(IntegerVector const & v)109 IntegerVector vectormultiply(IntegerVector const &v)const 110 { 111 assert(v.size()==width); 112 IntegerVector ret(height); 113 for(int i=0;i<height;i++) 114 ret[i]=dotRow(v,i); 115 return ret; 116 } 117 /** 118 * Decides if v is in the kernel of the matrix. 119 */ inKernel(IntegerVector const & v)120 bool inKernel(IntegerVector const &v)const 121 { 122 assert(v.size()==width); 123 for(int i=0;i<height;i++) 124 if(dotRowLong(v,i)!=0)return false; 125 return true; 126 } 127 inline friend Matrix operator*(typ s, const Matrix& q) 128 { 129 Matrix p=q; 130 for(int i=0;i<q.data.size();i++)p.data[i]=s*q.data[i]; 131 return p; 132 } tropicalProduct(const Matrix & a,const Matrix & b)133 friend Matrix tropicalProduct(const Matrix& a, const Matrix& b) 134 { 135 int neutral=-888888888; 136 assert(a.width==b.height); 137 Matrix ret(a.height,b.width); 138 for(int i=0;i<a.height;i++) 139 for(int j=0;j<b.width;j++) 140 { 141 int sum=neutral; 142 for(int k=0;k<b.height;k++) 143 { 144 if(a[i][k]+b[k][j]>sum)sum=a[i][k]+b[k][j]; 145 } 146 ret[i][j]=sum; 147 } 148 return ret; 149 } 150 friend Matrix operator*(const Matrix& a, const Matrix& b) 151 { 152 assert(a.width==b.height); 153 Matrix ret(b.width,a.height); 154 for(int i=0;i<b.width;i++) 155 ret.setRow(i,a.vectormultiply(b.column(i))); 156 return ret.transposed(); 157 } 158 Matrix operator-()const 159 { 160 Matrix ret(height,width); 161 for(int i=0;i<height*width;i++) 162 ret.data[i]=-((*this).data[i]); 163 return ret; 164 } 165 /* template<class T> 166 Matrix<T>(const Matrix<T>& c):v(c.size()){ 167 for(int i=0;i<size();i++)v[i]=typ(c[i]);} 168 */ 169 170 /** 171 Returns the specified submatrix. The endRow and endColumn are not included. 172 */ submatrix(int startRow,int startColumn,int endRow,int endColumn)173 Matrix submatrix(int startRow, int startColumn, int endRow, int endColumn)const 174 { 175 assert(startRow>=0); 176 assert(startColumn>=0); 177 assert(endRow>=startRow); 178 assert(endColumn>=startColumn); 179 assert(endRow<=height); 180 assert(endColumn<=width); 181 Matrix ret(endRow-startRow,endColumn-startColumn); 182 for(int i=startRow;i<endRow;i++) 183 for(int j=startColumn;j<endColumn;j++) 184 // ret[i-startRow][j-startColumn]=rows[i][j]; 185 ret.data[(i-startRow)*ret.width+j-startColumn]=data[i*width+j]; 186 return ret; 187 } submatrixRows(int startRow,int endRow)188 Matrix submatrixRows(int startRow, int endRow)const 189 { 190 return submatrix(startRow,0,endRow,getWidth()); 191 } 192 /** 193 Returns the specified submatrix. The height is the height of *this, but the returned matrix only contains 194 those columns i, for which subset[i]==1. All other entries of subset must be zero. 195 */ submatrixColumnSubsetBoolean(IntegerVector const & subset)196 Matrix submatrixColumnSubsetBoolean(IntegerVector const &subset)const 197 { 198 assert(subset.size()==this->getWidth()); 199 Matrix ret(getHeight(),subset.sum()); 200 for(int i=0;i<getHeight();i++) 201 { 202 int J=0; 203 for(int j=0;j<subset.size();j++) 204 if(subset[j])ret[i][J++]=data[i*getWidth()+j]; 205 } 206 return ret; 207 } 208 209 class RowRef; 210 class const_RowRef{ 211 int rowNumM; 212 Matrix const &matrix; 213 friend class RowRef; 214 public: const_RowRef(const Matrix & matrix_,int rowNum_)215 inline const_RowRef(const Matrix &matrix_, int rowNum_): 216 rowNumM(rowNum_*matrix_.width), 217 matrix(matrix_) 218 { 219 } 220 inline typ const &operator[](int j)const 221 { 222 assert(j>=0); 223 assert(j<matrix.width); 224 return matrix.data[rowNumM+j]; 225 } UNCHECKEDACCESS(int j)226 inline typ const &UNCHECKEDACCESS(int j)const 227 { 228 return matrix.data[rowNumM+j]; 229 } toVector()230 const Vektor<typ> toVector()const 231 { 232 Vektor<typ> ret(matrix.width); 233 for(int j=0;j<matrix.width;j++) 234 ret[j]=matrix.data[rowNumM+j]; 235 return ret; 236 } 237 // operator const IntegerVector()const 238 operator Vektor<typ>()const 239 { 240 return toVector(); 241 } 242 bool operator==(Vektor<typ> const &b)const 243 { 244 return toVector()==b; 245 } dot(Vektor<typ> const & b)246 typ dot(Vektor<typ> const &b)const 247 { 248 return dot(toVector(),b); 249 } 250 Vektor<typ> operator-()const 251 { 252 return -toVector(); 253 } 254 /* bool isZero()const 255 { 256 for(int j=0;j<matrix.width;j++)if(!matrix.isZero(matrix.data[matrix.width*rowNum+j]))return false; 257 return true; 258 }*/ 259 }; 260 261 class RowRef{ 262 int rowNumM; 263 Matrix &matrix; 264 public: RowRef(Matrix & matrix_,int rowNum_)265 inline RowRef(Matrix &matrix_, int rowNum_): 266 rowNumM(rowNum_*matrix_.width), 267 matrix(matrix_) 268 { 269 } 270 inline typ &operator[](int j) 271 { 272 assert(j>=0); 273 assert(j<matrix.width); 274 return matrix.data[rowNumM+j]; 275 } UNCHECKEDACCESS(int j)276 inline typ &UNCHECKEDACCESS(int j) 277 { 278 return matrix.data[rowNumM+j]; 279 } 280 RowRef &operator=(Vektor<typ> const &v) 281 { 282 assert(v.size()==matrix.width); 283 for(int j=0;j<matrix.width;j++) 284 matrix.data[rowNumM+j]=v[j]; 285 286 return *this; 287 } 288 RowRef &operator=(RowRef const &v) 289 { 290 assert(v.matrix.width==matrix.width); 291 for(int j=0;j<matrix.width;j++) 292 matrix.data[rowNumM+j]=v.matrix.data[v.rowNumM+j]; 293 294 return *this; 295 } 296 RowRef &operator+=(Vektor<typ> const &v) 297 { 298 assert(v.size()==matrix.width); 299 for(int j=0;j<matrix.width;j++) 300 matrix.data[rowNumM+j]+=v.v[j]; 301 302 return *this; 303 } 304 RowRef &operator+=(RowRef const &v) 305 { 306 assert(v.matrix.width==matrix.width); 307 for(int j=0;j<matrix.width;j++) 308 matrix.data[rowNumM+j]+=v.matrix.data[v.rowNumM+j]; 309 310 return *this; 311 } 312 RowRef &operator+=(const_RowRef const &v) 313 { 314 assert(v.matrix.width==matrix.width); 315 for(int j=0;j<matrix.width;j++) 316 matrix.data[rowNumM+j]+=v.matrix.data[v.rowNumM+j]; 317 318 return *this; 319 } 320 RowRef &operator=(const_RowRef const &v) 321 { 322 assert(v.matrix.width==matrix.width); 323 for(int j=0;j<matrix.width;j++) 324 matrix.data[rowNumM+j]=v.matrix.data[v.rowNumM+j]; 325 326 return *this; 327 } toVector()328 const Vektor<typ> toVector()const 329 { 330 Vektor<typ> ret(matrix.width); 331 for(int j=0;j<matrix.width;j++) 332 ret[j]=matrix.data[rowNumM+j]; 333 return ret; 334 } 335 // operator const IntegerVector()const 336 operator Vektor<typ>()const 337 { 338 return toVector(); 339 } dot(Vektor<typ> const & b)340 typ dot(Vektor<typ> const &b)const 341 { 342 return dot(toVector(),b); 343 } 344 /* Vector toVector() 345 { 346 Vector ret(matrix.width); 347 for(int j=0;j<matrix.width;j++) 348 ret[j]=matrix.data[matrix.width*rowNum+j]; 349 return ret; 350 } 351 void set(Vector const &v) 352 { 353 assert(v.size()==matrix.width); 354 for(int j=0;j<matrix.width;j++) 355 matrix.data[matrix.width*rowNum+j]=v[j]; 356 } 357 bool isZero()const 358 { 359 for(int j=0;j<matrix.width;j++)if(!matrix.isZero(matrix.data[matrix.width*rowNum+j]))return false; 360 return true; 361 }*/ 362 }; 363 364 // const Vektor<typ>& operator[](int n)const{assert(n>=0 && n<getHeight());return (rows[n]);} 365 // Vektor<typ>& operator[](int n){assert(n>=0 && n<getHeight());return (rows[n]);} 366 // Bugfix for gcc4.5 (passing assertion to the above operator: 367 // Vektor<typ>& operator[](int n){if(!(n>=0 && n<getHeight()))(*(const Matrix<typ>*)(this))[n];return (rows[n]);} 368 369 inline RowRef operator[](int i) 370 { 371 assert(i>=0); 372 assert(i<height); 373 return RowRef(*this,i); 374 } 375 inline const_RowRef operator[](int i)const 376 { 377 assert(i>=0); 378 assert(i<height); 379 return const_RowRef(*this,i); 380 } 381 382 /** 383 Takes two matrices with the same number of columns and construct 384 a new matrix which has the rows of the matrix top on the top and 385 the rows of the matrix bottom at the bottom. The return value is 386 the constructed matrix. 387 */ combineOnTop(Matrix const & top,Matrix const & bottom)388 friend Matrix combineOnTop(Matrix const &top, Matrix const &bottom) 389 { 390 assert(bottom.getWidth()==top.getWidth()); 391 Matrix ret(top.getHeight()+bottom.getHeight(),top.getWidth()); 392 for(int i=0;i<top.getHeight()*top.getWidth();i++)ret.data[i]=top.data[i]; 393 for(int i=0;i<bottom.getHeight()*top.getWidth();i++)ret.data[i+top.getHeight()*top.getWidth()]=bottom.data[i]; 394 395 return ret; 396 } 397 /** 398 Takes two matrices with the same number of rows and construct 399 a new matrix which has the columns of the matrix left on the left and 400 the columns of the matrix right on the right. The return value is 401 the constructed matrix. 402 */ combineLeftRight(Matrix const & left,Matrix const & right)403 friend Matrix combineLeftRight(Matrix const &left, Matrix const &right) 404 { 405 assert(left.getHeight()==right.getHeight()); 406 Matrix ret(left.getHeight(),left.getWidth()+right.getWidth()); 407 for(int i=0;i<left.getHeight();i++) 408 { 409 for(int j=0;j<left.getWidth();j++)ret.data[i*ret.getWidth()+j]=left.data[i*left.getWidth()+j]; 410 for(int j=0;j<right.getWidth();j++)ret.data[i*ret.getWidth()+j+left.getWidth()]=right.data[i*right.getWidth()+j]; 411 } 412 return ret; 413 } 414 friend Printer& operator<<(Printer &p, Matrix const &m) 415 { 416 Vektor<int> widths(m.getWidth()); 417 for(int i=0;i<m.getHeight();i++) 418 for(int j=0;j<m.getWidth();j++) 419 { 420 stringstream s; 421 s<<m[i][j]; 422 if(s.str().length()>widths[j])widths[j]=s.str().length(); 423 } 424 stringstream s; 425 s<<"{"; 426 for(int i=0;i<m.getHeight();i++) 427 { 428 if(i)s<<","; 429 s<<"\n"; 430 s<<"("; 431 for(int j=0;j<m.getWidth();j++) 432 { 433 if(j)s<<","; 434 stringstream s2; 435 s2<<m[i][j]; 436 for(int k=s2.str().length();k<widths[j];k++)s<<" "; 437 s<<s2.str(); 438 } 439 s<<")"; 440 } 441 s<<"}\n"; 442 443 p<<s.str(); 444 return p; 445 } UNCHECKEDACCESS(int i,int j)446 const typ& UNCHECKEDACCESS(int i,int j)const{ 447 /* assert(i>=0); 448 assert(i<height); 449 assert(j>=0); 450 assert(j<width);*/ 451 return (data[i*getWidth()+j]);} UNCHECKEDACCESS(int i,int j)452 typ& UNCHECKEDACCESS(int i,int j){ 453 /* assert(i>=0); 454 assert(i<height); 455 assert(j>=0); 456 assert(j<width);*/ 457 return (data[i*getWidth()+j]);} UNCHECKEDACCESSRESTRICT(int i,int j)458 typ* __restrict__ UNCHECKEDACCESSRESTRICT(int i,int j){return & (data[i*getWidth()+j]);} swapColumns(int a,int b)459 void swapColumns(int a, int b) 460 { 461 if(a!=b) 462 for(int i=0;i<getHeight();i++) 463 swap(data[i*getWidth()+a],data[i*getWidth()+b]); 464 } 465 }; 466 467 typedef Matrix<int> IntegerMatrix; 468 typedef Matrix<double> FloatMatrix; 469 470 IntegerMatrix rowsToIntegerMatrix(IntegerVectorList const &rows, int width=-1);//width specifies the matrix width. If no width is specied the width is found by looking at the length of the rows. The function "asserts" if the length of the rows does not match the matrix size or if the width was not specified and could not be read off from the rows. 471 IntegerMatrix rowToIntegerMatrix(IntegerVector const &row); 472 473 FloatMatrix integerToFloatMatrix(IntegerMatrix const &m); 474 IntegerVector flattenMatrix(IntegerMatrix const &m); 475 int rank(IntegerMatrix const &m); 476 477 #endif 478