libStatGen Software  1
MathMatrix.cpp
1 /*
2  * Copyright (C) 2010 Regents of the University of Michigan
3  *
4  * This program is free software: you can redistribute it and/or modify
5  * it under the terms of the GNU General Public License as published by
6  * the Free Software Foundation, either version 3 of the License, or
7  * (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program. If not, see <http://www.gnu.org/licenses/>.
16  */
17 
18 #include "MathMatrix.h"
19 #include "MathVector.h"
20 #include "MathConstant.h"
21 #include "Sort.h"
22 #include "Error.h"
23 
24 #include <string.h>
25 #include <math.h>
26 #include <stdio.h>
27 
28 int Matrix::alloc = 2;
29 
30 Matrix::~Matrix()
31 {
32  // printf("Deleting Matrix %s...\n", (const char *) label);
33 
34  for (int i=0; i<size; i++)
35  delete data[i];
36 
37  if (size)
38  delete [] data;
39 
40  if (extraSize)
41  delete [] extras;
42 }
43 
44 void Matrix::Init()
45 {
46  rows = cols = extraSize = size = 0;
47  data = NULL;
48  extras = NULL;
49  label = "[Matrix]";
50 }
51 
52 void Matrix::SetLabel(const char * name)
53 {
54  label = '[';
55  label += name;
56  label += ']';
57 }
58 
59 void Matrix::Dimension(int m, int n)
60 {
61  if (n == cols && m == rows)
62  return;
63 
64  if (n > extraSize)
65  {
66  int newSize = (n + alloc) / alloc * alloc;
67  ColumnExtras * newExtras = new ColumnExtras [newSize];
68 
69  if (extras != NULL)
70  for (int i = 0; i < extraSize; i++)
71  newExtras[i] = extras[i];
72 
73  if (extraSize)
74  delete [] extras;
75 
76  extraSize = newSize;
77  extras = newExtras;
78  }
79 
80  if (m > size)
81  {
82  int newSize = (m + alloc) / alloc * alloc;
83  Vector ** newData = new Vector * [newSize];
84 
85  if (data != NULL)
86  for (int i = 0; i < size; i++)
87  newData[i] = data[i];
88 
89  for (int i = size; i < newSize; i++)
90  newData[i] = new Vector(n);
91 
92  if (size)
93  delete [] data;
94 
95  size = newSize;
96  data = newData;
97  }
98 
99  if (cols != n)
100  for (int i = 0; i < rows; i++)
101  data[i]->Dimension(n);
102 
103  if (rows != m)
104  for (int i = rows; i < m; i++)
105  data[i]->Dimension(n);
106 
107  rows = m;
108  cols = n;
109 }
110 
111 void Matrix::Dimension(int m, int n, double value)
112 {
113  int originalRows = rows;
114  int originalColumns = cols;
115 
116  Dimension(m, n);
117 
118  if (rows > originalRows)
119  for (int i = originalRows; i < rows; i++)
120  data[i]->Set(value);
121 
122  if (cols > originalColumns)
123  for (int i = 0; i < originalRows; i++)
124  for (int j = originalColumns; j < cols; j++)
125  data[i]->data[j] = value;
126 }
127 
128 void Matrix::Zero()
129 {
130  for (int i = 0; i < rows; i++)
131  for (int j = 0; j < cols; j++)
132  (*(data[i]))[j] = 0.0;
133 }
134 
135 void Matrix::Identity()
136 {
137  if (rows != cols)
138  error("Matrix.Identity - Identity matrices must be square");
139 
140  for (int i = 0; i < rows; i++)
141  for (int j = 0; j < cols; j++)
142  if (i == j)
143  (*(data[i]))[j] = 1.0;
144  else
145  (*(data[i]))[j] = 0.0;
146 }
147 
148 void Matrix::Set(double k)
149 {
150  for (int i = 0; i < rows; i++)
151  for (int j = 0; j < cols; j++)
152  (*(data[i]))[j] = k;
153 }
154 
155 void Matrix::Negate()
156 {
157  for (int i = 0; i < rows; i++)
158  for (int j = 0; j < cols; j++)
159  (*(data[i]))[j] = -(*(data[i]))[j];
160 }
161 
162 void Matrix::Copy(const Matrix & m)
163 {
164  Dimension(m.rows, m.cols);
165 
166  if (m.data != NULL)
167  for (int i = 0; i < rows; i++)
168  for (int j = 0; j < cols; j++)
169  (*this)[i][j] = m[i][j];
170 }
171 
172 void Matrix::Transpose(const Matrix & m)
173 {
174  Dimension(m.cols, m.rows);
175 
176  for (int i = 0; i < rows; i++)
177  for (int j = 0; j < cols; j++)
178  (*(data[i]))[j] = m[j][i];
179 }
180 
181 void Matrix::Add(double k)
182 {
183  for (int i = 0; i < rows; i++)
184  for (int j = 0; j < cols; j++)
185  (*(data[i]))[j] += k;
186 }
187 
188 void Matrix::Multiply(double k)
189 {
190  for (int i = 0; i < rows; i++)
191  for (int j = 0; j < cols; j++)
192  (*(data[i]))[j] *= k;
193 }
194 
195 void Matrix::Add(const Matrix & m)
196 {
197  if ((rows != m.rows) && (cols != m.cols))
198  error("Matrix.Add - Attempted to add incompatible matrices\n"
199  "Matrices - %s [%d, %d] + %s [%d, %d]\n",
200  (const char *) label, rows, cols,
201  (const char *) m.label, m.rows, m.cols);
202 
203  for (int i = 0; i < rows; i++)
204  for (int j = 0; j < cols; j++)
205  (*(data[i]))[j] += m[i][j];
206 }
207 
208 void Matrix::AddMultiple(double k, const Matrix & m)
209 {
210  if ((rows != m.rows) && (cols != m.cols))
211  error("Matrix.AddMultiple - Attempted to add incompatible matrices\n"
212  "Matrices - %s [%d, %d] + k * %s [%d, %d]\n",
213  (const char *) label, rows, cols,
214  (const char *) m.label, m.rows, m.cols);
215 
216  for (int i = 0; i < rows; i++)
217  for (int j = 0; j < cols; j++)
218  (*(data[i]))[j] += k * m[i][j];
219 }
220 
221 
222 void Matrix::Product(const Matrix & l, const Matrix & r)
223 {
224  if (l.cols != r.rows)
225  error("Matrix.Multiply - Attempted to multiply incompatible matrices\n"
226  "Matrices - %s [%d, %d] + %s [%d, %d]\n",
227  (const char *) l.label, l.rows, l.cols,
228  (const char *) r.label, r.rows, r.cols);
229 
230  Dimension(l.rows, r.cols);
231  Zero();
232 
233  for (int k = 0; k < l.cols; k++)
234  for (int i = 0; i < rows; i++)
235  for (int j = 0; j < cols; j++)
236  (*(data[i]))[j] += l[i][k] * r[k][j];
237 }
238 
239 void Matrix::AddRows(double k, int r1, int r2)
240 {
241  Vector v(*(data[r1]));
242 
243  v.Multiply(k);
244  data[r2]->Add(v);
245 }
246 
247 void Matrix::MultiplyRow(int r1, double k)
248 {
249  data[r1]->Multiply(k);
250 }
251 
252 void Matrix::AddRows(int r1, int r2)
253 {
254  data[r2]->Add(*(data[r1]));
255 }
256 
257 void Matrix::Reduce(double tol)
258 {
259  double pivot;
260  int pivotr = 0; // Initializing pivotr is not necessary, but avoids warnings
261  int r = 0; // the row we are currently reducing
262 
263  for (int j = 0; j < cols; j++)
264  {
265  if (r > rows)
266  return;
267 
268  pivot = 0.0;
269  for (int i = r; i < rows; i++)
270  if (fabs((*this)[i][j]) > pivot)
271  {
272  pivot = fabs((*this)[i][j]);
273  pivotr = i;
274  }
275 
276  if (pivot <= tol)
277  {
278  for (int i = r; i < rows; i++)
279  (*this)[i][j] = 0.0;
280  continue;
281  }
282 
283  SwapRows(pivotr, r);
284 
285  double scale = (*this)[r][j];
286 
287  (*this)[r][j] = 1.0;
288  for (int k = j+1; k < cols; k++)
289  (*this)[r][k] /= scale;
290 
291  for (int i = r + 1; r < rows; i++)
292  {
293  scale = (*this)[i][j];
294  (*this)[i][j] = 0.0;
295  for (int k = j+1; k < cols; k++)
296  (*this)[i][k] -= (*this)[r][k] * scale;
297  }
298 
299  r++;
300  }
301 }
302 
303 void Matrix::DeleteRow(int r)
304 {
305  Vector * temp = data[r];
306 
307  for (int i = r + 1; i < rows; i++)
308  data[i-1] = data[i];
309 
310  data[rows - 1] = temp;
311  rows--;
312 }
313 
314 void Matrix::DeleteColumn(int c)
315 {
316  for (int i = 0; i < rows; i++)
317  data[i] -> DeleteDimension(c);
318 
319  for (int i = c + 1; i < cols; i++)
320  extras[i-1] = extras[i];
321 
322  cols--;
323 }
324 
325 void Matrix::SwapColumns(int c1, int c2)
326 {
327  double temp;
328 
329  for (int i = 0; i < rows; i++)
330  {
331  temp = (*data[i])[c1];
332  (*data[i])[c1] = (*data[i])[c2];
333  (*data[i])[c2] = temp;
334  }
335 
336  extras[c1].Swap(extras[c2]);
337 }
338 
339 void Matrix::Read(FILE * f)
340 {
341  int r, c;
342  char buffer[100];
343  int numItems = 0;
344 
345  numItems = fscanf(f, " %s =", buffer);
346  if(numItems != 1) { }
347  buffer[strlen(buffer) - 1] = 0;
348  SetLabel(buffer);
349 
350  numItems = fscanf(f, " [ %d x %d ]", &r, &c);
351  if(numItems != 2) { }
352  Dimension(r, c);
353 
354  for (int c = 0; c < cols; c++)
355  {
356  numItems = fscanf(f, " %s", buffer);
357  if(numItems != 1) { }
358  SetColumnLabel(c, buffer);
359  }
360 
361  for (int r = 0; r < rows; r++)
362  for (int c = 0; c < cols; c++)
363  {
364  numItems = fscanf(f, " %lf", &((*this)[r][c]));
365  if(numItems != 1) { }
366  }
367 }
368 
369 
370 void Matrix::Print(FILE * f, int r, int c)
371 {
372  if (r == -1 || r > rows) r = rows;
373  if (c == -1 || c > cols) c = cols;
374 
375  char dimensions[30];
376 
377  sprintf(dimensions, "[%d x %d]", r, c);
378 
379  int columnZero = label.Length() > 15 ? label.Length() : 15;
380 
381  fprintf(f, "\n%*s =\n%*s ", columnZero, (const char *) label,
382  columnZero, dimensions);
383 
384  int * precision = new int [c + 1];
385  int * width = new int [c + 1];
386 
387  for (int j = 0; j < c; j++)
388  {
389  precision[j] = extras[j].GetPrecision();
390  width[j] = extras[j].GetWidth();
391  fprintf(f, "%*s ", width[j], (const char *) extras[j].label);
392  }
393 
394  for (int i = 0; i < r; i++)
395  {
396  fprintf(f, "\n%*s ", columnZero, (const char *) data[i]->label);
397  for (int j = 0; j < c; j++)
398  fprintf(f, "%*.*f ", width[j], precision[j], (*this)[i][j]);
399  }
400 
401  fprintf(f, "\n");
402 
403  delete [] precision;
404  delete [] width;
405 }
406 
407 void Matrix::CopyLabels(Matrix & M)
408 {
409  for (int i = 0; i < rows; i++)
410  if (i < M.rows)
411  data[i]->SetLabel(M[i].label);
412 
413  for (int i = 0; i < cols; i++)
414  if (i < M.cols)
415  SetColumnLabel(i, M.GetColumnLabel(i));
416 }
417 
418 // ColumnExtras class
419 //
420 
421 void ColumnExtras::Init()
422 {
423  label = "column";
424  dirty = true;
425  precision = 3;
426  width = 7;
427 }
428 
429 ColumnExtras::~ColumnExtras()
430 { }
431 
432 void ColumnExtras::SetLabel(const char * name)
433 {
434  label = name;
435 }
436 
437 int ColumnExtras::GetWidth()
438 {
439  if (dirty)
440  {
441  if (precision + 2 > width)
442  width = precision + 2;
443  if (label.Length() > width)
444  width = label.Length();
445  dirty = false;
446  }
447  return width;
448 }
449 
450 void ColumnExtras::Copy(ColumnExtras & c)
451 {
452  width = c.width;
453  precision = c.precision;
454  dirty = c.dirty;
455  label = c.label;
456 }
457 
458 #define SWAP(a,b) {int swap=(a); (a)=(b); (b)=swap;}
459 #define SWAPBOOL(a,b) {bool swap=(a); (a)=(b); (b)=swap;}
460 
461 void ColumnExtras::Swap(ColumnExtras & c)
462 {
463  SWAP(c.width, width);
464  SWAP(c.precision, precision);
465  SWAPBOOL(c.dirty, dirty);
466  c.label.Swap(label);
467 }
468 
469 int Matrix::CompareRows(Vector ** row1, Vector ** row2)
470 {
471  if ((**row1)[0] < (**row2)[0]) return -1;
472  if ((**row1)[0] > (**row2)[0]) return 1;
473  return 0;
474 }
475 
476 void Matrix::Sort()
477 {
478  QuickSort(data, rows, sizeof(Vector *), COMPAREFUNC CompareRows);
479 }
480 
481 bool Matrix::operator == (const Matrix & rhs) const
482 {
483  if (rhs.rows != rows || rhs.cols != cols) return false;
484 
485  for (int i = 0; i < rows; i++)
486  if ((*this)[i] != rhs[i])
487  return false;
488  return true;
489 }
490 
491 void Matrix::StackBottom(const Matrix & m)
492 {
493  if (m.cols != cols)
494  error("Attempted to stack matrices with different number of columns");
495 
496  int end = rows;
497 
498  Dimension(rows + m.rows, cols);
499 
500  for (int i = 0; i < m.rows; i++)
501  *(data[i + end]) = m[i];
502 }
503 
504 void Matrix::StackLeft(const Matrix & m)
505 {
506  if (m.rows != rows)
507  error("Attempted to stack matrics with different numbers of rows");
508 
509  for (int i = 0; i < rows; i++)
510  data[i]->Stack(m[i]);
511 
512  Dimension(rows, cols + m.cols);
513 }
514 
515 void Matrix::Swap(Matrix & m)
516 {
517  label.Swap(m.label);
518 
519  ColumnExtras * tmpExtras = extras;
520  extras = m.extras;
521  m.extras = tmpExtras;
522 
523  int swap;
524  swap = rows;
525  rows = m.rows;
526  m.rows = swap;
527  swap = cols;
528  cols = m.cols;
529  m.cols = swap;
530  swap = size;
531  size = m.size;
532  m.size = swap;
533  swap = extraSize;
534  extraSize = m.extraSize;
535  m.extraSize = swap;
536 
537  Vector ** tmpData = data;
538  data = m.data;
539  m.data = tmpData;
540 }
541 
542 double Matrix::Min() const
543 {
544  if (rows == 0 || cols == 0)
545  return 0.0;
546 
547  double minimum = data[0]->Min();
548 
549  for (int i = 1; i < rows; i++)
550  minimum = min(data[i]->Min(), minimum);
551 
552  return minimum;
553 }
554 
555 double Matrix::Max() const
556 {
557  if (rows == 0 || cols == 0)
558  return 0.0;
559 
560  double maximum = data[0]->Max();
561 
562  for (int i = 1; i < rows; i++)
563  maximum = max(data[i]->Max(), maximum);
564 
565  return maximum;
566 }
567 
568 double Matrix::Mean() const
569 {
570  if (rows == 0 || cols == 0)
571  return 0.0;
572 
573  double sum = data[0]->Sum();
574 
575  for (int i = 1; i < rows; i++)
576  sum += data[i]->Sum();
577 
578  return sum / (rows * cols);
579 }
580 
581 double Matrix::SafeMin() const
582 {
583  double lo = (rows > 0 && cols > 0) ? _NAN_ : 0.0;
584 
585  int i, j;
586 
587  for (i = 0; i < rows; i++)
588  {
589  for (j = 0; j < cols; j++)
590  if (data[i]->data[j] != _NAN_)
591  {
592  lo = data[i]->data[j];
593  break;
594  }
595  if (j != cols) break;
596  }
597 
598  for (; i < rows; i++, j = 0)
599  for (; j < cols; j++)
600  if (data[i]->data[j] < lo && data[i]->data[j] != _NAN_)
601  lo = data[i]->data[j];
602 
603  return lo;
604 }
605 
606 double Matrix::SafeMax() const
607 {
608  double hi = (rows > 0 && cols > 0) ? _NAN_ : 0.0;
609 
610  int i, j;
611 
612  for (i = 0; i < rows; i++)
613  {
614  for (j = 0; j < cols; j++)
615  if (data[i]->data[j] != _NAN_)
616  {
617  hi = data[i]->data[j];
618  break;
619  }
620  if (j != cols) break;
621  }
622 
623  for (; i < rows; i++, j = 0)
624  for (; j < cols; j++)
625  if (data[i]->data[j] > hi && data[i]->data[j] != _NAN_)
626  hi = data[i]->data[j];
627 
628  return hi;
629 }
630 
631 double Matrix::SafeMean() const
632 {
633  double sum = 0.0;
634  int count = 0;
635 
636  for (int i = 0; i < rows; i++)
637  for (int j = 0; j < cols; j++)
638  if ((*this)[i][j] != _NAN_)
639  {
640  sum += (*this)[i][j];
641  count ++;
642  }
643 
644  return (count) ? sum / count : 0.0;
645 }
646 
647 int Matrix::SafeCount() const
648 {
649  int total = 0;
650 
651  for (int i = 0; i < rows; i++)
652  total += data[i]->SafeCount();
653 
654  return total;
655 }
656 
657 void Matrix::PrintUpper(FILE * f, int r, int c, bool print_diag)
658 {
659  int columnZero;
660  int * precision = NULL, * width = NULL; // Initialization avoids compiler warnings
661 
662  SetupPrint(f, r, c, columnZero, precision, width);
663 
664  int upper = print_diag ? 0 : 1;
665  for (int i = 0; i < r ; i++)
666  {
667  fprintf(f, "\n%*s ", columnZero, (const char *) data[i]->label);
668 
669  for (int j = 0; j < upper; j++)
670  fprintf(f, "%*.*s ", width[j], precision[j], " ");
671  for (int j = upper; j < c; j++)
672  fprintf(f, "%*.*f ", width[j], precision[j], (*this)[i][j]);
673 
674  upper++;
675  }
676 
677  fprintf(f, "\n");
678 
679  delete [] precision;
680  delete [] width;
681 }
682 
683 void Matrix::PrintLower(FILE * f, int r, int c, bool print_diag)
684 {
685  if (r == -1 || r > rows) r = rows;
686  if (c == -1 || c > cols) c = cols;
687 
688  String dimensions;
689  dimensions.printf("[%d x %d]", r, c);
690 
691  int columnZero = label.Length() > 15 ? label.Length() : 15;
692 
693  fprintf(f, "\n%*s =\n%*s ", columnZero, (const char *) label,
694  columnZero, (const char *) dimensions);
695 
696  int * precision = new int [c + 1];
697  int * width = new int [c + 1];
698 
699  for (int j = 0; j < c; j++)
700  {
701  precision[j] = extras[j].GetPrecision();
702  width[j] = extras[j].GetWidth();
703  fprintf(f, "%*s ", width[j], (const char *) extras[j].label);
704  }
705 
706  int upper = print_diag ? 1 : 0;
707 
708  for (int i = 0; i < r ; i++)
709  {
710  fprintf(f, "\n%*s ", columnZero, (const char *) data[i]->label);
711  for (int j = 0; j < upper; j++)
712  fprintf(f, "%*.*f ", width[j], precision[j],(*this)[i][j]);
713  for (int j = upper; j < c; j++)
714  fprintf(f, "%*.*s ", width[j], precision[j]," ");
715 
716  upper++;
717  }
718 
719  fprintf(f, "\n");
720 
721  delete [] precision;
722  delete [] width;
723 }
724 
725 
726 void Matrix::SetupPrint(FILE *f, int r, int c, int & column_zero, int * precision, int * width)
727 {
728  if (r == -1 || r > rows) r = rows;
729  if (c == -1 || c > cols) c = cols;
730 
731  String dimensions;
732  dimensions.printf("[%d x %d]", r, c);
733 
734  column_zero = label.Length() > 15 ? label.Length() : 15;
735 
736  fprintf(f, "\n%*s =\n%*s ", column_zero, (const char *) label,
737  column_zero, (const char *) dimensions);
738 
739  precision = new int [c + 1];
740  width = new int [c + 1];
741 
742  for (int j = 0; j < c; j++)
743  {
744  precision[j] = extras[j].GetPrecision();
745  width[j] = extras[j].GetWidth();
746  fprintf(f, "%*s ", width[j], (const char *) extras[j].label);
747  }
748 }