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