00001 #ifndef CxMatrixTem_h
00002 #define CxMatrixTem_h
00003
00004 #include "Data/Basis/CxStream.h"
00005 #include "Data/Element/CxArithTypes.h"
00006 #include "CxVector.h"
00007
00008
00009
00010
00011
00015 template<class C>
00016 class CxMatrixTem {
00017 public:
00018
00021
00023 CxMatrixTem()
00024 {
00025 _nr = 0;
00026 _nc = 0;
00027 _data = 0;
00028 }
00029
00031 CxMatrixTem(int nRow, int nCol)
00032 {
00033 _nr = nRow;
00034 _nc = nCol;
00035 _data = new C[_nr * _nc];
00036 }
00037
00039 CxMatrixTem(int nRow, int nCol, C a)
00040 {
00041 _nr = nRow;
00042 _nc = nCol;
00043 _data = new C[_nr * _nc];
00044
00045 C* t = _data;
00046 int i = nElem();
00047 while (--i >= 0)
00048 *t++ = a;
00049 }
00050
00052 CxMatrixTem(int nRow, int nCol, C* data)
00053 {
00054 _nr = nRow;
00055 _nc = nCol;
00056 _data = data;
00057 }
00058
00060 CxMatrixTem(const CxMatrixTem &m)
00061 {
00062 _nr = m.nRow();
00063 _nc = m.nCol();
00064 _data = new C[_nr * _nc];
00065
00066 C* t = m._data;
00067 C* u = _data;
00068 int i = m.nElem();
00069 while (--i >= 0)
00070 *u++ = *t++;
00071 }
00072
00074
00075
00077
00078
00079 ~CxMatrixTem()
00080 {
00081
00082 }
00083
00086
00088 int
00089 nRow() const
00090 {
00091 return _nr;
00092 }
00093
00095 int
00096 nCol() const
00097 {
00098 return _nc;
00099 }
00100
00102 int
00103 nElem() const
00104 {
00105 return _nr * _nc;
00106 }
00107
00109 int
00110 valid() const
00111 {
00112 return ((_nc != 0) && (_nr != 0));
00113 }
00114
00115
00116 C*
00117 dataPtr() const
00118 {
00119 return _data;
00120 }
00121
00123
00126
00127
00129 CxMatrixTem<C>&
00130 operator=(const CxMatrixTem<C>& m)
00131 {
00132 if (this != &m) {
00133 delete [] _data;
00134 _nr = m.nRow();
00135 _nc = m.nCol();
00136 _data = new C [_nr * _nc];
00137 C *t = _data;
00138 C *u = m._data;
00139 int i = m.nElem();
00140 while (--i >= 0)
00141 *t++ = *u++;
00142 }
00143 return *this;
00144 }
00145
00147 C*
00148 operator[](int i) const
00149 {
00150 return &_data[i*_nc];
00151 }
00152
00154 CxMatrixTem<C>
00155 operator*(double val) const
00156 {
00157 CxMatrixTem<C> m(*this);
00158 C* t = m._data;
00159 int i = nElem();
00160 while (--i >= 0)
00161 *t++ = *t * val ;
00162 return m;
00163 }
00164
00166 CxMatrixTem<C>
00167 operator*(const CxMatrixTem<C>& a) const
00168 {
00169 if (nCol() != a.nRow()) {
00170 CX_CERR << "nonconformant CxMatrixTem * CxMatrixTem operands."
00171 << CX_ENDL;
00172 return CxMatrixTem<C>(0,0);
00173 }
00174 CxMatrixTem<C> m((*this).nRow(), a.nCol(), 0.0);
00175 double sum;
00176 int i, j, k;
00177 for (i=0 ; i<nRow() ; i++) {
00178 for (j=0 ; j<a.nCol() ; j++) {
00179 sum = 0;
00180 for (k=0 ; k<nCol() ; k++)
00181 sum += (*this)[i][k] * a[k][j];
00182 m[i][j] = sum;
00183 }
00184 }
00185 return m;
00186 }
00187
00189 template<class CC>
00190 friend CxVector<CC>
00191 operator*(const CxVector<CC>& , const CxMatrixTem<CC>& );
00192
00194 template<class CC>
00195 friend CxVector<CC>
00196 operator*(const CxMatrixTem<CC>& , const CxVector<CC>& );
00197
00199 CxMatrixTem<C>
00200 operator+(const CxMatrixTem<C>& a) const
00201 {
00202 if ((a.nCol() != nCol()) || (a.nRow() != nRow()) ) {
00203 CX_CERR << "nonconformant CxMatrixTem + CxMatrixTem operands." << CX_ENDL;
00204 return CxMatrixTem<C>(0,0);
00205 }
00206 CxMatrixTem<C> m(nRow(), nCol());
00207 int i, j;
00208 for (i=0 ; i<nRow() ; i++) {
00209 for (j=0 ; j<nCol() ; j++)
00210 m[i][j] = a[i][j] + (*this)[i][j];
00211 }
00212 return m;
00213 }
00214
00216 CxMatrixTem<C> operator+(double val) const
00217 {
00218
00219 CxMatrixTem<C> m(*this);
00220 C* t = m._data;
00221 int i = nElem();
00222 while (--i >= 0)
00223 *t++ = val + *t;
00224 return m;
00225 }
00226
00228
00231
00232 CxVector<C>
00233 getRow(long r) const
00234 {
00235
00236 if(r>= _nr || r<0){
00237 CX_CERR << "getRow: rownr exceeded the number of rows in CxMatrixTem."
00238 << CX_ENDL;
00239 return CxVector<C>(0);
00240 }
00241
00242 CxVector<C> tmp(_nc);
00243
00244 for(int i = 0; i<_nc; i++)
00245 tmp[i] = (*this)[r][i] ;
00246
00247 return tmp ;
00248 }
00249
00250 CxVector<C>
00251 getCol(long c) const
00252 {
00253
00254 if(c>= _nc || c<0){
00255 CX_CERR << "getCol: Colnr exceeded number of columns in CxMatrixTem."
00256 << CX_ENDL;
00257 return CxVector<C>(0);
00258 }
00259
00260 CxVector<C> tmp(_nr);
00261 for(int i = 0; i<_nr; i++)
00262 {
00263 tmp[i] = (*this)[i][c] ;
00264 }
00265 return tmp;
00266
00267 }
00268
00269
00270 CxMatrixTem<C>
00271 setRow(long r, CxVector<C> v, bool makeCopy = false)
00272 {
00273
00274 if(r>= _nr || r<0){
00275 CX_CERR << "setRow: rownr exceeded the number of rows in CxMatrixTem."
00276 << CX_ENDL;
00277 return CxMatrixTem<C>(0,0) ;
00278 }
00279 if(_nc != v.Size()){
00280 CX_CERR << "setRow: v.nElem not equal to matrix dimension" << CX_ENDL;
00281 return CxMatrixTem<C>(0,0) ;
00282 }
00283
00284 CxMatrixTem<C> tmp(0,0);
00285
00286 if(makeCopy) {
00287
00288 tmp = CxMatrixTem<C>(_nr,_nc);
00289
00290 for(int row = 0; row<_nr; row++)
00291 for(int col = 0; col<_nc; col++)
00292 if(row!=r) tmp[row][col] = (*this)[row][col] ;
00293 else tmp[row][col] = v[col] ;
00294 }
00295 else {
00296
00297 for(int col = 0; col<_nc; col++)
00298 (*this)[r][col] = v[col] ;
00299 }
00300
00301
00302 return tmp;
00303
00304 }
00305
00306 CxMatrixTem<C>
00307 setCol(long c, CxVector<C> v, bool makeCopy = false)
00308 {
00309
00310 if(c>= _nc || c<0){
00311 CX_CERR << "setCol: Colnr exceeded number of columns in CxMatrixTem."
00312 << CX_ENDL;
00313 return CxMatrixTem<C>(0,0) ;
00314 }
00315 if(_nr != v.Size()){
00316 CX_CERR << "setCol: v.nElem not equal to matrix dimension" << CX_ENDL;
00317 return CxMatrixTem<C>(0,0) ;
00318 }
00319
00320 CxMatrixTem<C> tmp(0,0);
00321
00322 if(makeCopy) {
00323
00324 tmp = CxMatrixTem<C>(_nr,_nc);
00325
00326 for(int row = 0; row<_nr; row++)
00327 for(int col = 0; col<_nc; col++)
00328 if(col!=c) tmp[row][col] = (*this)[row][col] ;
00329 else tmp[row][col] = v[row] ;
00330 }
00331 else {
00332
00333 for(int row = 0; row<_nr; row++)
00334 (*this)[row][c] = v[row] ;
00335 }
00336
00337 return tmp;
00338
00339 }
00340
00342 CxMatrixTem<C> i() const;
00343
00344 CxMatrixTem<C>
00345 t() const
00346 {
00347 CxMatrixTem<C> m(nCol(), nRow());
00348 int i, j;
00349 for (i=0 ; i<nRow() ; i++)
00350 for (j=0 ; j<nCol() ; j++)
00351 m[j][i] = (*this)[i][j];
00352 return m;
00353 }
00354
00355
00357 CxMatrixTem<C>
00358 covariance(CxVector<C> &mean) const ;
00359
00363 CxMatrixTem<C>
00364 svd(CxVector<C>& W, CxMatrixTem<C>& V) const;
00365
00369 CxMatrixTem<C>
00370 klm(int newDim) const;
00371
00373
00374
00375 friend class CxVector<C>;
00376
00377 int _nr;
00378 int _nc;
00379 C* _data;
00380
00381 };
00382
00383
00384 typedef CxMatrixTem<CxReal64> CxMatrix ;
00385
00386 template<class C>
00387 inline CxVector<C>
00388 operator*(const CxVector<C> &a, const CxMatrixTem<C> &b)
00389 {
00390 if (a.Size() != b.nRow()) {
00391 CX_CERR << "nonconformant CxVector<C> * CxMatrixTem operands." << CX_ENDL;
00392 return CxVector<C>(0);
00393 }
00394 CxVector<C> v(b.nCol());
00395 double sum;
00396 int i, j;
00397 for (i=0 ; i<b.nCol() ; i++) {
00398 sum = 0;
00399 for (j=0 ; j<b.nRow() ; j++)
00400 sum += a[j] * b[j][i];
00401 v[i] = sum;
00402 }
00403 return v;
00404 }
00405
00406
00407 template<class C>
00408 inline CxVector<C>
00409 operator*(const CxMatrixTem<C> &a, const CxVector<C> &b)
00410 {
00411 if (b.Size() != a.nCol()) {
00412 CX_CERR << "nonconformant CxMatrixTem * CxVector<C> operands." << CX_ENDL;
00413 return CxVector<C>(0);
00414 }
00415 CxVector<C> v(a.nRow());
00416 double sum;
00417 int i, j;
00418 for (i=0 ; i<a.nRow() ; i++) {
00419 sum = 0;
00420 for (j=0 ; j<a.nCol() ; j++)
00421 sum += a[i][j] * b[j];
00422 v[i] = sum;
00423 }
00424 return v;
00425 }
00426
00427
00428 template<class C>
00429 static void svdcmp(C* a, CxVector<C> &w, C* v, int n, int m);
00430
00431 template<class C>
00432 static int ludcmp(C* a, int n, short *indx, double *d);
00433
00434 template<class C>
00435 static void lubksb(C* a, int n, short *indx, double *b);
00436
00437 template<class C> CxMatrixTem<C>
00438 CxMatrixTem<C>::i() const
00439 {
00440 if (nRow() != nCol()) {
00441 CX_CERR << "Inverse: matrix is not square!." << CX_ENDL;
00442 return CxMatrixTem<C>(0, 0);
00443 }
00444
00445 int size = nRow();
00446 CxMatrixTem<C> m(size,size), tmp(*this);
00447 short* idx = new short [size];
00448 double d;
00449
00450 if (!ludcmp(tmp._data, size, idx, &d)) {
00451 CX_CERR << "Inverse: singular matrix can't be inverted." << CX_ENDL;
00452 delete [] idx;
00453 return CxMatrixTem<C>(0,0);
00454 }
00455
00456 double * res = new double [size];
00457
00458 for (int j = 0; j < size; j++) {
00459 for (int i = 0; i < size; i++)
00460 res[i] = 0.0;
00461 res[j] = 1.0;
00462 lubksb(tmp._data, size, idx, res);
00463 for (int i = 0; i < size; i++)
00464 m[i][j] = res[i];
00465 }
00466
00467 delete [] res;
00468 delete [] idx;
00469
00470 return m;
00471 }
00472
00473 #ifndef HxMatrix_EPS_DEF
00474 #define HxMatrix_EPS_DEF
00475 const double HxMatrix_EPS = 1e-12;
00476 #endif
00477
00478
00479
00480 template<class C>
00481 static int ludcmp(C* a, int n, short *indx, double *d)
00482 {
00483 double *vv = new double [n];
00484
00485 *d=1.0;
00486 int i, j, k;
00487 for ( i=0;i<n;i++) {
00488 double big=0.0;
00489 for (j=0;j<n;j++) {
00490 double temp = fabs(a[i*n+j]);
00491 if (temp > big)
00492 big=temp;
00493 }
00494 if (big == 0.0)
00495 return(0);
00496 vv[i]=1.0/big;
00497 }
00498 for (j=0;j<n;j++) {
00499 for (i=0;i<j;i++) {
00500 double sum=a[i*n+j];
00501 for (k=0;k<i;k++)
00502 sum -= a[i*n+k]*a[k*n+j];
00503 a[i*n+j]=sum;
00504 }
00505 double big=0.0;
00506 int imax = 0;
00507 for (i=j;i<n;i++) {
00508 double sum=a[i*n+j];
00509 for (k=0;k<j;k++)
00510 sum -= a[i*n+k]*a[k*n+j];
00511 a[i*n+j]=sum;
00512
00513 double dum = vv[i]*fabs(sum);
00514 if ( dum >= big) {
00515 big=dum;
00516 imax=i;
00517 }
00518 }
00519 if (j != imax) {
00520 for (k=0;k<n;k++) {
00521 double dum=a[imax*n+k];
00522 a[imax*n+k]=a[j*n+k];
00523 a[j*n+k]=dum;
00524 }
00525 *d = -(*d);
00526 vv[imax]=vv[j];
00527 }
00528 indx[j]=imax;
00529 if (a[j*n+j] == 0.0)
00530 a[j*n+j]=HxMatrix_EPS;
00531 if (j != n-1) {
00532 double dum=1.0/a[j*n+j];
00533 for (i=j+1;i<n;i++)
00534 a[i*n+j] *= dum;
00535 }
00536 }
00537 delete [] vv;
00538
00539 return(1);
00540 }
00541
00542 template<class C>
00543 static void lubksb(C* a, int n, short *indx, double *b)
00544 {
00545 int ii = -1;
00546
00547 int i;
00548 for (i=0;i<n;i++) {
00549 int ip=indx[i];
00550 double sum=b[ip];
00551 b[ip]=b[i];
00552 if (ii>=0)
00553 for (int j=ii;j<=i-1;j++)
00554 sum -= a[i*n+j]*b[j];
00555 else if (sum) ii=i;
00556 b[i]=sum;
00557 }
00558 for (i=n-1;i>=0;i--) {
00559 double sum=b[i];
00560 for (int j=i+1;j<n;j++)
00561 sum -= a[i*n+j]*b[j];
00562 b[i]=sum/a[i*n+i];
00563 }
00564 }
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579 template<class C> CxMatrixTem<C>
00580 CxMatrixTem<C>::covariance(CxVector<C> &mean) const
00581 {
00582
00583 int i,j,k;
00584
00585 CxMatrixTem<C> COVxy(_nc,_nc, 0.0) ;
00586
00587 mean = CxVector<C>(_nc) ;
00588 for (i=0; i<_nc; i++) mean[i] = 0;
00589
00590 for (k=0; k<_nr; k++)
00591 for (i=0; i<_nc; i++) {
00592 mean[i] += (*this)[k][i];
00593 for (j=i; j<_nc; j++)
00594 COVxy[i][j] += ((*this)[k][i] * (*this)[k][j]) ;
00595 }
00596
00597 for (i=0; i<_nc; i++) mean[i] /= _nr;
00598
00599 for (i=0; i<_nc; i++)
00600 for (j=i; j<_nc; j++) {
00601 COVxy[i][j] = COVxy[i][j]/_nr - mean[i]*mean[j];
00602 COVxy[j][i] = COVxy[i][j];
00603 }
00604
00605
00606 return COVxy;
00607 }
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627 template<class C>
00628 void tred2(CxMatrixTem<C> &a, CxVector<C> &d, CxVector<C> &e)
00629 {
00630 int i,j,k,l;
00631 double scale,hh,h,g,f;
00632
00633 int n = d.Size() ;
00634
00635 for (i=n-1; i>0; i--)
00636 {
00637 l = i-1;
00638 h = scale = 0.0;
00639 if (l>0)
00640 {
00641 for (k=0; k<=l; k++) scale += fabs(a[i][k]);
00642 if (scale == 0.0) e[i] = a[i][l];
00643 else
00644 {
00645 for (k=0; k<=l; k++)
00646 {
00647 a[i][k] /= scale;
00648 h += a[i][k]*a[i][k];
00649 }
00650 f = a[i][l];
00651 g = (f >=0.0 ? -sqrt(h) : sqrt(h));
00652 e[i] = scale*g;
00653 h -= f*g;
00654 a[i][l] = f-g;
00655 f = 0.0;
00656 for (j=0; j<=l; j++)
00657 {
00658 a[j][i] = a[i][j]/h;
00659 g = 0.0;
00660 for (k=0; k<=j; k++) g += a[j][k]*a[i][k];
00661 for (k=j+1; k<=l; k++) g += a[k][j]*a[i][k];
00662 e[j] = g/h;
00663 f += e[j]*a[i][j];
00664 }
00665 hh = f/(h+h);
00666 for (j=0; j<=l; j++)
00667 {
00668 f = a[i][j];
00669 e[j] = g = e[j]-hh*f;
00670 for (k=0; k<=j; k++) a[j][k] -= (f*e[k]+g*a[i][k]);
00671 }
00672 }
00673 }
00674 else e[i] = a[i][l];
00675 d[i] = h;
00676 }
00677
00678 d[0] = 0.0;
00679 e[0] = 0.0;
00680
00681 for (i=0; i<n; i++)
00682 {
00683 l = i-1;
00684 if (d[i])
00685 {
00686 for (j=0; j<=l; j++)
00687 {
00688 g = 0.0;
00689 for (k=0; k<=l; k++) g += a[i][k]*a[k][j];
00690 for (k=0; k<=l; k++) a[k][j] -= g*a[k][i];
00691 }
00692 }
00693 d[i] = a[i][i];
00694 a[i][i] = 1.0;
00695 for (j=0; j<=l; j++) a[j][i] = a[i][j] = 0.0;
00696 }
00697
00698 }
00699
00700 #define SIGN(a, b) ((b) < 0.0 ? -fabs(a): fabs(a))
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714 template<class C>
00715 void tqli(CxVector<C> &d, CxVector<C> &e, CxMatrixTem<C> &z)
00716 {
00717
00718 int m,l,iter,i,k;
00719 double s,r,p,g,f,dd,c,b;
00720
00721 int n = e.Size() ;
00722
00723 for (i=1; i<n; i++) e[i-1]=e[i];
00724 e[n-1] = 0.0;
00725
00726
00727 for (l=0; l<n; l++)
00728 {
00729 iter=0;
00730 do
00731 {
00732 for (m=l; m<n-1; m++)
00733 {
00734 dd = fabs(d[m]) + fabs(d[m+1]);
00735 if ( (fabs(e[m]) + dd) == dd) break;
00736 }
00737 if (m != l)
00738 {
00739 if (iter++ == 30) printf("Too many iterations in tqli\n");
00740 g = (d[l+1]-d[l])/(2.0*e[l]);
00741 r = sqrt(g*g + 1.0);
00742 g = d[m] - d[l] + e[l]/(g+SIGN(r,g));
00743 s = c = 1.0;
00744 p = 0.0;
00745 for (i=m-1; i>=l; i--)
00746 {
00747 f = s*e[i];
00748 b = c*e[i];
00749 e[i+1] = (r=sqrt(f*f+g*g));
00750 if (r == 0.0)
00751 {
00752 d[i+1] -= p;
00753 e[m] = 0.0;
00754 break;
00755 }
00756 s = f/r;
00757 c = g/r;
00758 g = d[i+1]-p;
00759 r = (d[i]-g)*s + 2.0*c*b;
00760 d[i+1] = g+(p=s*r);
00761 g = c*r-b;
00762
00763 for (k=0; k<n; k++)
00764 {
00765 f = z[k][i+1];
00766 z[k][i+1] = s*z[k][i] + c*f;
00767 z[k][i] = c*z[k][i] - s*f;
00768 }
00769 }
00770 if (r==0.0 && i>=l) continue;
00771 d[l] -= p;
00772 e[l] = g;
00773 e[m] = 0.0;
00774 }
00775 } while (m != l);
00776 }
00777
00778
00779 int j;
00780 for (i=0; i<n-1; i++)
00781 {
00782 p = d[k=i];
00783 for (j=i+1; j<n; j++) if (d[j]>=p) p = d[k=j];
00784 if (k!=i) {
00785 d[k] = d[i];
00786 d[i] = p;
00787
00788 for (j=0; j<n; j++)
00789 {
00790 p = z[j][i];
00791 z[j][i] = z[j][k];
00792 z[j][k] = p;
00793 }
00794 }
00795 }
00796 }
00797
00798
00799
00800
00801
00802 template<class C> CxMatrixTem<C>
00803 CxMatrixTem<C>::klm(int newDim) const
00804 {
00805 if (newDim<1 || newDim>_nc)
00806 {
00807 printf("klm: new dimensionality is invalid, dimsize is set as 1...\n");
00808 newDim = 1;
00809 }
00810
00811 CxMatrixTem<C> newfvec(_nr,newDim);
00812
00813 CxMatrixTem<C> a = this->covariance() ;
00814
00815 CxVector<C> d(_nc) ;
00816
00817 CxVector<C> e(_nc) ;
00818
00819 tred2(a,d,e);
00820
00821 tqli(d,e,a);
00822
00823
00824
00825 for (int k=0; k<_nr; k++)
00826 for (int i=0; i<newDim; i++) {
00827 double vl = 0.0;
00828 for (int j=0; j<_nc; j++) vl += ((*this)[k][j]*a[j][i]);
00829 newfvec[k][i] = vl;
00830 }
00831
00832 return newfvec;
00833 }
00834
00835 #endif