CiftiLib
A C++ library for CIFTI-2 and CIFTI-1 files
MatrixFunctions.h
1 #ifndef __MATRIX_UTILITIES_H__
2 #define __MATRIX_UTILITIES_H__
3 
4 /*LICENSE_START*/
5 /*
6  * Copyright (c) 2014, Washington University School of Medicine
7  * All rights reserved.
8  *
9  * Redistribution and use in source and binary forms, with or without modification,
10  * are permitted provided that the following conditions are met:
11  *
12  * 1. Redistributions of source code must retain the above copyright notice,
13  * this list of conditions and the following disclaimer.
14  *
15  * 2. Redistributions in binary form must reproduce the above copyright notice,
16  * this list of conditions and the following disclaimer in the documentation
17  * and/or other materials provided with the distribution.
18  *
19  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
21  * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
23  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
24  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
25  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
26  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
27  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
28  * EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29  */
30 
31 #include <vector>
32 #include <cmath>
33 #include "stdint.h"
34 
35 using namespace std;
36 
37 //NOTE: this is not intended to be used outside of FloatMatrix.cxx, error condition is a 0x0 matrix result
38 //NOTE: if a matrix has a row shorter than the first row, expect problems. Calling checkDim will look for this, but it is not used internally, as FloatMatrix maintains this invariant
39 
40 namespace cifti {
41 
43  {
44  typedef int64_t msize_t;//NOTE: must be signed due to using -1 as a sentinel
45 
46  public:
50  template <typename T1, typename T2, typename T3, typename A>
51  static void multiply(const vector<vector<T1> > &left, const vector<vector<T2> > &right, vector<vector<T3> > &result);
52 
56  template <typename T1, typename T2, typename T3>
57  static void multiply(const vector<vector<T1> > &left, const T2 right, vector<vector<T3> > &result);
58 
62  template <typename T>
63  static void rref(vector<vector<T> > &inout);
64 
68  template <typename T>
69  static void inverse(const vector<vector<T> > &in, vector<vector<T> > &result);
70 
74  template <typename T1, typename T2, typename T3>
75  static void add(const vector<vector<T1> > &left, const vector<vector<T2> > &right, vector<vector<T3> > &result);
76 
80  template <typename T1, typename T2, typename T3>
81  static void add(const vector<vector<T1> > &left, const T2 right, vector<vector<T3> > &result);
82 
86  template <typename T1, typename T2, typename T3>
87  static void subtract(const vector<vector<T1> > &left, const vector<vector<T2> > &right, vector<vector<T3> > &result);
88 
92  template <typename T>
93  static void transpose(const vector<vector<T> > &in, vector<vector<T> > &result);
94 
98  template <typename T>
99  static bool checkDim(const vector<vector<T> > &in);
100 
104  template <typename T>
105  static void resize(const msize_t rows, const msize_t columns, vector<vector<T> > &result, bool destructive = false);
106 
110  template <typename T>
111  static void zeros(const msize_t rows, const msize_t columns, vector<vector<T> > &result);
112 
116  template <typename T>
117  static void ones(const msize_t rows, const msize_t columns, vector<vector<T> > &result);
118 
122  template <typename T>
123  static void identity(const msize_t size, vector<vector<T> > &result);
124 
128  template <typename T1, typename T2, typename T3>
129  static void horizCat(const vector<vector<T1> > &left, const vector<vector<T2> > &right, vector<vector<T3> > &result);
130 
134  template <typename T1, typename T2, typename T3>
135  static void vertCat(const vector<vector<T1> > &top, const vector<vector<T2> > &bottom, vector<vector<T3> > &result);
136 
140  template <typename T>
141  static void getChunk(const msize_t firstrow, const msize_t lastrow, const msize_t firstcol, const msize_t lastcol, const vector<vector<T> > &in, vector<vector<T> > &result);
142 
143  private:
147  template <typename T>
148  static void rref_big(vector<vector<T> > &inout);
149  };
150 
151  template <typename T1, typename T2, typename T3, typename A>
152  void MatrixFunctions::multiply(const vector<vector<T1> >& left, const vector<vector<T2> >& right, vector<vector<T3> >& result)
153  {//the stupid multiply O(n^3) - the O(n^2.78) version might not be that hard to implement with the other functions here, but not as stable
154  msize_t leftrows = (msize_t)left.size(), rightrows = (msize_t)right.size(), leftcols, rightcols;
155  vector<vector<T3> > tempstorage, *tresult = &result;//pointer because you can't change a reference
156  bool copyout = false;
157  if (&left == &result || &right == &result)
158  {
159  copyout = true;
160  tresult = &tempstorage;
161  }
162  if (leftrows && rightrows)
163  {
164  leftcols = (msize_t)left[0].size();
165  rightcols = (msize_t)right[0].size();
166  if (leftcols && rightcols && (rightrows == leftcols))
167  {
168  resize(leftrows, rightcols, (*tresult), true);//could use zeros(), but common index last lets us zero at the same time
169  msize_t i, j, k;
170  for (i = 0; i < leftrows; ++i)
171  {
172  for (j = 0; j < rightcols; ++j)
173  {
174  A accum = 0;
175  for (k = 0; k < leftcols; ++k)
176  {
177  accum += left[i][k] * right[k][j];
178  }
179  (*tresult)[i][j] = accum;
180  }
181  }
182  } else {
183  result.resize(0);
184  return;
185  }
186  } else {
187  result.resize(0);
188  return;
189  }
190  if (copyout)
191  {
192  result = tempstorage;
193  }
194  }
195 
196  template <typename T1, typename T2, typename T3>
197  void MatrixFunctions::multiply(const vector<vector<T1> > &left, const T2 right, vector<vector<T3> > &result)
198  {
199  msize_t leftrows = (msize_t)left.size(), leftcols;
200  bool doresize = true;
201  if (&left == &result)
202  {
203  doresize = false;//don't resize if an input is an output
204  }
205  if (leftrows)
206  {
207  leftcols = (msize_t)left[0].size();
208  if (leftcols)
209  {
210  if (doresize) resize(leftrows, leftcols, result, true);
211  msize_t i, j;
212  for (i = 0; i < leftrows; ++i)
213  {
214  for (j = 0; j < leftcols; ++j)
215  {
216  result[i][j] = left[i][j] * right;
217  }
218  }
219  } else {
220  result.resize(0);
221  return;
222  }
223  } else {
224  result.resize(0);
225  return;
226  }
227  }
228 
229  template<typename T>
230  void MatrixFunctions::rref_big(vector<vector<T> > &inout)
231  {
232  msize_t rows = (msize_t)inout.size(), cols;
233  if (rows > 0)
234  {
235  cols = (msize_t)inout[0].size();
236  if (cols > 0)
237  {
238  vector<msize_t> pivots(rows, -1), missingPivots;
239  msize_t i, j, k, myrow = 0;
240  msize_t pivotrow;
241  T tempval;
242  for (i = 0; i < cols; ++i)
243  {
244  if (myrow >= rows) break;//no pivots left
245  tempval = 0;
246  pivotrow = -1;
247  for (j = myrow; j < rows; ++j)
248  {//only search below for new pivot
249  if (abs(inout[j][i]) > tempval)
250  {
251  pivotrow = (msize_t)j;
252  tempval = abs(inout[j][i]);
253  }
254  }
255  if (pivotrow == -1)
256  {//naively expect linearly dependence to show as an exact zero
257  missingPivots.push_back(i);//record the missing pivot
258  continue;//move to the next column
259  }
260  inout[pivotrow].swap(inout[myrow]);//STL swap via pointers for constant time row swap
261  pivots[myrow] = i;//save the pivot location for back substitution
262  tempval = inout[myrow][i];
263  inout[myrow][i] = (T)1;
264  for (j = i + 1; j < cols; ++j)
265  {
266  inout[myrow][j] /= tempval;//divide row by pivot
267  }
268  for (j = myrow + 1; j < rows; ++j)
269  {//zero ONLY below pivot for now
270  tempval = inout[j][i];
271  inout[j][i] = (T)0;
272  for (k = i + 1; k < cols; ++k)
273  {
274  inout[j][k] -= tempval * inout[myrow][k];
275  }
276  }
277  ++myrow;//increment row on successful pivot
278  }
279  msize_t numMissing = (msize_t)missingPivots.size();
280  if (myrow > 1)//if there is only 1 pivot, there is no back substitution to do
281  {
282  msize_t lastPivotCol = pivots[myrow - 1];
283  for (i = myrow - 1; i > 0; --i)//loop through pivots, can't zero above the top pivot so exclude it
284  {
285  msize_t pivotCol = pivots[i];
286  for (j = i - 1; j >= 0; --j)//loop through rows above pivot
287  {
288  tempval = inout[j][pivotCol];
289  inout[j][pivotCol] = (T)0;//flat zero the entry above the pivot
290  for (k = numMissing - 1; k >= 0; --k)//back substitute within pivot range where pivots are missing
291  {
292  msize_t missingCol = missingPivots[k];
293  if (missingCol <= pivotCol) break;//equals will never trip, but whatever
294  inout[j][missingCol] -= tempval * inout[i][missingCol];
295  }
296  for (k = lastPivotCol + 1; k < cols; ++k)//loop through elements that are outside the pivot area
297  {
298  inout[j][k] -= tempval * inout[i][k];
299  }
300  }
301  }
302  }
303  } else {
304  inout.resize(0);
305  return;
306  }
307  } else {
308  inout.resize(0);
309  return;
310  }
311  }
312 
313  template<typename T>
314  void MatrixFunctions::rref(vector<vector<T> > &inout)
315  {
316  msize_t rows = (msize_t)inout.size(), cols;
317  if (rows)
318  {
319  cols = (msize_t)inout[0].size();
320  if (cols)
321  {
322  if (rows > 7 || cols > 7)//when the matrix has this many rows/columns, it is faster to allocate storage for tracking pivots, and back substitute
323  {
324  rref_big(inout);
325  return;
326  }
327  msize_t i, j, k, myrow = 0;
328  msize_t pivotrow;
329  T tempval;
330  for (i = 0; i < cols; ++i)
331  {
332  if (myrow >= rows) break;//no pivots left
333  tempval = 0;
334  pivotrow = -1;
335  for (j = myrow; j < rows; ++j)
336  {//only search below for new pivot
337  if (abs(inout[j][i]) > tempval)
338  {
339  pivotrow = (msize_t)j;
340  tempval = abs(inout[j][i]);
341  }
342  }
343  if (pivotrow == -1)//it may be a good idea to include a "very small value" check here, but it could mess up if used on a matrix with all values very small
344  {//naively expect linearly dependence to show as an exact zero
345  continue;//move to the next column
346  }
347  inout[pivotrow].swap(inout[myrow]);//STL swap via pointers for constant time row swap
348  tempval = inout[myrow][i];
349  inout[myrow][i] = 1;
350  for (j = i + 1; j < cols; ++j)
351  {
352  inout[myrow][j] /= tempval;//divide row by pivot
353  }
354  for (j = 0; j < myrow; ++j)
355  {//zero above pivot
356  tempval = inout[j][i];
357  inout[j][i] = 0;
358  for (k = i + 1; k < cols; ++k)
359  {
360  inout[j][k] -= tempval * inout[myrow][k];
361  }
362  }
363  for (j = myrow + 1; j < rows; ++j)
364  {//zero below pivot
365  tempval = inout[j][i];
366  inout[j][i] = 0;
367  for (k = i + 1; k < cols; ++k)
368  {
369  inout[j][k] -= tempval * inout[myrow][k];
370  }
371  }
372  ++myrow;//increment row on successful pivot
373  }
374  } else {
375  inout.resize(0);
376  return;
377  }
378  } else {
379  inout.resize(0);
380  return;
381  }
382  }
383 
384  template<typename T>
385  void MatrixFunctions::inverse(const vector<vector<T> > &in, vector<vector<T> > &result)
386  {//rref implementation, there are faster (more complicated) ways - if it isn't invertible, it will hand back something strange
387  msize_t inrows = (msize_t)in.size(), incols;
388  if (inrows)
389  {
390  incols = (msize_t)in[0].size();
391  if (incols == inrows)
392  {
393  vector<vector<T> > inter, inter2;
394  identity(incols, inter2);
395  horizCat(in, inter2, inter);
396  rref(inter);
397  getChunk(0, inrows, incols, incols * 2, inter, result);//already using a local variable, doesn't need to check for reference duplicity
398  } else {
399  result.resize(0);
400  return;
401  }
402  } else {
403  result.resize(0);
404  return;
405  }
406  }
407 
408  template <typename T1, typename T2, typename T3>
409  void MatrixFunctions::add(const vector<vector<T1> >& left, const vector<vector<T2> >& right, vector<vector<T3> >& result)
410  {
411  msize_t inrows = (msize_t)left.size(), incols;
412  bool doresize = true;
413  if (&left == &result || &right == &result)
414  {
415  doresize = false;//don't resize if an input is an output - this is ok for addition, don't need a copy
416  }
417  if (inrows)
418  {
419  incols = (msize_t)left[0].size();
420  if (inrows == (msize_t)right.size() && incols == (msize_t)right[0].size())
421  {
422  if (doresize) resize(inrows, incols, result, true);
423  for (msize_t i = 0; i < inrows; ++i)
424  {
425  for (msize_t j = 0; j < incols; ++j)
426  {
427  result[i][j] = left[i][j] + right[i][j];
428  }
429  }
430  } else {
431  result.resize(0);//use empty matrix for error condition
432  return;
433  }
434  } else {
435  result.resize(0);
436  return;
437  }
438  }
439 
440  template <typename T1, typename T2, typename T3>
441  void MatrixFunctions::add(const vector<vector<T1> >& left, const T2 right, vector<vector<T3> >& result)
442  {
443  msize_t inrows = (msize_t)left.size(), incols;
444  bool doresize = true;
445  if (&left == &result)
446  {
447  doresize = false;//don't resize if an input is an output - this is ok for addition, don't need a copy
448  }
449  if (inrows)
450  {
451  incols = (msize_t)left[0].size();
452  if (doresize) resize(inrows, incols, result, true);
453  for (msize_t i = 0; i < inrows; ++i)
454  {
455  for (msize_t j = 0; j < incols; ++j)
456  {
457  result[i][j] = left[i][j] + right;
458  }
459  }
460  } else {
461  result.resize(0);
462  return;
463  }
464  }
465 
466  template <typename T1, typename T2, typename T3>
467  void MatrixFunctions::subtract(const vector<vector<T1> >& left, const vector<vector<T2> >& right, vector<vector<T3> >& result)
468  {
469  msize_t inrows = (msize_t)left.size(), incols;
470  bool doresize = true;
471  if (&left == &result || &right == &result)
472  {
473  doresize = false;//don't resize if an input is an output
474  }
475  if (inrows)
476  {
477  incols = (msize_t)left[0].size();
478  if (inrows == (msize_t)right.size() && incols == (msize_t)right[0].size())
479  {
480  if (doresize) resize(inrows, incols, result, true);
481  for (msize_t i = 0; i < inrows; ++i)
482  {
483  for (msize_t j = 0; j < incols; ++j)
484  {
485  result[i][j] = left[i][j] - right[i][j];
486  }
487  }
488  } else {
489  result.resize(0);
490  return;
491  }
492  } else {
493  result.resize(0);
494  return;
495  }
496  }
497 
498  template<typename T>
499  void MatrixFunctions::transpose(const vector<vector<T> > &in, vector<vector<T> > &result)
500  {
501  msize_t inrows = (msize_t)in.size(), incols;
502  vector<vector<T> > tempstorage, *tresult = &result;
503  bool copyout = false;
504  if (&in == &result)
505  {
506  copyout = true;
507  tresult = &tempstorage;
508  }
509  if (inrows)
510  {
511  incols = (msize_t)in[0].size();
512  resize(incols, inrows, (*tresult), true);
513  for (msize_t i = 0; i < inrows; ++i)
514  {
515  for (msize_t j = 0; j < incols; ++j)
516  {
517  (*tresult)[j][i] = in[i][j];
518  }
519  }
520  } else {
521  result.resize(0);
522  }
523  if (copyout)
524  {
525  result = tempstorage;
526  }
527  }
528 
529  template<typename T>
530  bool MatrixFunctions::checkDim(const vector<vector<T> > &in)
531  {
532  bool ret = true;
533  msize_t rows = (msize_t)in.size(), columns;
534  if (rows)
535  {
536  columns = (msize_t)in[0].size();
537  for (msize_t i = 1; i < rows; ++i)
538  {
539  if (in[i].size() != columns)
540  {
541  ret = false;
542  }
543  }
544  }
545  return ret;
546  }
547 
548  template<typename T>
549  void MatrixFunctions::resize(const msize_t rows, const msize_t columns, vector<vector<T> >& result, bool destructive)
550  {
551  if (destructive && result.size() && ((msize_t)result.capacity() < rows || (msize_t)result[0].capacity() < columns))
552  {//for large matrices, copying to preserve contents is slow
553  result.resize(0);//not intended to dealloc, just to set number of items to copy to zero
554  }//default is nondestructive resize, copies everything
555  result.resize(rows);
556  for (msize_t i = 0; i < (const msize_t)rows; ++i)
557  {//naive method, may end up copying everything twice if both row and col resizes require realloc
558  result[i].resize(columns);
559  }
560  }
561 
562  template<typename T>
563  void MatrixFunctions::zeros(const msize_t rows, const msize_t columns, vector<vector<T> >& result)
564  {
565  resize(rows, columns, result, true);
566  for (msize_t i = 0; i < rows; ++i)
567  {
568  for (msize_t j = 0; j < columns; ++j)
569  {
570  result[i][j] = 0;//should cast to float or double fine
571  }
572  }
573  }
574 
575  template<typename T>
576  void MatrixFunctions::ones(const msize_t rows, const msize_t columns, vector<vector<T> >& result)
577  {
578  resize(rows, columns, result, true);
579  for (msize_t i = 0; i < rows; ++i)
580  {
581  for (msize_t j = 0; j < columns; ++j)
582  {
583  result[i][j] = 1;//should cast to float or double fine
584  }
585  }
586  }
587 
588  template<typename T>
589  void MatrixFunctions::identity(const msize_t size, vector<vector<T> >& result)
590  {
591  resize(size, size, result, true);
592  for (msize_t i = 0; i < (const msize_t)size; ++i)
593  {
594  for (msize_t j = 0; j < (const msize_t)size; ++j)
595  {
596  result[i][j] = ((i == j) ? 1 : 0);//ditto, forgive the ternary
597  }
598  }
599  }
600 
601  template <typename T1, typename T2, typename T3>
602  void MatrixFunctions::horizCat(const vector<vector<T1> >& left, const vector<vector<T2> >& right, vector<vector<T3> >& result)
603  {
604  msize_t inrows = (msize_t)left.size(), leftcols, rightcols;
605  vector<vector<T3> > tempstorage, *tresult = &result;
606  bool copyout = false;
607  if (&left == &result || &right == &result)
608  {
609  copyout = true;
610  tresult = &tempstorage;
611  }
612  if (inrows && inrows == (msize_t)right.size())
613  {
614  leftcols = (msize_t)left[0].size();
615  rightcols = (msize_t)right[0].size();
616  (*tresult) = left;//use STL copy to start
617  resize(inrows, leftcols + rightcols, (*tresult));//values survive nondestructive resize
618  for (msize_t i = 0; i < inrows; ++i)
619  {
620  for (msize_t j = 0; j < rightcols; ++j)
621  {
622  (*tresult)[i][j + leftcols] = right[i][j];
623  }
624  }
625  } else {
626  result.resize(0);
627  return;
628  }
629  if (copyout)
630  {
631  result = tempstorage;
632  }
633  }
634 
635  template <typename T1, typename T2, typename T3>
636  void MatrixFunctions::vertCat(const vector<vector<T1> >& top, const vector<vector<T2> >& bottom, vector<vector<T3> >& result)
637  {
638  msize_t toprows = (msize_t)top.size(), botrows = (msize_t)bottom.size(), incols;
639  vector<vector<T3> > tempstorage, *tresult = &result;
640  bool copyout = false;
641  if (&top == &result || &bottom == &result)
642  {
643  copyout = true;
644  tresult = &tempstorage;
645  }
646  if (toprows && botrows)
647  {
648  incols = (msize_t)top[0].size();
649  if (incols == (msize_t)bottom[0].size())
650  {
651  (*tresult) = top;
652  resize(toprows + botrows, incols, (*tresult));//nondestructive resize
653  for (msize_t i = 0; i < botrows; ++i)
654  {
655  for (msize_t j = 0; j < incols; ++j)
656  {
657  (*tresult)[i + toprows][j] = bottom[i][j];
658  }
659  }
660  } else {
661  result.resize(0);
662  return;
663  }
664  } else {
665  result.resize(0);
666  return;
667  }
668  if (copyout)
669  {
670  result = tempstorage;
671  }
672  }
673 
674  template<typename T>
675  void MatrixFunctions::getChunk(const msize_t firstrow, const msize_t lastrow, const msize_t firstcol, const msize_t lastcol, const vector<vector<T> >& in, vector<vector<T> >& result)
676  {
677  msize_t outrows = lastrow - firstrow;
678  msize_t outcols = lastcol - firstcol;
679  if (lastrow <= firstrow || lastcol <= firstcol || firstrow < 0 || firstcol < 0 || lastrow > (msize_t)in.size() || lastcol > (msize_t)in[0].size())
680  {
681  result.resize(0);
682  return;
683  }
684  vector<vector<T> > tempstorage, *tresult = &result;
685  bool copyout = false;
686  if (&in == &result)
687  {
688  copyout = true;
689  tresult = &tempstorage;
690  }
691  resize(outrows, outcols, (*tresult), true);
692  for (msize_t i = 0; i < outrows; ++i)
693  {
694  for (msize_t j = 0; j < outcols; ++j)
695  {
696  (*tresult)[i][j] = in[i + firstrow][j + firstcol];
697  }
698  }
699  if (copyout)
700  {
701  result = tempstorage;
702  }
703  }
704 
705 }
706 
707 #endif
708 
namespace for all CiftiLib functionality
Definition: CiftiBrainModelsMap.h:41
Definition: MatrixFunctions.h:42