SHOGUN  4.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
SGMatrix.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2011-2013 Heiko Strathmann
8  * Written (W) 2012 Fernando Jose Iglesias Garcia
9  * Written (W) 2010,2012 Soeren Sonnenburg
10  * Copyright (C) 2010 Berlin Institute of Technology
11  * Copyright (C) 2012 Soeren Sonnenburg
12  */
13 
14 #include <shogun/lib/config.h>
15 #include <shogun/io/SGIO.h>
16 #include <shogun/io/File.h>
17 #include <shogun/lib/SGMatrix.h>
18 #include <shogun/lib/SGVector.h>
21 #include <limits>
22 
23 #ifdef HAVE_EIGEN3
25 #endif
26 
27 namespace shogun {
28 
29 template <class T>
31 {
32  init_data();
33 }
34 
35 template <class T>
36 SGMatrix<T>::SGMatrix(bool ref_counting) : SGReferencedData(ref_counting)
37 {
38  init_data();
39 }
40 
41 template <class T>
42 SGMatrix<T>::SGMatrix(T* m, index_t nrows, index_t ncols, bool ref_counting)
43  : SGReferencedData(ref_counting), matrix(m),
44  num_rows(nrows), num_cols(ncols) { }
45 
46 template <class T>
47 SGMatrix<T>::SGMatrix(index_t nrows, index_t ncols, bool ref_counting)
48  : SGReferencedData(ref_counting), num_rows(nrows), num_cols(ncols)
49 {
50  matrix=SG_MALLOC(T, ((int64_t) nrows)*ncols);
51 }
52 
53 template <class T>
55 {
56  REQUIRE(vec.vector, "Vector not initialized!\n");
57  matrix=vec.vector;
58  num_rows=vec.vlen;
59  num_cols=1;
60 }
61 
62 template <class T>
64 : SGReferencedData(vec)
65 {
66  REQUIRE(vec.vector, "Vector not initialized!\n");
67  REQUIRE(nrows>0, "Number of rows (%d) has to be a positive integer!\n", nrows);
68  REQUIRE(ncols>0, "Number of cols (%d) has to be a positive integer!\n", ncols);
69  REQUIRE(vec.vlen==nrows*ncols, "Number of elements in the matrix (%d) must "
70  "be the same as the number of elements in the vector (%d)!\n",
71  nrows*ncols, vec.vlen);
72  matrix=vec.vector;
73  num_rows=nrows;
74  num_cols=ncols;
75 }
76 
77 template <class T>
79 {
80  copy_data(orig);
81 }
82 
83 #ifdef HAVE_EIGEN3
84 template <class T>
86 : SGReferencedData(false), matrix(mat.data()),
87  num_rows(mat.rows()), num_cols(mat.cols())
88 {
89 
90 }
91 
92 template <class T>
94 {
95  return EigenMatrixXtMap(matrix, num_rows, num_cols);
96 }
97 #endif
98 
99 template <class T>
101 {
102  unref();
103 }
104 
105 template <class T>
107 {
108  if (num_rows!=other.num_rows || num_cols!=other.num_cols)
109  return false;
110 
111  if (matrix!=other.matrix)
112  return false;
113 
114  return true;
115 }
116 
117 template <class T>
119 {
120  if (num_rows!=other.num_rows || num_cols!=other.num_cols)
121  return false;
122 
123  for (int64_t i=0; i<int64_t(num_rows)*num_cols; ++i)
124  {
125  if (matrix[i]!=other.matrix[i])
126  return false;
127  }
128 
129  return true;
130 }
131 
132 template <class T>
133 void SGMatrix<T>::set_const(T const_elem)
134 {
135  for (int64_t i=0; i<int64_t(num_rows)*num_cols; i++)
136  matrix[i]=const_elem ;
137 }
138 
139 template <class T>
141 {
142  if (matrix && (int64_t(num_rows)*num_cols))
143  set_const(0);
144 }
145 
146 template <>
148 {
149  if (matrix && (int64_t(num_rows)*num_cols))
150  set_const(complex128_t(0.0));
151 }
152 
153 template <class T>
155 {
156  if (num_rows!=num_cols)
157  return false;
158  for (int i=0; i<num_rows; ++i)
159  {
160  for (int j=i+1; j<num_cols; ++j)
161  {
162  if (matrix[j*num_rows+i]!=matrix[i*num_rows+j])
163  return false;
164  }
165  }
166  return true;
167 }
168 
169 template <>
171 {
172  if (num_rows!=num_cols)
173  return false;
174  for (int i=0; i<num_rows; ++i)
175  {
176  for (int j=i+1; j<num_cols; ++j)
177  {
178  if (!CMath::fequals<float32_t>(matrix[j*num_rows+i],
179  matrix[i*num_rows+j], FLT_EPSILON))
180  return false;
181  }
182  }
183  return true;
184 }
185 
186 template <>
188 {
189  if (num_rows!=num_cols)
190  return false;
191  for (int i=0; i<num_rows; ++i)
192  {
193  for (int j=i+1; j<num_cols; ++j)
194  {
195  if (!CMath::fequals<float64_t>(matrix[j*num_rows+i],
196  matrix[i*num_rows+j], DBL_EPSILON))
197  return false;
198  }
199  }
200  return true;
201 }
202 
203 template <>
205 {
206  if (num_rows!=num_cols)
207  return false;
208  for (int i=0; i<num_rows; ++i)
209  {
210  for (int j=i+1; j<num_cols; ++j)
211  {
212  if (!CMath::fequals<floatmax_t>(matrix[j*num_rows+i],
213  matrix[i*num_rows+j], LDBL_EPSILON))
214  return false;
215  }
216  }
217  return true;
218 }
219 
220 template <>
222 {
223  if (num_rows!=num_cols)
224  return false;
225  for (int i=0; i<num_rows; ++i)
226  {
227  for (int j=i+1; j<num_cols; ++j)
228  {
229  if (!(CMath::fequals<float64_t>(matrix[j*num_rows+i].real(),
230  matrix[i*num_rows+j].real(), DBL_EPSILON) &&
231  CMath::fequals<float64_t>(matrix[j*num_rows+i].imag(),
232  matrix[i*num_rows+j].imag(), DBL_EPSILON)))
233  return false;
234  }
235  }
236  return true;
237 }
238 
239 template <class T>
241 {
242  T max=matrix[0];
243  for (int64_t i=1; i<int64_t(num_rows)*num_cols; ++i)
244  {
245  if (matrix[i]>max)
246  max=matrix[i];
247  }
248 
249  return max;
250 }
251 
252 template <>
254 {
255  SG_SERROR("SGMatrix::max_single():: Not supported for complex128_t\n");
256  return complex128_t(0.0);
257 }
258 
259 template <class T>
261 {
262  return SGMatrix<T>(clone_matrix(matrix, num_rows, num_cols),
263  num_rows, num_cols);
264 }
265 
266 template <class T>
267 T* SGMatrix<T>::clone_matrix(const T* matrix, int32_t nrows, int32_t ncols)
268 {
269  T* result = SG_MALLOC(T, int64_t(nrows)*ncols);
270  for (int64_t i=0; i<int64_t(nrows)*ncols; i++)
271  result[i]=matrix[i];
272 
273  return result;
274 }
275 
276 template <class T>
278  T*& matrix, int32_t& num_feat, int32_t& num_vec)
279 {
280  /* this should be done in-place! Heiko */
281  T* transposed=SG_MALLOC(T, int64_t(num_vec)*num_feat);
282  for (int64_t i=0; i<num_vec; i++)
283  {
284  for (int64_t j=0; j<num_feat; j++)
285  transposed[i+j*num_vec]=matrix[i*num_feat+j];
286  }
287 
288  SG_FREE(matrix);
289  matrix=transposed;
290 
291  CMath::swap(num_feat, num_vec);
292 }
293 
294 template <class T>
295 void SGMatrix<T>::create_diagonal_matrix(T* matrix, T* v,int32_t size)
296 {
297  for(int64_t i=0;i<size;i++)
298  {
299  for(int64_t j=0;j<size;j++)
300  {
301  if(i==j)
302  matrix[j*size+i]=v[i];
303  else
304  matrix[j*size+i]=0;
305  }
306  }
307 }
308 
309 template <class T>
311  float64_t* mat, int32_t cols, int32_t rows)
312 {
313  float64_t trace=0;
314  for (int64_t i=0; i<rows; i++)
315  trace+=mat[i*cols+i];
316  return trace;
317 }
318 
319 template <class T>
320 T* SGMatrix<T>::get_row_sum(T* matrix, int32_t m, int32_t n)
321 {
322  T* rowsums=SG_CALLOC(T, n);
323 
324  for (int64_t i=0; i<n; i++)
325  {
326  for (int64_t j=0; j<m; j++)
327  rowsums[i]+=matrix[j+i*m];
328  }
329  return rowsums;
330 }
331 
332 template <class T>
333 T* SGMatrix<T>::get_column_sum(T* matrix, int32_t m, int32_t n)
334 {
335  T* colsums=SG_CALLOC(T, m);
336 
337  for (int64_t i=0; i<n; i++)
338  {
339  for (int64_t j=0; j<m; j++)
340  colsums[j]+=matrix[j+i*m];
341  }
342  return colsums;
343 }
344 
345 template <class T>
347 {
348  center_matrix(matrix, num_rows, num_cols);
349 }
350 
351 template <class T>
352 void SGMatrix<T>::center_matrix(T* matrix, int32_t m, int32_t n)
353 {
354  float64_t num_data=n;
355 
356  T* colsums=get_column_sum(matrix, m,n);
357  T* rowsums=get_row_sum(matrix, m,n);
358 
359  for (int32_t i=0; i<m; i++)
360  colsums[i]/=num_data;
361  for (int32_t j=0; j<n; j++)
362  rowsums[j]/=num_data;
363 
364  T s=SGVector<T>::sum(rowsums, n)/num_data;
365 
366  for (int64_t i=0; i<n; i++)
367  {
368  for (int64_t j=0; j<m; j++)
369  matrix[i*m+j]+=s-colsums[j]-rowsums[i];
370  }
371 
372  SG_FREE(rowsums);
373  SG_FREE(colsums);
374 }
375 
376 template <class T>
378 {
379  /* compute "row" sums (which I would call column sums), i.e. sum of all
380  * elements in a fixed column */
381  T* means=get_row_sum(matrix, num_rows, num_cols);
382 
383  /* substract column mean from every corresponding entry */
384  for (int64_t i=0; i<num_cols; ++i)
385  {
386  means[i]/=num_rows;
387  for (int64_t j=0; j<num_rows; ++j)
388  matrix[i*num_rows+j]-=means[i];
389  }
390 
391  SG_FREE(means);
392 }
393 
394 template<class T> void SGMatrix<T>::display_matrix(const char* name) const
395 {
396  display_matrix(matrix, num_rows, num_cols, name);
397 }
398 
399 template <class T>
401  const SGMatrix<T> matrix, const char* name,
402  const char* prefix)
403 {
404  matrix.display_matrix();
405 }
406 
407 template <>
409  const bool* matrix, int32_t rows, int32_t cols, const char* name,
410  const char* prefix)
411 {
412  ASSERT(rows>=0 && cols>=0)
413  SG_SPRINT("%s%s=[\n", prefix, name)
414  for (int64_t i=0; i<rows; i++)
415  {
416  SG_SPRINT("%s[", prefix)
417  for (int64_t j=0; j<cols; j++)
418  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i] ? 1 : 0,
419  j==cols-1? "" : ",");
420  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
421  }
422  SG_SPRINT("%s]\n", prefix)
423 }
424 
425 template <>
427  const char* matrix, int32_t rows, int32_t cols, const char* name,
428  const char* prefix)
429 {
430  ASSERT(rows>=0 && cols>=0)
431  SG_SPRINT("%s%s=[\n", prefix, name)
432  for (int64_t i=0; i<rows; i++)
433  {
434  SG_SPRINT("%s[", prefix)
435  for (int64_t j=0; j<cols; j++)
436  SG_SPRINT("%s\t%c%s", prefix, matrix[j*rows+i],
437  j==cols-1? "" : ",");
438  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
439  }
440  SG_SPRINT("%s]\n", prefix)
441 }
442 
443 template <>
445  const int8_t* matrix, int32_t rows, int32_t cols, const char* name,
446  const char* prefix)
447 {
448  ASSERT(rows>=0 && cols>=0)
449  SG_SPRINT("%s%s=[\n", prefix, name)
450  for (int64_t i=0; i<rows; i++)
451  {
452  SG_SPRINT("%s[", prefix)
453  for (int64_t j=0; j<cols; j++)
454  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
455  j==cols-1? "" : ",");
456  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
457  }
458  SG_SPRINT("%s]\n", prefix)
459 }
460 
461 template <>
463  const uint8_t* matrix, int32_t rows, int32_t cols, const char* name,
464  const char* prefix)
465 {
466  ASSERT(rows>=0 && cols>=0)
467  SG_SPRINT("%s%s=[\n", prefix, name)
468  for (int64_t i=0; i<rows; i++)
469  {
470  SG_SPRINT("%s[", prefix)
471  for (int64_t j=0; j<cols; j++)
472  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
473  j==cols-1? "" : ",");
474  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
475  }
476  SG_SPRINT("%s]\n", prefix)
477 }
478 
479 template <>
481  const int16_t* matrix, int32_t rows, int32_t cols, const char* name,
482  const char* prefix)
483 {
484  ASSERT(rows>=0 && cols>=0)
485  SG_SPRINT("%s%s=[\n", prefix, name)
486  for (int64_t i=0; i<rows; i++)
487  {
488  SG_SPRINT("%s[", prefix)
489  for (int64_t j=0; j<cols; j++)
490  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
491  j==cols-1? "" : ",");
492  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
493  }
494  SG_SPRINT("%s]\n", prefix)
495 }
496 
497 template <>
499  const uint16_t* matrix, int32_t rows, int32_t cols, const char* name,
500  const char* prefix)
501 {
502  ASSERT(rows>=0 && cols>=0)
503  SG_SPRINT("%s%s=[\n", prefix, name)
504  for (int64_t i=0; i<rows; i++)
505  {
506  SG_SPRINT("%s[", prefix)
507  for (int64_t j=0; j<cols; j++)
508  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
509  j==cols-1? "" : ",");
510  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
511  }
512  SG_SPRINT("%s]\n", prefix)
513 }
514 
515 
516 template <>
518  const int32_t* matrix, int32_t rows, int32_t cols, const char* name,
519  const char* prefix)
520 {
521  ASSERT(rows>=0 && cols>=0)
522  SG_SPRINT("%s%s=[\n", prefix, name)
523  for (int64_t i=0; i<rows; i++)
524  {
525  SG_SPRINT("%s[", prefix)
526  for (int64_t j=0; j<cols; j++)
527  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
528  j==cols-1? "" : ",");
529  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
530  }
531  SG_SPRINT("%s]\n", prefix)
532 }
533 
534 template <>
536  const uint32_t* matrix, int32_t rows, int32_t cols, const char* name,
537  const char* prefix)
538 {
539  ASSERT(rows>=0 && cols>=0)
540  SG_SPRINT("%s%s=[\n", prefix, name)
541  for (int64_t i=0; i<rows; i++)
542  {
543  SG_SPRINT("%s[", prefix)
544  for (int64_t j=0; j<cols; j++)
545  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
546  j==cols-1? "" : ",");
547  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
548  }
549  SG_SPRINT("%s]\n", prefix)
550 }
551 template <>
553  const int64_t* matrix, int32_t rows, int32_t cols, const char* name,
554  const char* prefix)
555 {
556  ASSERT(rows>=0 && cols>=0)
557  SG_SPRINT("%s%s=[\n", prefix, name)
558  for (int64_t i=0; i<rows; i++)
559  {
560  SG_SPRINT("%s[", prefix)
561  for (int64_t j=0; j<cols; j++)
562  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
563  j==cols-1? "" : ",");
564  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
565  }
566  SG_SPRINT("%s]\n", prefix)
567 }
568 
569 template <>
571  const uint64_t* matrix, int32_t rows, int32_t cols, const char* name,
572  const char* prefix)
573 {
574  ASSERT(rows>=0 && cols>=0)
575  SG_SPRINT("%s%s=[\n", prefix, name)
576  for (int64_t i=0; i<rows; i++)
577  {
578  SG_SPRINT("%s[", prefix)
579  for (int64_t j=0; j<cols; j++)
580  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
581  j==cols-1? "" : ",");
582  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
583  }
584  SG_SPRINT("%s]\n", prefix)
585 }
586 
587 template <>
589  const float32_t* matrix, int32_t rows, int32_t cols, const char* name,
590  const char* prefix)
591 {
592  ASSERT(rows>=0 && cols>=0)
593  SG_SPRINT("%s%s=[\n", prefix, name)
594  for (int64_t i=0; i<rows; i++)
595  {
596  SG_SPRINT("%s[", prefix)
597  for (int64_t j=0; j<cols; j++)
598  SG_SPRINT("%s\t%.18g%s", prefix, (float) matrix[j*rows+i],
599  j==cols-1? "" : ",");
600  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
601  }
602  SG_SPRINT("%s]\n", prefix)
603 }
604 
605 template <>
607  const float64_t* matrix, int32_t rows, int32_t cols, const char* name,
608  const char* prefix)
609 {
610  ASSERT(rows>=0 && cols>=0)
611  SG_SPRINT("%s%s=[\n", prefix, name)
612  for (int64_t i=0; i<rows; i++)
613  {
614  SG_SPRINT("%s[", prefix)
615  for (int64_t j=0; j<cols; j++)
616  SG_SPRINT("%s\t%.18g%s", prefix, (double) matrix[j*rows+i],
617  j==cols-1? "" : ",");
618  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
619  }
620  SG_SPRINT("%s]\n", prefix)
621 }
622 
623 template <>
625  const floatmax_t* matrix, int32_t rows, int32_t cols, const char* name,
626  const char* prefix)
627 {
628  ASSERT(rows>=0 && cols>=0)
629  SG_SPRINT("%s%s=[\n", prefix, name)
630  for (int64_t i=0; i<rows; i++)
631  {
632  SG_SPRINT("%s[", prefix)
633  for (int64_t j=0; j<cols; j++)
634  SG_SPRINT("%s\t%.18g%s", prefix, (double) matrix[j*rows+i],
635  j==cols-1? "" : ",");
636  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
637  }
638  SG_SPRINT("%s]\n", prefix)
639 }
640 
641 template <>
643  const complex128_t* matrix, int32_t rows, int32_t cols, const char* name,
644  const char* prefix)
645 {
646  ASSERT(rows>=0 && cols>=0)
647  SG_SPRINT("%s%s=[\n", prefix, name)
648  for (int64_t i=0; i<rows; i++)
649  {
650  SG_SPRINT("%s[", prefix)
651  for (int64_t j=0; j<cols; j++)
652  SG_SPRINT("%s\t(%.18g+i%.18g)%s", prefix, matrix[j*rows+i].real(),
653  matrix[j*rows+i].imag(), j==cols-1? "" : ",");
654  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
655  }
656  SG_SPRINT("%s]\n", prefix)
657 }
658 
659 template <>
661 {
663  return SGMatrix<char>();
664 }
665 
666 template <>
668 {
669  SGMatrix<int8_t> I(size, size);
670  for (index_t i=0; i<size; ++i)
671  {
672  for (index_t j=0; j<size; ++j)
673  I(i,j)=i==j ? scale : 0.0;
674  }
675 
676  return I;
677 }
678 
679 template <>
681 {
682  SGMatrix<uint8_t> I(size, size);
683  for (index_t i=0; i<size; ++i)
684  {
685  for (index_t j=0; j<size; ++j)
686  I(i,j)=i==j ? scale : 0.0;
687  }
688 
689  return I;
690 }
691 
692 template <>
694 {
695  SGMatrix<bool> I(size, size);
696  for (index_t i=0; i<size; ++i)
697  {
698  for (index_t j=0; j<size; ++j)
699  I(i,j)=i==j ? scale : (!scale);
700  }
701 
702  return I;
703 }
704 
705 template <>
707 {
708  SGMatrix<int16_t> I(size, size);
709  for (index_t i=0; i<size; ++i)
710  {
711  for (index_t j=0; j<size; ++j)
712  I(i,j)=i==j ? scale : 0.0;
713  }
714 
715  return I;
716 }
717 
718 template <>
720 {
721  SGMatrix<uint16_t> I(size, size);
722  for (index_t i=0; i<size; ++i)
723  {
724  for (index_t j=0; j<size; ++j)
725  I(i,j)=i==j ? scale : 0.0;
726  }
727 
728  return I;
729 }
730 
731 template <>
733 {
734  SGMatrix<int32_t> I(size, size);
735  for (index_t i=0; i<size; ++i)
736  {
737  for (index_t j=0; j<size; ++j)
738  I(i,j)=i==j ? scale : 0.0;
739  }
740 
741  return I;
742 }
743 
744 template <>
746 {
747  SGMatrix<uint32_t> I(size, size);
748  for (index_t i=0; i<size; ++i)
749  {
750  for (index_t j=0; j<size; ++j)
751  I(i,j)=i==j ? scale : 0.0;
752  }
753 
754  return I;
755 }
756 
757 template <>
759 {
760  SGMatrix<int64_t> I(size, size);
761  for (index_t i=0; i<size; ++i)
762  {
763  for (index_t j=0; j<size; ++j)
764  I(i,j)=i==j ? scale : 0.0;
765  }
766 
767  return I;
768 }
769 
770 template <>
772 {
773  SGMatrix<uint64_t> I(size, size);
774  for (index_t i=0; i<size; ++i)
775  {
776  for (index_t j=0; j<size; ++j)
777  I(i,j)=i==j ? scale : 0.0;
778  }
779 
780  return I;
781 }
782 
783 template <>
785 {
786  SGMatrix<float32_t> I(size, size);
787  for (index_t i=0; i<size; ++i)
788  {
789  for (index_t j=0; j<size; ++j)
790  I(i,j)=i==j ? scale : 0.0;
791  }
792 
793  return I;
794 }
795 
796 template <>
798 {
799  SGMatrix<float64_t> I(size, size);
800  for (index_t i=0; i<size; ++i)
801  {
802  for (index_t j=0; j<size; ++j)
803  I(i,j)=i==j ? scale : 0.0;
804  }
805 
806  return I;
807 }
808 
809 template <>
811 {
812  SGMatrix<floatmax_t> I(size, size);
813  for (index_t i=0; i<size; ++i)
814  {
815  for (index_t j=0; j<size; ++j)
816  I(i,j)=i==j ? scale : 0.0;
817  }
818 
819  return I;
820 }
821 
822 template <>
824 {
825  SGMatrix<complex128_t> I(size, size);
826  for (index_t i=0; i<size; ++i)
827  {
828  for (index_t j=0; j<size; ++j)
829  I(i,j)=i==j ? scale : complex128_t(0.0);
830  }
831 
832  return I;
833 }
834 
835 //Howto construct the pseudo inverse (from "The Matrix Cookbook")
836 //
837 //Assume A does not have full rank, i.e. A is n \times m and rank(A) = r < min(n;m).
838 //
839 //The matrix A+ known as the pseudo inverse is unique and does always exist.
840 //
841 //The pseudo inverse A+ can be constructed from the singular value
842 //decomposition A = UDV^T , by A^+ = V(D+)U^T.
843 
844 #ifdef HAVE_LAPACK
845 template <class T>
847  float64_t* matrix, int32_t rows, int32_t cols, float64_t* target)
848 {
849  if (!target)
850  target=SG_MALLOC(float64_t, rows*cols);
851 
852  char jobu='A';
853  char jobvt='A';
854  int m=rows; /* for calling external lib */
855  int n=cols; /* for calling external lib */
856  int lda=m; /* for calling external lib */
857  int ldu=m; /* for calling external lib */
858  int ldvt=n; /* for calling external lib */
859  int info=-1; /* for calling external lib */
860  int32_t lsize=CMath::min((int32_t) m, (int32_t) n);
861  double* s=SG_MALLOC(double, lsize);
862  double* u=SG_MALLOC(double, m*m);
863  double* vt=SG_MALLOC(double, n*n);
864 
865  wrap_dgesvd(jobu, jobvt, m, n, matrix, lda, s, u, ldu, vt, ldvt, &info);
866  ASSERT(info==0)
867 
868  for (int64_t i=0; i<n; i++)
869  {
870  for (int64_t j=0; j<lsize; j++)
871  vt[i*n+j]=vt[i*n+j]/s[j];
872  }
873 
874  cblas_dgemm(CblasColMajor, CblasTrans, CblasTrans, m, n, m, 1.0, vt, ldvt, u, ldu, 0, target, m);
875 
876  SG_FREE(u);
877  SG_FREE(vt);
878  SG_FREE(s);
879 
880  return target;
881 }
882 
884 template <class T>
886 {
887  ASSERT(matrix.num_cols==matrix.num_rows)
888  int32_t* ipiv = SG_MALLOC(int32_t, matrix.num_cols);
889  clapack_dgetrf(CblasColMajor,matrix.num_cols,matrix.num_cols,matrix.matrix,matrix.num_cols,ipiv);
890  clapack_dgetri(CblasColMajor,matrix.num_cols,matrix.matrix,matrix.num_cols,ipiv);
891  SG_FREE(ipiv);
892 }
893 
894 template <class T>
896 {
897  if (matrix.num_rows!=matrix.num_cols)
898  {
899  SG_SERROR("SGMatrix::compute_eigenvectors(SGMatrix<float64_t>): matrix"
900  " rows and columns are not equal!\n");
901  }
902 
903  /* use reference counting for SGVector */
904  SGVector<float64_t> result(NULL, 0, true);
905  result.vlen=matrix.num_rows;
906  result.vector=compute_eigenvectors(matrix.matrix, matrix.num_rows,
907  matrix.num_rows);
908  return result;
909 }
910 
911 template <class T>
912 double* SGMatrix<T>::compute_eigenvectors(double* matrix, int n, int m)
913 {
914  ASSERT(n == m)
915 
916  char V='V';
917  char U='U';
918  int info;
919  int ord=n;
920  int lda=n;
921  double* eigenvalues=SG_CALLOC(float64_t, n+1);
922 
923  // lapack sym matrix eigenvalues+vectors
924  wrap_dsyev(V, U, ord, matrix, lda,
925  eigenvalues, &info);
926 
927  if (info!=0)
928  SG_SERROR("DSYEV failed with code %d\n", info)
929 
930  return eigenvalues;
931 }
932 
933 template <class T>
934 void SGMatrix<T>::compute_few_eigenvectors(double* matrix_, double*& eigenvalues, double*& eigenvectors,
935  int n, int il, int iu)
936 {
937  eigenvalues = SG_MALLOC(double, n);
938  eigenvectors = SG_MALLOC(double, (iu-il+1)*n);
939  int status = 0;
940  wrap_dsyevr('V','U',n,matrix_,n,il,iu,eigenvalues,eigenvectors,&status);
941  ASSERT(status==0)
942 }
943 
944 #endif //HAVE_LAPACK
945 
946 template <class T>
949  bool transpose_A, bool transpose_B, float64_t scale)
950 {
951  /* these variables store size of transposed matrices*/
952  index_t cols_A=transpose_A ? A.num_rows : A.num_cols;
953  index_t rows_A=transpose_A ? A.num_cols : A.num_rows;
954  index_t rows_B=transpose_B ? B.num_cols : B.num_rows;
955  index_t cols_B=transpose_B ? B.num_rows : B.num_cols;
956 
957  /* do a dimension check */
958  if (cols_A!=rows_B)
959  {
960  SG_SERROR("SGMatrix::matrix_multiply(): Dimension mismatch: "
961  "A(%dx%d)*B(%dx%D)\n", rows_A, cols_A, rows_B, cols_B);
962  }
963 
964  /* allocate result matrix */
965  SGMatrix<float64_t> C(rows_A, cols_B);
966  C.zero();
967 #ifdef HAVE_LAPACK
968  /* multiply */
969  cblas_dgemm(CblasColMajor,
970  transpose_A ? CblasTrans : CblasNoTrans,
971  transpose_B ? CblasTrans : CblasNoTrans,
972  rows_A, cols_B, cols_A, scale,
973  A.matrix, A.num_rows, B.matrix, B.num_rows,
974  0.0, C.matrix, C.num_rows);
975 #else
976  /* C(i,j) = scale * \Sigma A(i,k)*B(k,j) */
977  for (int32_t i=0; i<rows_A; i++)
978  {
979  for (int32_t j=0; j<cols_B; j++)
980  {
981  for (int32_t k=0; k<cols_A; k++)
982  {
983  float64_t x1=transpose_A ? A(k,i):A(i,k);
984  float64_t x2=transpose_B ? B(j,k):B(k,j);
985  C(i,j)+=x1*x2;
986  }
987 
988  C(i,j)*=scale;
989  }
990  }
991 #endif //HAVE_LAPACK
992 
993  return C;
994 }
995 
996 template<class T>
998  index_t num_cols, SGMatrix<T> pre_allocated)
999 {
1000  SGMatrix<T> result;
1001 
1002  /* evtl use pre-allocated space */
1003  if (pre_allocated.matrix)
1004  {
1005  result=pre_allocated;
1006 
1007  /* check dimension consistency */
1008  if (pre_allocated.num_rows!=num_rows ||
1009  pre_allocated.num_cols!=num_cols)
1010  {
1011  SG_SERROR("SGMatrix<T>::get_allocated_matrix(). Provided target"
1012  "matrix dimensions (%dx%d) do not match passed data "
1013  "dimensions (%dx%d)!\n", pre_allocated.num_rows,
1014  pre_allocated.num_cols, num_rows, num_cols);
1015  }
1016  }
1017  else
1018  {
1019  /* otherwise, allocate space */
1020  result=SGMatrix<T>(num_rows, num_cols);
1021  }
1022 
1023  return result;
1024 }
1025 
1026 template<class T>
1028 {
1029  matrix=((SGMatrix*)(&orig))->matrix;
1030  num_rows=((SGMatrix*)(&orig))->num_rows;
1031  num_cols=((SGMatrix*)(&orig))->num_cols;
1032 }
1033 
1034 template<class T>
1036 {
1037  matrix=NULL;
1038  num_rows=0;
1039  num_cols=0;
1040 }
1041 
1042 template<class T>
1044 {
1045  SG_FREE(matrix);
1046  matrix=NULL;
1047  num_rows=0;
1048  num_cols=0;
1049 }
1050 
1051 template<class T>
1053 {
1054  ASSERT(loader)
1055  unref();
1056 
1058  SGMatrix<T> mat;
1059  loader->get_matrix(mat.matrix, mat.num_rows, mat.num_cols);
1060  copy_data(mat);
1061  copy_refcount(mat);
1062  ref();
1064 }
1065 
1066 template<>
1068 {
1069  SG_SERROR("SGMatrix::load():: Not supported for complex128_t\n");
1070 }
1071 
1072 template<class T>
1074 {
1075  ASSERT(writer)
1077  writer->set_matrix(matrix, num_rows, num_cols);
1079 }
1080 
1081 template<>
1083 {
1084  SG_SERROR("SGMatrix::save():: Not supported for complex128_t\n");
1085 }
1086 
1087 template<class T>
1089 {
1090  SGVector<T> rowv(num_cols);
1091  for (int64_t i = 0; i < num_cols; i++)
1092  {
1093  rowv[i] = matrix[i*num_rows+row];
1094  }
1095  return rowv;
1096 }
1097 
1098 template<class T>
1100 {
1101  index_t diag_vlen=CMath::min(num_cols, num_rows);
1102  SGVector<T> diag(diag_vlen);
1103 
1104  for (int64_t i=0; i<diag_vlen; i++)
1105  {
1106  diag[i]=matrix[i*num_rows+i];
1107  }
1108 
1109  return diag;
1110 }
1111 
1112 template class SGMatrix<bool>;
1113 template class SGMatrix<char>;
1114 template class SGMatrix<int8_t>;
1115 template class SGMatrix<uint8_t>;
1116 template class SGMatrix<int16_t>;
1117 template class SGMatrix<uint16_t>;
1118 template class SGMatrix<int32_t>;
1119 template class SGMatrix<uint32_t>;
1120 template class SGMatrix<int64_t>;
1121 template class SGMatrix<uint64_t>;
1122 template class SGMatrix<float32_t>;
1123 template class SGMatrix<float64_t>;
1124 template class SGMatrix<floatmax_t>;
1125 template class SGMatrix<complex128_t>;
1126 }
void display_matrix(const char *name="matrix") const
Definition: SGMatrix.cpp:394
void wrap_dsyevr(char jobz, char uplo, int n, double *a, int lda, int il, int iu, double *eigenvalues, double *eigenvectors, int *info)
Definition: lapack.cpp:307
#define SG_RESET_LOCALE
Definition: SGIO.h:86
std::complex< float64_t > complex128_t
Definition: common.h:67
virtual void set_matrix(const bool *matrix, int32_t num_feat, int32_t num_vec)
Definition: File.cpp:126
int32_t index_t
Definition: common.h:62
#define REQUIRE(x,...)
Definition: SGIO.h:206
index_t num_cols
Definition: SGMatrix.h:378
#define SG_SNOTIMPLEMENTED
Definition: SGIO.h:198
virtual ~SGMatrix()
Definition: SGMatrix.cpp:100
bool operator==(SGMatrix< T > &other)
Definition: SGMatrix.cpp:106
#define SG_SET_LOCALE_C
Definition: SGIO.h:85
index_t num_rows
Definition: SGMatrix.h:376
shogun matrix
Definition: Parameter.h:26
void wrap_dgesvd(char jobu, char jobvt, int m, int n, double *a, int lda, double *sing, double *u, int ldu, double *vt, int ldvt, int *info)
Definition: lapack.cpp:253
index_t vlen
Definition: SGVector.h:494
#define SG_SPRINT(...)
Definition: SGIO.h:180
#define ASSERT(x)
Definition: SGIO.h:201
shogun vector
Definition: Parameter.h:28
virtual void get_matrix(bool *&matrix, int32_t &num_feat, int32_t &num_vec)
Definition: File.cpp:109
shogun reference count managed data
double float64_t
Definition: common.h:50
long double floatmax_t
Definition: common.h:51
A File access base class.
Definition: File.h:34
void wrap_dsyev(char jobz, char uplo, int n, double *a, int lda, double *w, int *info)
Definition: lapack.cpp:235
virtual void copy_data(const SGReferencedData &orig)
Definition: SGMatrix.cpp:1027
float float32_t
Definition: common.h:49
bool is_symmetric()
Definition: SGMatrix.cpp:154
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
#define SG_SERROR(...)
Definition: SGIO.h:179
void scale(Matrix A, Matrix B, typename Matrix::Scalar alpha)
Definition: Core.h:93
Matrix::Scalar max(Matrix m)
Definition: Redux.h:66
bool equals(SGMatrix< T > &other)
Definition: SGMatrix.cpp:118
void set_const(T const_elem)
Definition: SGMatrix.cpp:133
virtual void init_data()
Definition: SGMatrix.cpp:1035

SHOGUN Machine Learning Toolbox - Documentation