18#include "MathMatrix.h"
19#include "MathVector.h"
20#include "MathConstant.h"
34 for (
int i=0; i<size; i++)
46 rows = cols = extraSize = size = 0;
52void Matrix::SetLabel(
const char * name)
59void Matrix::Dimension(
int m,
int n)
61 if (n == cols && m == rows)
66 int newSize = (n + alloc) / alloc * alloc;
70 for (
int i = 0; i < extraSize; i++)
71 newExtras[i] = extras[i];
82 int newSize = (m + alloc) / alloc * alloc;
86 for (
int i = 0; i < size; i++)
89 for (
int i = size; i < newSize; i++)
90 newData[i] =
new Vector(n);
100 for (
int i = 0; i < rows; i++)
101 data[i]->Dimension(n);
104 for (
int i = rows; i < m; i++)
105 data[i]->Dimension(n);
111void Matrix::Dimension(
int m,
int n,
double value)
113 int originalRows = rows;
114 int originalColumns = cols;
118 if (rows > originalRows)
119 for (
int i = originalRows; i < rows; i++)
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;
130 for (
int i = 0; i < rows; i++)
131 for (
int j = 0; j < cols; j++)
132 (*(data[i]))[j] = 0.0;
135void Matrix::Identity()
138 error(
"Matrix.Identity - Identity matrices must be square");
140 for (
int i = 0; i < rows; i++)
141 for (
int j = 0; j < cols; j++)
143 (*(data[i]))[j] = 1.0;
145 (*(data[i]))[j] = 0.0;
148void Matrix::Set(
double k)
150 for (
int i = 0; i < rows; i++)
151 for (
int j = 0; j < cols; j++)
157 for (
int i = 0; i < rows; i++)
158 for (
int j = 0; j < cols; j++)
159 (*(data[i]))[j] = -(*(data[i]))[j];
162void Matrix::Copy(
const Matrix & m)
164 Dimension(m.rows, m.cols);
167 for (
int i = 0; i < rows; i++)
168 for (
int j = 0; j < cols; j++)
169 (*
this)[i][j] = m[i][j];
172void Matrix::Transpose(
const Matrix & m)
174 Dimension(m.cols, m.rows);
176 for (
int i = 0; i < rows; i++)
177 for (
int j = 0; j < cols; j++)
178 (*(data[i]))[j] = m[j][i];
181void Matrix::Add(
double k)
183 for (
int i = 0; i < rows; i++)
184 for (
int j = 0; j < cols; j++)
185 (*(data[i]))[j] += k;
188void Matrix::Multiply(
double k)
190 for (
int i = 0; i < rows; i++)
191 for (
int j = 0; j < cols; j++)
192 (*(data[i]))[j] *= k;
195void Matrix::Add(
const Matrix & m)
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);
203 for (
int i = 0; i < rows; i++)
204 for (
int j = 0; j < cols; j++)
205 (*(data[i]))[j] += m[i][j];
208void Matrix::AddMultiple(
double k,
const Matrix & m)
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);
216 for (
int i = 0; i < rows; i++)
217 for (
int j = 0; j < cols; j++)
218 (*(data[i]))[j] += k * m[i][j];
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);
230 Dimension(l.rows, r.cols);
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];
239void Matrix::AddRows(
double k,
int r1,
int r2)
247void Matrix::MultiplyRow(
int r1,
double k)
249 data[r1]->Multiply(k);
252void Matrix::AddRows(
int r1,
int r2)
254 data[r2]->Add(*(data[r1]));
257void Matrix::Reduce(
double tol)
263 for (
int j = 0; j < cols; j++)
269 for (
int i = r; i < rows; i++)
270 if (fabs((*
this)[i][j]) > pivot)
272 pivot = fabs((*
this)[i][j]);
278 for (
int i = r; i < rows; i++)
285 double scale = (*this)[r][j];
288 for (
int k = j+1; k < cols; k++)
289 (*
this)[r][k] /= scale;
291 for (
int i = r + 1; r < rows; i++)
293 scale = (*this)[i][j];
295 for (
int k = j+1; k < cols; k++)
296 (*
this)[i][k] -= (*this)[r][k] * scale;
303void Matrix::DeleteRow(
int r)
307 for (
int i = r + 1; i < rows; i++)
310 data[rows - 1] = temp;
314void Matrix::DeleteColumn(
int c)
316 for (
int i = 0; i < rows; i++)
317 data[i] -> DeleteDimension(c);
319 for (
int i = c + 1; i < cols; i++)
320 extras[i-1] = extras[i];
325void Matrix::SwapColumns(
int c1,
int c2)
329 for (
int i = 0; i < rows; i++)
331 temp = (*data[i])[c1];
332 (*data[i])[c1] = (*data[i])[c2];
333 (*data[i])[c2] = temp;
336 extras[c1].Swap(extras[c2]);
339void Matrix::Read(FILE * f)
345 numItems = fscanf(f,
" %s =", buffer);
346 if(numItems != 1) { }
347 buffer[strlen(buffer) - 1] = 0;
350 numItems = fscanf(f,
" [ %d x %d ]", &r, &c);
351 if(numItems != 2) { }
354 for (
int c = 0; c < cols; c++)
356 numItems = fscanf(f,
" %s", buffer);
357 if(numItems != 1) { }
358 SetColumnLabel(c, buffer);
361 for (
int r = 0; r < rows; r++)
362 for (
int c = 0; c < cols; c++)
364 numItems = fscanf(f,
" %lf", &((*
this)[r][c]));
365 if(numItems != 1) { }
370void Matrix::Print(FILE * f,
int r,
int c)
372 if (r == -1 || r > rows) r = rows;
373 if (c == -1 || c > cols) c = cols;
377 sprintf(dimensions,
"[%d x %d]", r, c);
379 int columnZero = label.Length() > 15 ? label.Length() : 15;
381 fprintf(f,
"\n%*s =\n%*s ", columnZero, (
const char *) label,
382 columnZero, dimensions);
384 int * precision =
new int [c + 1];
385 int * width =
new int [c + 1];
387 for (
int j = 0; j < c; j++)
389 precision[j] = extras[j].GetPrecision();
390 width[j] = extras[j].GetWidth();
391 fprintf(f,
"%*s ", width[j], (
const char *) extras[j].label);
394 for (
int i = 0; i < r; i++)
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]);
407void Matrix::CopyLabels(
Matrix & M)
409 for (
int i = 0; i < rows; i++)
411 data[i]->SetLabel(M[i].label);
413 for (
int i = 0; i < cols; i++)
415 SetColumnLabel(i, M.GetColumnLabel(i));
421void ColumnExtras::Init()
429ColumnExtras::~ColumnExtras()
432void ColumnExtras::SetLabel(
const char * name)
437int ColumnExtras::GetWidth()
441 if (precision + 2 > width)
442 width = precision + 2;
443 if (label.Length() > width)
444 width = label.Length();
453 precision = c.precision;
458#define SWAP(a,b) {int swap=(a); (a)=(b); (b)=swap;}
459#define SWAPBOOL(a,b) {bool swap=(a); (a)=(b); (b)=swap;}
463 SWAP(c.width, width);
464 SWAP(c.precision, precision);
465 SWAPBOOL(c.dirty, dirty);
471 if ((**row1)[0] < (**row2)[0])
return -1;
472 if ((**row1)[0] > (**row2)[0])
return 1;
478 QuickSort(data, rows,
sizeof(
Vector *), COMPAREFUNC CompareRows);
481bool Matrix::operator == (
const Matrix & rhs)
const
483 if (rhs.rows != rows || rhs.cols != cols)
return false;
485 for (
int i = 0; i < rows; i++)
486 if ((*
this)[i] != rhs[i])
491void Matrix::StackBottom(
const Matrix & m)
494 error(
"Attempted to stack matrices with different number of columns");
498 Dimension(rows + m.rows, cols);
500 for (
int i = 0; i < m.rows; i++)
501 *(data[i + end]) = m[i];
504void Matrix::StackLeft(
const Matrix & m)
507 error(
"Attempted to stack matrics with different numbers of rows");
509 for (
int i = 0; i < rows; i++)
510 data[i]->Stack(m[i]);
512 Dimension(rows, cols + m.cols);
515void Matrix::Swap(
Matrix & m)
521 m.extras = tmpExtras;
534 extraSize = m.extraSize;
542double Matrix::Min()
const
544 if (rows == 0 || cols == 0)
547 double minimum = data[0]->Min();
549 for (
int i = 1; i < rows; i++)
550 minimum = min(data[i]->Min(), minimum);
555double Matrix::Max()
const
557 if (rows == 0 || cols == 0)
560 double maximum = data[0]->Max();
562 for (
int i = 1; i < rows; i++)
563 maximum = max(data[i]->Max(), maximum);
568double Matrix::Mean()
const
570 if (rows == 0 || cols == 0)
573 double sum = data[0]->Sum();
575 for (
int i = 1; i < rows; i++)
576 sum += data[i]->Sum();
578 return sum / (rows * cols);
581double Matrix::SafeMin()
const
583 double lo = (rows > 0 && cols > 0) ? _NAN_ : 0.0;
587 for (i = 0; i < rows; i++)
589 for (j = 0; j < cols; j++)
590 if (data[i]->data[j] != _NAN_)
592 lo = data[i]->data[j];
595 if (j != cols)
break;
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];
606double Matrix::SafeMax()
const
608 double hi = (rows > 0 && cols > 0) ? _NAN_ : 0.0;
612 for (i = 0; i < rows; i++)
614 for (j = 0; j < cols; j++)
615 if (data[i]->data[j] != _NAN_)
617 hi = data[i]->data[j];
620 if (j != cols)
break;
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];
631double Matrix::SafeMean()
const
636 for (
int i = 0; i < rows; i++)
637 for (
int j = 0; j < cols; j++)
638 if ((*
this)[i][j] != _NAN_)
640 sum += (*this)[i][j];
644 return (count) ? sum / count : 0.0;
647int Matrix::SafeCount()
const
651 for (
int i = 0; i < rows; i++)
652 total += data[i]->SafeCount();
657void Matrix::PrintUpper(FILE * f,
int r,
int c,
bool print_diag)
660 int * precision = NULL, * width = NULL;
662 SetupPrint(f, r, c, columnZero, precision, width);
664 int upper = print_diag ? 0 : 1;
665 for (
int i = 0; i < r ; i++)
667 fprintf(f,
"\n%*s ", columnZero, (
const char *) data[i]->label);
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]);
683void Matrix::PrintLower(FILE * f,
int r,
int c,
bool print_diag)
685 if (r == -1 || r > rows) r = rows;
686 if (c == -1 || c > cols) c = cols;
689 dimensions.printf(
"[%d x %d]", r, c);
691 int columnZero = label.Length() > 15 ? label.Length() : 15;
693 fprintf(f,
"\n%*s =\n%*s ", columnZero, (
const char *) label,
694 columnZero, (
const char *) dimensions);
696 int * precision =
new int [c + 1];
697 int * width =
new int [c + 1];
699 for (
int j = 0; j < c; j++)
701 precision[j] = extras[j].GetPrecision();
702 width[j] = extras[j].GetWidth();
703 fprintf(f,
"%*s ", width[j], (
const char *) extras[j].label);
706 int upper = print_diag ? 1 : 0;
708 for (
int i = 0; i < r ; i++)
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],
" ");
726void Matrix::SetupPrint(FILE *f,
int r,
int c,
int & column_zero,
int * precision,
int * width)
728 if (r == -1 || r > rows) r = rows;
729 if (c == -1 || c > cols) c = cols;
732 dimensions.printf(
"[%d x %d]", r, c);
734 column_zero = label.Length() > 15 ? label.Length() : 15;
736 fprintf(f,
"\n%*s =\n%*s ", column_zero, (
const char *) label,
737 column_zero, (
const char *) dimensions);
739 precision =
new int [c + 1];
740 width =
new int [c + 1];
742 for (
int j = 0; j < c; j++)
744 precision[j] = extras[j].GetPrecision();
745 width[j] = extras[j].GetWidth();
746 fprintf(f,
"%*s ", width[j], (
const char *) extras[j].label);