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

SHOGUN Machine Learning Toolbox - Documentation