00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include "MathMatrix.h"
00019 #include "MathVector.h"
00020 #include "MathConstant.h"
00021 #include "Sort.h"
00022 #include "Error.h"
00023
00024 #include <string.h>
00025 #include <math.h>
00026 #include <stdio.h>
00027
00028 int Matrix::alloc = 2;
00029
00030 Matrix::~Matrix()
00031 {
00032
00033
00034 for (int i=0; i<size; i++)
00035 delete data[i];
00036
00037 if (size)
00038 delete [] data;
00039
00040 if (extraSize)
00041 delete [] extras;
00042 }
00043
00044 void Matrix::Init()
00045 {
00046 rows = cols = extraSize = size = 0;
00047 data = NULL;
00048 extras = NULL;
00049 label = "[Matrix]";
00050 }
00051
00052 void Matrix::SetLabel(const char * name)
00053 {
00054 label = '[';
00055 label += name;
00056 label += ']';
00057 }
00058
00059 void Matrix::Dimension(int m, int n)
00060 {
00061 if (n == cols && m == rows)
00062 return;
00063
00064 if (n > extraSize)
00065 {
00066 int newSize = (n + alloc) / alloc * alloc;
00067 ColumnExtras * newExtras = new ColumnExtras [newSize];
00068
00069 if (extras != NULL)
00070 for (int i = 0; i < extraSize; i++)
00071 newExtras[i] = extras[i];
00072
00073 if (extraSize)
00074 delete [] extras;
00075
00076 extraSize = newSize;
00077 extras = newExtras;
00078 }
00079
00080 if (m > size)
00081 {
00082 int newSize = (m + alloc) / alloc * alloc;
00083 Vector ** newData = new Vector * [newSize];
00084
00085 if (data != NULL)
00086 for (int i = 0; i < size; i++)
00087 newData[i] = data[i];
00088
00089 for (int i = size; i < newSize; i++)
00090 newData[i] = new Vector(n);
00091
00092 if (size)
00093 delete [] data;
00094
00095 size = newSize;
00096 data = newData;
00097 }
00098
00099 if (cols != n)
00100 for (int i = 0; i < rows; i++)
00101 data[i]->Dimension(n);
00102
00103 if (rows != m)
00104 for (int i = rows; i < m; i++)
00105 data[i]->Dimension(n);
00106
00107 rows = m;
00108 cols = n;
00109 }
00110
00111 void Matrix::Dimension(int m, int n, double value)
00112 {
00113 int originalRows = rows;
00114 int originalColumns = cols;
00115
00116 Dimension(m, n);
00117
00118 if (rows > originalRows)
00119 for (int i = originalRows; i < rows; i++)
00120 data[i]->Set(value);
00121
00122 if (cols > originalColumns)
00123 for (int i = 0; i < originalRows; i++)
00124 for (int j = originalColumns; j < cols; j++)
00125 data[i]->data[j] = value;
00126 }
00127
00128 void Matrix::Zero()
00129 {
00130 for (int i = 0; i < rows; i++)
00131 for (int j = 0; j < cols; j++)
00132 (*(data[i]))[j] = 0.0;
00133 }
00134
00135 void Matrix::Identity()
00136 {
00137 if (rows != cols)
00138 error("Matrix.Identity - Identity matrices must be square");
00139
00140 for (int i = 0; i < rows; i++)
00141 for (int j = 0; j < cols; j++)
00142 if (i == j)
00143 (*(data[i]))[j] = 1.0;
00144 else
00145 (*(data[i]))[j] = 0.0;
00146 }
00147
00148 void Matrix::Set(double k)
00149 {
00150 for (int i = 0; i < rows; i++)
00151 for (int j = 0; j < cols; j++)
00152 (*(data[i]))[j] = k;
00153 }
00154
00155 void Matrix::Negate()
00156 {
00157 for (int i = 0; i < rows; i++)
00158 for (int j = 0; j < cols; j++)
00159 (*(data[i]))[j] = -(*(data[i]))[j];
00160 }
00161
00162 void Matrix::Copy(const Matrix & m)
00163 {
00164 Dimension(m.rows, m.cols);
00165
00166 if (m.data != NULL)
00167 for (int i = 0; i < rows; i++)
00168 for (int j = 0; j < cols; j++)
00169 (*this)[i][j] = m[i][j];
00170 }
00171
00172 void Matrix::Transpose(const Matrix & m)
00173 {
00174 Dimension(m.cols, m.rows);
00175
00176 for (int i = 0; i < rows; i++)
00177 for (int j = 0; j < cols; j++)
00178 (*(data[i]))[j] = m[j][i];
00179 }
00180
00181 void Matrix::Add(double k)
00182 {
00183 for (int i = 0; i < rows; i++)
00184 for (int j = 0; j < cols; j++)
00185 (*(data[i]))[j] += k;
00186 }
00187
00188 void Matrix::Multiply(double k)
00189 {
00190 for (int i = 0; i < rows; i++)
00191 for (int j = 0; j < cols; j++)
00192 (*(data[i]))[j] *= k;
00193 }
00194
00195 void Matrix::Add(const Matrix & m)
00196 {
00197 if ((rows != m.rows) && (cols != m.cols))
00198 error("Matrix.Add - Attempted to add incompatible matrices\n"
00199 "Matrices - %s [%d, %d] + %s [%d, %d]\n",
00200 (const char *) label, rows, cols,
00201 (const char *) m.label, m.rows, m.cols);
00202
00203 for (int i = 0; i < rows; i++)
00204 for (int j = 0; j < cols; j++)
00205 (*(data[i]))[j] += m[i][j];
00206 }
00207
00208 void Matrix::AddMultiple(double k, const Matrix & m)
00209 {
00210 if ((rows != m.rows) && (cols != m.cols))
00211 error("Matrix.AddMultiple - Attempted to add incompatible matrices\n"
00212 "Matrices - %s [%d, %d] + k * %s [%d, %d]\n",
00213 (const char *) label, rows, cols,
00214 (const char *) m.label, m.rows, m.cols);
00215
00216 for (int i = 0; i < rows; i++)
00217 for (int j = 0; j < cols; j++)
00218 (*(data[i]))[j] += k * m[i][j];
00219 }
00220
00221
00222 void Matrix::Product(const Matrix & l, const Matrix & r)
00223 {
00224 if (l.cols != r.rows)
00225 error("Matrix.Multiply - Attempted to multiply incompatible matrices\n"
00226 "Matrices - %s [%d, %d] + %s [%d, %d]\n",
00227 (const char *) l.label, l.rows, l.cols,
00228 (const char *) r.label, r.rows, r.cols);
00229
00230 Dimension(l.rows, r.cols);
00231 Zero();
00232
00233 for (int k = 0; k < l.cols; k++)
00234 for (int i = 0; i < rows; i++)
00235 for (int j = 0; j < cols; j++)
00236 (*(data[i]))[j] += l[i][k] * r[k][j];
00237 }
00238
00239 void Matrix::AddRows(double k, int r1, int r2)
00240 {
00241 Vector v(*(data[r1]));
00242
00243 v.Multiply(k);
00244 data[r2]->Add(v);
00245 }
00246
00247 void Matrix::MultiplyRow(int r1, double k)
00248 {
00249 data[r1]->Multiply(k);
00250 }
00251
00252 void Matrix::AddRows(int r1, int r2)
00253 {
00254 data[r2]->Add(*(data[r1]));
00255 }
00256
00257 void Matrix::Reduce(double tol)
00258 {
00259 double pivot;
00260 int pivotr = 0;
00261 int r = 0;
00262
00263 for (int j = 0; j < cols; j++)
00264 {
00265 if (r > rows)
00266 return;
00267
00268 pivot = 0.0;
00269 for (int i = r; i < rows; i++)
00270 if (fabs((*this)[i][j]) > pivot)
00271 {
00272 pivot = fabs((*this)[i][j]);
00273 pivotr = i;
00274 }
00275
00276 if (pivot <= tol)
00277 {
00278 for (int i = r; i < rows; i++)
00279 (*this)[i][j] = 0.0;
00280 continue;
00281 }
00282
00283 SwapRows(pivotr, r);
00284
00285 double scale = (*this)[r][j];
00286
00287 (*this)[r][j] = 1.0;
00288 for (int k = j+1; k < cols; k++)
00289 (*this)[r][k] /= scale;
00290
00291 for (int i = r + 1; r < rows; i++)
00292 {
00293 scale = (*this)[i][j];
00294 (*this)[i][j] = 0.0;
00295 for (int k = j+1; k < cols; k++)
00296 (*this)[i][k] -= (*this)[r][k] * scale;
00297 }
00298
00299 r++;
00300 }
00301 }
00302
00303 void Matrix::DeleteRow(int r)
00304 {
00305 Vector * temp = data[r];
00306
00307 for (int i = r + 1; i < rows; i++)
00308 data[i-1] = data[i];
00309
00310 data[rows - 1] = temp;
00311 rows--;
00312 }
00313
00314 void Matrix::DeleteColumn(int c)
00315 {
00316 for (int i = 0; i < rows; i++)
00317 data[i] -> DeleteDimension(c);
00318
00319 for (int i = c + 1; i < cols; i++)
00320 extras[i-1] = extras[i];
00321
00322 cols--;
00323 }
00324
00325 void Matrix::SwapColumns(int c1, int c2)
00326 {
00327 double temp;
00328
00329 for (int i = 0; i < rows; i++)
00330 {
00331 temp = (*data[i])[c1];
00332 (*data[i])[c1] = (*data[i])[c2];
00333 (*data[i])[c2] = temp;
00334 }
00335
00336 extras[c1].Swap(extras[c2]);
00337 }
00338
00339 void Matrix::Read(FILE * f)
00340 {
00341 int r, c;
00342 char buffer[100];
00343 int numItems = 0;
00344
00345 numItems = fscanf(f, " %s =", buffer);
00346 if(numItems != 1) { }
00347 buffer[strlen(buffer) - 1] = 0;
00348 SetLabel(buffer);
00349
00350 numItems = fscanf(f, " [ %d x %d ]", &r, &c);
00351 if(numItems != 2) { }
00352 Dimension(r, c);
00353
00354 for (int c = 0; c < cols; c++)
00355 {
00356 numItems = fscanf(f, " %s", buffer);
00357 if(numItems != 1) { }
00358 SetColumnLabel(c, buffer);
00359 }
00360
00361 for (int r = 0; r < rows; r++)
00362 for (int c = 0; c < cols; c++)
00363 {
00364 numItems = fscanf(f, " %lf", &((*this)[r][c]));
00365 if(numItems != 1) { }
00366 }
00367 }
00368
00369
00370 void Matrix::Print(FILE * f, int r, int c)
00371 {
00372 if (r == -1 || r > rows) r = rows;
00373 if (c == -1 || c > cols) c = cols;
00374
00375 char dimensions[30];
00376
00377 sprintf(dimensions, "[%d x %d]", r, c);
00378
00379 int columnZero = label.Length() > 15 ? label.Length() : 15;
00380
00381 fprintf(f, "\n%*s =\n%*s ", columnZero, (const char *) label,
00382 columnZero, dimensions);
00383
00384 int * precision = new int [c + 1];
00385 int * width = new int [c + 1];
00386
00387 for (int j = 0; j < c; j++)
00388 {
00389 precision[j] = extras[j].GetPrecision();
00390 width[j] = extras[j].GetWidth();
00391 fprintf(f, "%*s ", width[j], (const char *) extras[j].label);
00392 }
00393
00394 for (int i = 0; i < r; i++)
00395 {
00396 fprintf(f, "\n%*s ", columnZero, (const char *) data[i]->label);
00397 for (int j = 0; j < c; j++)
00398 fprintf(f, "%*.*f ", width[j], precision[j], (*this)[i][j]);
00399 }
00400
00401 fprintf(f, "\n");
00402
00403 delete [] precision;
00404 delete [] width;
00405 }
00406
00407 void Matrix::CopyLabels(Matrix & M)
00408 {
00409 for (int i = 0; i < rows; i++)
00410 if (i < M.rows)
00411 data[i]->SetLabel(M[i].label);
00412
00413 for (int i = 0; i < cols; i++)
00414 if (i < M.cols)
00415 SetColumnLabel(i, M.GetColumnLabel(i));
00416 }
00417
00418
00419
00420
00421 void ColumnExtras::Init()
00422 {
00423 label = "column";
00424 dirty = true;
00425 precision = 3;
00426 width = 7;
00427 }
00428
00429 ColumnExtras::~ColumnExtras()
00430 { }
00431
00432 void ColumnExtras::SetLabel(const char * name)
00433 {
00434 label = name;
00435 }
00436
00437 int ColumnExtras::GetWidth()
00438 {
00439 if (dirty)
00440 {
00441 if (precision + 2 > width)
00442 width = precision + 2;
00443 if (label.Length() > width)
00444 width = label.Length();
00445 dirty = false;
00446 }
00447 return width;
00448 }
00449
00450 void ColumnExtras::Copy(ColumnExtras & c)
00451 {
00452 width = c.width;
00453 precision = c.precision;
00454 dirty = c.dirty;
00455 label = c.label;
00456 }
00457
00458 #define SWAP(a,b) {int swap=(a); (a)=(b); (b)=swap;}
00459 #define SWAPBOOL(a,b) {bool swap=(a); (a)=(b); (b)=swap;}
00460
00461 void ColumnExtras::Swap(ColumnExtras & c)
00462 {
00463 SWAP(c.width, width);
00464 SWAP(c.precision, precision);
00465 SWAPBOOL(c.dirty, dirty);
00466 c.label.Swap(label);
00467 }
00468
00469 int Matrix::CompareRows(Vector ** row1, Vector ** row2)
00470 {
00471 if ((**row1)[0] < (**row2)[0]) return -1;
00472 if ((**row1)[0] > (**row2)[0]) return 1;
00473 return 0;
00474 }
00475
00476 void Matrix::Sort()
00477 {
00478 QuickSort(data, rows, sizeof(Vector *), COMPAREFUNC CompareRows);
00479 }
00480
00481 bool Matrix::operator == (const Matrix & rhs) const
00482 {
00483 if (rhs.rows != rows || rhs.cols != cols) return false;
00484
00485 for (int i = 0; i < rows; i++)
00486 if ((*this)[i] != rhs[i])
00487 return false;
00488 return true;
00489 }
00490
00491 void Matrix::StackBottom(const Matrix & m)
00492 {
00493 if (m.cols != cols)
00494 error("Attempted to stack matrices with different number of columns");
00495
00496 int end = rows;
00497
00498 Dimension(rows + m.rows, cols);
00499
00500 for (int i = 0; i < m.rows; i++)
00501 *(data[i + end]) = m[i];
00502 }
00503
00504 void Matrix::StackLeft(const Matrix & m)
00505 {
00506 if (m.rows != rows)
00507 error("Attempted to stack matrics with different numbers of rows");
00508
00509 for (int i = 0; i < rows; i++)
00510 data[i]->Stack(m[i]);
00511
00512 Dimension(rows, cols + m.cols);
00513 }
00514
00515 void Matrix::Swap(Matrix & m)
00516 {
00517 label.Swap(m.label);
00518
00519 ColumnExtras * tmpExtras = extras;
00520 extras = m.extras;
00521 m.extras = tmpExtras;
00522
00523 int swap;
00524 swap = rows;
00525 rows = m.rows;
00526 m.rows = swap;
00527 swap = cols;
00528 cols = m.cols;
00529 m.cols = swap;
00530 swap = size;
00531 size = m.size;
00532 m.size = swap;
00533 swap = extraSize;
00534 extraSize = m.extraSize;
00535 m.extraSize = swap;
00536
00537 Vector ** tmpData = data;
00538 data = m.data;
00539 m.data = tmpData;
00540 }
00541
00542 double Matrix::Min() const
00543 {
00544 if (rows == 0 || cols == 0)
00545 return 0.0;
00546
00547 double minimum = data[0]->Min();
00548
00549 for (int i = 1; i < rows; i++)
00550 minimum = min(data[i]->Min(), minimum);
00551
00552 return minimum;
00553 }
00554
00555 double Matrix::Max() const
00556 {
00557 if (rows == 0 || cols == 0)
00558 return 0.0;
00559
00560 double maximum = data[0]->Max();
00561
00562 for (int i = 1; i < rows; i++)
00563 maximum = max(data[i]->Max(), maximum);
00564
00565 return maximum;
00566 }
00567
00568 double Matrix::Mean() const
00569 {
00570 if (rows == 0 || cols == 0)
00571 return 0.0;
00572
00573 double sum = data[0]->Sum();
00574
00575 for (int i = 1; i < rows; i++)
00576 sum += data[i]->Sum();
00577
00578 return sum / (rows * cols);
00579 }
00580
00581 double Matrix::SafeMin() const
00582 {
00583 double lo = (rows > 0 && cols > 0) ? _NAN_ : 0.0;
00584
00585 int i, j;
00586
00587 for (i = 0; i < rows; i++)
00588 {
00589 for (j = 0; j < cols; j++)
00590 if (data[i]->data[j] != _NAN_)
00591 {
00592 lo = data[i]->data[j];
00593 break;
00594 }
00595 if (j != cols) break;
00596 }
00597
00598 for (; i < rows; i++, j = 0)
00599 for (; j < cols; j++)
00600 if (data[i]->data[j] < lo && data[i]->data[j] != _NAN_)
00601 lo = data[i]->data[j];
00602
00603 return lo;
00604 }
00605
00606 double Matrix::SafeMax() const
00607 {
00608 double hi = (rows > 0 && cols > 0) ? _NAN_ : 0.0;
00609
00610 int i, j;
00611
00612 for (i = 0; i < rows; i++)
00613 {
00614 for (j = 0; j < cols; j++)
00615 if (data[i]->data[j] != _NAN_)
00616 {
00617 hi = data[i]->data[j];
00618 break;
00619 }
00620 if (j != cols) break;
00621 }
00622
00623 for (; i < rows; i++, j = 0)
00624 for (; j < cols; j++)
00625 if (data[i]->data[j] > hi && data[i]->data[j] != _NAN_)
00626 hi = data[i]->data[j];
00627
00628 return hi;
00629 }
00630
00631 double Matrix::SafeMean() const
00632 {
00633 double sum = 0.0;
00634 int count = 0;
00635
00636 for (int i = 0; i < rows; i++)
00637 for (int j = 0; j < cols; j++)
00638 if ((*this)[i][j] != _NAN_)
00639 {
00640 sum += (*this)[i][j];
00641 count ++;
00642 }
00643
00644 return (count) ? sum / count : 0.0;
00645 }
00646
00647 int Matrix::SafeCount() const
00648 {
00649 int total = 0;
00650
00651 for (int i = 0; i < rows; i++)
00652 total += data[i]->SafeCount();
00653
00654 return total;
00655 }
00656
00657 void Matrix::PrintUpper(FILE * f, int r, int c, bool print_diag)
00658 {
00659 int columnZero;
00660 int * precision = NULL, * width = NULL;
00661
00662 SetupPrint(f, r, c, columnZero, precision, width);
00663
00664 int upper = print_diag ? 0 : 1;
00665 for (int i = 0; i < r ; i++)
00666 {
00667 fprintf(f, "\n%*s ", columnZero, (const char *) data[i]->label);
00668
00669 for (int j = 0; j < upper; j++)
00670 fprintf(f, "%*.*s ", width[j], precision[j], " ");
00671 for (int j = upper; j < c; j++)
00672 fprintf(f, "%*.*f ", width[j], precision[j], (*this)[i][j]);
00673
00674 upper++;
00675 }
00676
00677 fprintf(f, "\n");
00678
00679 delete [] precision;
00680 delete [] width;
00681 }
00682
00683 void Matrix::PrintLower(FILE * f, int r, int c, bool print_diag)
00684 {
00685 if (r == -1 || r > rows) r = rows;
00686 if (c == -1 || c > cols) c = cols;
00687
00688 String dimensions;
00689 dimensions.printf("[%d x %d]", r, c);
00690
00691 int columnZero = label.Length() > 15 ? label.Length() : 15;
00692
00693 fprintf(f, "\n%*s =\n%*s ", columnZero, (const char *) label,
00694 columnZero, (const char *) dimensions);
00695
00696 int * precision = new int [c + 1];
00697 int * width = new int [c + 1];
00698
00699 for (int j = 0; j < c; j++)
00700 {
00701 precision[j] = extras[j].GetPrecision();
00702 width[j] = extras[j].GetWidth();
00703 fprintf(f, "%*s ", width[j], (const char *) extras[j].label);
00704 }
00705
00706 int upper = print_diag ? 1 : 0;
00707
00708 for (int i = 0; i < r ; i++)
00709 {
00710 fprintf(f, "\n%*s ", columnZero, (const char *) data[i]->label);
00711 for (int j = 0; j < upper; j++)
00712 fprintf(f, "%*.*f ", width[j], precision[j],(*this)[i][j]);
00713 for (int j = upper; j < c; j++)
00714 fprintf(f, "%*.*s ", width[j], precision[j]," ");
00715
00716 upper++;
00717 }
00718
00719 fprintf(f, "\n");
00720
00721 delete [] precision;
00722 delete [] width;
00723 }
00724
00725
00726 void Matrix::SetupPrint(FILE *f, int r, int c, int & column_zero, int * precision, int * width)
00727 {
00728 if (r == -1 || r > rows) r = rows;
00729 if (c == -1 || c > cols) c = cols;
00730
00731 String dimensions;
00732 dimensions.printf("[%d x %d]", r, c);
00733
00734 column_zero = label.Length() > 15 ? label.Length() : 15;
00735
00736 fprintf(f, "\n%*s =\n%*s ", column_zero, (const char *) label,
00737 column_zero, (const char *) dimensions);
00738
00739 precision = new int [c + 1];
00740 width = new int [c + 1];
00741
00742 for (int j = 0; j < c; j++)
00743 {
00744 precision[j] = extras[j].GetPrecision();
00745 width[j] = extras[j].GetWidth();
00746 fprintf(f, "%*s ", width[j], (const char *) extras[j].label);
00747 }
00748 }