#include <doublematrix.h>
Inheritance diagram for matvec::doubleMatrix:
Public Types | |
| typedef int | size_type |
| typedef double | value_type |
| typedef double | element_type |
| typedef double * | pointer |
| typedef double * | iterator |
| typedef double & | reference |
| typedef const double * | const_iterator |
| typedef const double * | const_pointer |
| typedef const double & | const_reference |
| typedef std::reverse_iterator< iterator > | reverse_iterator |
| typedef std::reverse_iterator< const_iterator > | const_reverse_iterator |
Public Methods | |
| doubleMatrix (void) | |
| doubleMatrix (const size_type m, const size_type n) | |
| doubleMatrix (const size_type m, const size_type n, const double **a) | |
| doubleMatrix (const doubleMatrix &a) | |
| doubleMatrix (const Matrix< double > &a) | |
| doubleMatrix & | operator= (const Matrix< double > &a) |
| doubleMatrix & | operator= (const Matrix< bool > &a) |
| bool | psd (void) const |
| void | svd (doubleMatrix &u, Vector< double > &s, doubleMatrix &v) |
| int | rank (void) |
| doubleMatrix & | inv () |
| doubleMatrix | ginv0 (void) const |
| doubleMatrix | splines (const doubleMatrix knots, const unsigned type=0) const |
| doubleMatrix & | ginv1 (unsigned *irank=0) |
| doubleMatrix | mat_log (double tol=0) const |
| doubleMatrix | mat_exp (double tol=0) const |
| void | mat_exp_der (Vector< doubleMatrix > &der, double tol=0) const |
| doubleMatrix | covariance (const doubleMatrix &B) const |
| Vector< double > | variance (orientation orien=COLUMN) const |
| doubleMatrix | lu_solve (const doubleMatrix &rhs) |
| Vector< double > | lu_solve (const Vector< double > &rhs) |
| void | gs_solve (const doubleMatrix &v, doubleMatrix &solmat, const double relax=1.0, const double stopval=0.001, const int mxiter=1000) |
| void | gs_solve (const Vector< double > &v, Vector< double > &solvec, const double relax=1.0, const double stopval=0.001, const int mxiter=1000) |
| Vector< double > | eigen (const int job=1) |
| double | cond (void) |
| double | det (void) const |
| double | norm (const int p) const |
| double | norm (const std::string &s) const |
| double | logdet (void) |
| double | quadratic (const Vector< double > &a, const Vector< double > &b) const |
| doubleMatrix & | sqrtm (void) |
| doubleMatrix & | identity (const int m, const int n) |
| doubleMatrix & | identity (const int m) |
| doubleMatrix & | identity (void) |
| doubleMatrix & | sweep (const size_type i0, const size_type i1) |
| doubleMatrix & | sweep (void) |
| bool | symmetric (void) const |
| double ** | begin (void) |
| const double ** | begin (void) const |
| void | copy (const Matrix &a) |
| Matrix & | assign (const Matrix &a) |
| Matrix & | assign (const double &x) |
| Matrix & | resize (const size_type n, const size_type m, const double &val=double()) |
| Matrix & | resize (const size_type n, const size_type m, const double *a) |
| Matrix & | resize (const Matrix &a) |
| Matrix & | reserve (const size_type n, const size_type m) |
| Matrix | operator+ (const Matrix &a) const |
| Matrix | operator+ (const double &x) const |
| Matrix & | operator+ (void) |
| Matrix | operator- (const Matrix &a) const |
| Matrix | operator- (const double &x) const |
| Matrix | operator- (void) const |
| Matrix | operator * (const Matrix &a) const |
| Vector< double > | operator * (const Vector< double > &a) const |
| Matrix | operator * (const double &x) const |
| Matrix | operator/ (const Matrix &a) const |
| Matrix | operator/ (const double &x) const |
| Matrix & | operator+= (const Matrix &a) |
| Matrix & | operator+= (const double &x) |
| Matrix & | operator-= (const Matrix &a) |
| Matrix & | operator-= (const double &x) |
| Matrix & | operator *= (const Matrix &a) |
| Matrix & | operator *= (const double &x) |
| Matrix & | operator/= (const Matrix &a) |
| Matrix & | operator/= (const double &x) |
| double * | operator[] (const size_type i) |
| const double * | operator[] (const size_type i) const |
| double & | operator() (const size_type i, const size_type j) |
| Matrix< bool > | operator! (void) const |
| Matrix< bool > | operator== (const Matrix &a) const |
| Matrix< bool > | operator== (const double &x) const |
| Matrix< bool > | operator< (const Matrix &a) const |
| Matrix< bool > | operator< (const double &x) const |
| Matrix< bool > | operator> (const Matrix &a) const |
| Matrix< bool > | operator> (const double &x) const |
| Matrix< bool > | operator!= (const Matrix &a) const |
| Matrix< bool > | operator!= (const double &x) const |
| Matrix< bool > | operator<= (const Matrix &a) const |
| Matrix< bool > | operator<= (const double &x) const |
| Matrix< bool > | operator>= (const Matrix &a) const |
| Matrix< bool > | operator>= (const double &x) const |
| bool | all (void) const |
| bool | any (void) const |
| bool | empty (void) const |
| bool | save (const std::string &fname, const int io_mode=std::ios::out) const |
| Matrix | apply (double(*f)(double)) const |
| Matrix | apply (double(*f)(double, double), const Matrix &b) const |
| Matrix | apply (double(*f)(double, double), const double &b) const |
| Matrix | submat (const size_type idx1, const size_type idx2, const int len1=0, const int len2=0) const |
| Matrix | submat (const Matrix< bool > &a) const |
| Matrix | transpose (void) const |
| Matrix | kron (const Matrix &b) const |
| Matrix | reshape (const size_type m, const size_type n, orientation orient=COLUMN) const |
| Matrix | diag (void) const |
| Vector< double > | diag (const size_type k) const |
| Vector< double > | vec (orientation orient=COLUMN) const |
| Vector< double > | sum (orientation orient=COLUMN) const |
| Vector< double > | sumsq (orientation orient=COLUMN) const |
| Vector< double > | max (orientation orient=COLUMN) const |
| Vector< double > | min (orientation orient=COLUMN) const |
| virtual std::ostream & | print (std::ostream &os=std::cout) const |
| double | trace (double init=double()) const |
| double & | at (const size_type i, const size_type j) |
| Matrix & | sortby (orientation orient, const size_type k) |
| int | size (void) const |
| int | num_rows (void) const |
| int | num_cols (void) const |
| void | clear (void) |
| void | input (const std::string &fname, const size_type m, const size_type n) |
| void | initialize (const size_type m, const size_type n, const double **a) |
Public Attributes | |
| double ** | me |
| int | nrow |
| int | ncol |
Friends | |
| class | Model |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Definition at line 44 of file matrix.h. Referenced by doubleMatrix(). |
|
|
|
|
|
Definition at line 31 of file doublematrix.h.
00031 : Matrix<double>() {}
|
|
||||||||||||
|
Definition at line 32 of file doublematrix.h. References matvec::Matrix< double >::size_type.
00032 : Matrix<double>(m,n) {}
|
|
||||||||||||||||
|
Definition at line 33 of file doublematrix.h. References matvec::Matrix< double >::size_type.
00033 : Matrix<double>(m,n,a) {}
|
|
|
Definition at line 34 of file doublematrix.h.
00034 : Matrix<double>(a) { }
|
|
|
Definition at line 35 of file doublematrix.h.
00035 : Matrix<double>(a) { }
|
|
|
|
|
|
|
|
||||||||||||
|
|
|
||||||||||||
|
|
|
|
|
|
|
|
|
|
Definition at line 81 of file matrix.h. Referenced by matvec::GLMM::AI_REML(), matvec::GLMM::Build_SinMat(), matvec::GLMM::contrast(), matvec::Population::fullsibs_prob(), matvec::NuFamily::fullsibs_prob(), matvec::NuFamily::multi_fullsibs_prob(), matvec::Model::multipoint_estimate_Ve(), matvec::Model::multipoint_Rec_table(), matvec::GLMM::residual(), matvec::GLMM::SSQCP(), and matvec::GLMM::Vec2Var().
00081 {copy(a); return *this;}
|
|
||||||||||||
|
|
|
|
Definition at line 76 of file matrix.h.
00076 { return (const T**)me; }
|
|
|
Definition at line 75 of file matrix.h. Referenced by matvec::Model::add_Ag(), matvec::GLMM::add_Ag(), matvec::Model::add_Ag_diag(), matvec::GLMM::add_Ag_old(), matvec::GLMM::add_AgSand(), matvec::Model::add_Ig(), matvec::Model::add_Ig_diag(), matvec::GLMM::add_IgSand(), matvec::Population::anterior(), matvec::NuFamily::anterior(), matvec::Model::contrast(), matvec::GLMM::contrast(), matvec::Model::CovMat(), det(), matvec::Model::estimate(), matvec::Population::full_cdist(), matvec::Individual::get_penetrance(), ginv1(), identity(), matvec::Model::inverse_residual_var(), matvec::Population::llhood_phenotype(), matvec::Model::out_lsmeans_to_stream(), matvec::Population::partial_cdist(), psd(), matvec::Pedigree::reld(), and sqrtm().
00075 { return me; }
|
|
|
Referenced by matvec::Model::covariance(), and matvec::Model::multipoint_setup(). |
|
|
Definition at line 370 of file doublematrix.cpp. References matvec::Vector< T >::begin(), matvec::Matrix< T >::begin(), matvec::Vector< T >::max(), matvec::Matrix< double >::me, matvec::Vector< T >::min(), matvec::Matrix< double >::ncol, matvec::Matrix< double >::nrow, matvec::svdcmp(), and matvec::warning().
00371 {
00372 int nrow = Matrix<double>::nrow;
00373 int ncol = Matrix<double>::ncol;
00374 double** me = Matrix<double>::me;
00375
00376 if (nrow < ncol) throw exception("doubleMatrix::cond(), M.nrow() >= M.ncol()");
00377
00378 Vector<double> w(ncol);
00379 Matrix<double> v(ncol,ncol);
00380 Matrix<double> u(nrow,ncol);
00381
00382 for (int i=0; i<nrow; i++) memcpy(u[i],me[i],sizeof(double)*ncol);
00383 svdcmp(u.begin(), nrow, ncol, w.begin(), v.begin()); // Singular value decomposition.
00384
00385 double smax = w.max();
00386 double smin = w.min();
00387 if (smin == 0.0) warning("Matrix::cond(): condition is infinite");
00388 return (smax/smin);
00389 }
|
|
|
Referenced by matvec::Model::covariance(), and operator=(). |
|
|
Definition at line 678 of file doublematrix.cpp. References matvec::Matrix< double >::me, matvec::Matrix< double >::ncol, matvec::Matrix< double >::nrow, matvec::Matrix< double >::num_cols(), and matvec::Matrix< double >::num_rows().
00679 {
00680 int nrow = Matrix<double>::nrow;
00681 int ncol = Matrix<double>::ncol;
00682 double** me = Matrix<double>::me;
00683
00684 assert (nrow == B.num_rows() && ncol == B.num_cols());
00685 doubleMatrix temp(ncol,ncol);
00686 int i,j,k;
00687 double ss,sx,sy;
00688 for (i=0; i<ncol; ++i) {
00689 for (j=0; j<=i; ++j) {
00690 for (ss=0.0,sx=0.0,sy=0.0,k=0; k<nrow; ++k) {
00691 ss += me[k][i]*B[k][j];
00692 sx += me[k][i];
00693 sy += B[k][j];
00694 }
00695 temp[i][j] = temp[j][i] = (ss - sx*sy/ncol)/(ncol - 1);
00696 }
00697 }
00698 return temp;
00699 }
|
|
|
Definition at line 556 of file doublematrix.cpp. References matvec::Vector< T >::begin(), matvec::Matrix< double >::begin(), matvec::Session::epsilon, matvec::ludcmp(), matvec::Matrix< double >::ncol, matvec::Matrix< double >::nrow, and matvec::SESSION.
00557 {
00558 int nrow = Matrix<double>::nrow;
00559 int ncol = Matrix<double>::ncol;
00560
00561 double detval = 0.0;
00562 if ( nrow != ncol ) throw exception("doubleMatrix::det(): matrix must be square");
00563 int i,d;
00564 Vector<int> indx(nrow);
00565 doubleMatrix temp(*this);
00566 ludcmp(temp.begin(),nrow,indx.begin(),d, SESSION.epsilon);
00567 detval = (double)d;
00568 for (i=0; i<nrow; i++) detval *= temp[i][i];
00569 return detval;
00570 }
|
|
|
|
|
|
|
|
|
Definition at line 509 of file doublematrix.cpp. References matvec::Matrix< T >::begin(), matvec::Vector< T >::begin(), matvec::Matrix< double >::me, matvec::Matrix< double >::ncol, matvec::nonsymm_eigen(), matvec::Matrix< double >::nrow, matvec::Matrix< T >::reserve(), matvec::symm_eigen(), symmetric(), and matvec::warning(). Referenced by matvec::GLMM::AI_REML(), matvec::bound_pd(), matvec::GLMM::Build_SinMat(), mat_exp(), mat_exp_der(), mat_log(), and matvec::Model::vce_emreml_multi_trait().
00510 {
00511 ///////////////////////////////////////////////////////////////////
00512 // job = 0, no eigenvectors are needed
00513 // 1, eigenvectors are also needed.
00514 // return eigenvalues
00515 // Note that matrix itself become eigenvectors.
00516 // all imaginary numbers in eigenvalus/eigenvectors are ignored
00517 // if there are some for non-symmetrix real matrix.
00518 ////////////////////////////////////////////////////////////////////
00519
00520 // int nrow = Matrix<double>::nrow;
00521 // int ncol = Matrix<double>::ncol;
00522 // double** me = Matrix<double>::me;
00523
00524 if (nrow != ncol) throw exception( "doubleMatrix::eigen(): matrix must be square" );
00525 Vector<double> vals(nrow);
00526 if (this->symmetric()) {
00527 symm_eigen(vals.begin(),me,nrow);
00528 }
00529 else {
00530 std::cout << *this;
00531 Vector<double> ri(nrow);
00532 Matrix<double> vr;
00533 Matrix<double> vi;
00534 if (job == 1) {
00535 vr.reserve(nrow,nrow);
00536 vi.reserve(nrow,nrow);
00537 }
00538 int i = nonsymm_eigen(job,me,nrow,vals.begin(),ri.begin(),vr.begin(),vi.begin());
00539 if (i == -2) {
00540 warning(" Matrix::eigen(): out of memory");
00541 }
00542 else if (i == -1) {
00543 warning(" Matrix::eigen(): not converged");
00544 }
00545 else if (i == 1) {
00546 warning(" Matrix::eigen(): eigenvalues/vectors are complex,\n "
00547 " and the imaginary numbers are all ignored");
00548 }
00549 if (job == 1) {
00550 for (i=0; i<nrow; i++) memcpy(me[i],vr[i],sizeof(double)*ncol);
00551 }
00552 }
00553 return vals;
00554 }
|
|
|
Definition at line 138 of file matrix.h. Referenced by matvec::GLMM::SSQCP().
|
|
|
Definition at line 319 of file doublematrix.cpp. References matvec::Vector< T >::begin(), matvec::Matrix< T >::begin(), matvec::Session::epsilon, matvec::Vector< T >::max(), matvec::Matrix< double >::me, matvec::Matrix< double >::ncol, matvec::Matrix< double >::nrow, matvec::Matrix< T >::reserve(), matvec::Vector< T >::reserve(), matvec::SESSION, and matvec::svdcmp(). Referenced by matvec::GLMM::AI_REML(), matvec::GLMM::Build_SinMat(), matvec::GLMM::contrast(), and matvec::BG::NRupdate().
00320 { // Penrose-Moore inverse of any matrix
00321 int nrow = Matrix<double>::nrow;
00322 int ncol = Matrix<double>::ncol;
00323 double** me = Matrix<double>::me;
00324
00325 doubleMatrix temp(ncol,nrow);
00326 if (ncol==0 || nrow==0) return temp;
00327 Vector<double> w; w.reserve(ncol);
00328 Matrix<double> v; v.reserve(ncol,ncol);
00329 int m = nrow;
00330 if (nrow < ncol) m=ncol; //u(m,col) where m>=col
00331 Matrix<double> u; u.reserve(m,ncol);
00332
00333 int i,j,k;
00334 double x;
00335 for (i=0; i<nrow; i++) memcpy(u[i],me[i],sizeof(double)*ncol);
00336 for (i=0; i<m-nrow; i++) memset(u[nrow+i],'\0',sizeof(double)*ncol);
00337 svdcmp(u.begin(), m, ncol, w.begin(), v.begin()); // Singular value decomposition.
00338 double wmax = w.max(); // Maximum singular value.
00339 double wmin = wmax*SESSION.epsilon;
00340 for (k = 0; k < ncol; k++) {
00341 if (w[k] < wmin) {
00342 w[k] = 0.0;
00343 }
00344 else {
00345 w[k] = 1.0/w[k];
00346 }
00347 }
00348 for (i = 0; i < ncol; i++) {
00349 for (k = 0; k < ncol; k++) v[i][k] *= w[k];
00350 for (j = 0; j < nrow; j++) {
00351 for (x=0.0,k=0; k<ncol; ++k) x += v[i][k]*u[j][k];
00352 temp[i][j] = x;
00353 }
00354 }
00355 return temp;
00356 }
|
|
|
Definition at line 358 of file doublematrix.cpp. References matvec::Matrix< double >::begin(), matvec::Session::epsilon, matvec::ginverse1(), matvec::Matrix< double >::num_rows(), matvec::SESSION, and symmetric(). Referenced by matvec::GLMM::contrast(), matvec::GLMM::glim(), matvec::normal_loginClasses(), matvec::GLMM::residual(), matvec::Population::setup_m_mme(), splines(), and matvec::GLMM::SSQCP().
00359 { // it only works for symmetric-positive-semidefinite matrix
00360 // Matrix itself becomes its inverse;
00361 if (!this->symmetric()) {
00362 throw exception("doubleMatrix::ginv1(): matrix must be symmetric");
00363 }
00364 double lgdet;
00365 unsigned myrank = ginverse1(begin(),num_rows(),lgdet,1,SESSION.epsilon);
00366 if (irank) *irank = myrank;
00367 return *this;
00368 }
|
|
||||||||||||||||||||||||
|
Definition at line 492 of file doublematrix.cpp. References gs_solve(), matvec::Vector< T >::resize(), and matvec::Vector< T >::size().
00494 {
00495 int n = rhs.size();
00496 doubleMatrix bmat(n,1),solmat(n,1);
00497 if (sol.size() != n) {
00498 sol.resize(n);
00499 }
00500 int i;
00501 for (i=0; i<n; i++) {
00502 bmat[i][0] = rhs[i];
00503 solmat[i][0] = sol[i];
00504 }
00505 gs_solve(bmat,solmat,relax,stopval,mxiter);
00506 for (i=0; i<n; i++) sol[i] = solmat[i][0];
00507 }
|
|
||||||||||||||||||||||||
|
Definition at line 426 of file doublematrix.cpp. References matvec::Matrix< double >::diag(), matvec::Session::epsilon, matvec::Matrix< double >::me, matvec::Matrix< double >::ncol, matvec::Matrix< double >::nrow, matvec::Matrix< double >::num_cols(), matvec::Matrix< double >::num_rows(), matvec::Matrix< double >::resize(), matvec::SESSION, and matvec::warning(). Referenced by gs_solve().
00428 {
00429 // relax relaxation coefficient, relax=1.2 seems good for some case
00430 // stopval the accuracy at which iteration stops.
00431 // mxiter maximum number of iterations allowed
00432
00433 int nrow = Matrix<double>::nrow;
00434 int ncol = Matrix<double>::ncol;
00435 double** me = Matrix<double>::me;
00436
00437 int m = rhs.num_rows();
00438 int n = rhs.num_cols();
00439 if (nrow != ncol) throw exception("doubleMatrix::gs_solve(): matrix must be square ");
00440 if (ncol != m ) throw exception("doubleMatrix::gs_solve(): matrix and rhs are not conformable");
00441 if (solmat.num_rows() != m || solmat.num_cols() != n) {
00442 solmat.resize(m,n);
00443 }
00444 int i,j,k;
00445 double diag;
00446 int niter = 0;
00447 double a,cmax;
00448 double tol = SESSION.epsilon;
00449
00450 Vector<double> oldsol(n);
00451 Vector<double> newsol(n);
00452 Vector<double> cval(n);
00453 Vector<double> local(n);
00454 do { // now iteration begins
00455 cmax = 0.0;
00456 for (i=0; i<m; i++) {
00457 for (k=0; k<n; k++) {
00458 oldsol[k] = solmat[i][k];
00459 local[k] = 0.0;
00460 }
00461 for (j=0; j<m; j++) {
00462 a = me[i][j];
00463 if (j != i) for (k=0; k<n; k++) local[k] += a * solmat[j][k];
00464 }
00465 diag = me[i][i];
00466 if (diag > tol) {
00467 for (k=0; k<n; k++) {
00468 newsol[k] = (rhs[i][k] - local[k])/diag;
00469 cval[k] = (newsol[k] - oldsol[k])*relax;
00470 if (fabs(cval[k]) > cmax) cmax = fabs(cval[k]);
00471 newsol[k] = solmat[i][k] = oldsol[k] + cval[k];
00472 }
00473 }
00474 else if (diag > -tol) {
00475 throw exception("doubleMatrix::gs_solve(): zero-diagonal found in Matrix::gs_solve()");
00476 }
00477 else {
00478 throw exception("doubleMatrix::gs_solve(): matrix is not psd");
00479 }
00480 }
00481 niter += 1;
00482 if ( (niter % 10) ==0) {
00483 std::cout << " GS: # of iter = " << niter << ", max_change = "
00484 << cmax << std::endl;
00485 }
00486 } while (niter <= mxiter && cmax > stopval);
00487 if (cmax > stopval) {
00488 warning("doubleMatrix::gs_solve(): not converge: %20.10f",cmax);
00489 }
00490 }
|
|
|
Definition at line 69 of file doublematrix.h. References matvec::Matrix< double >::num_cols(), and matvec::Matrix< double >::num_rows(). Referenced by identity().
|
|
|
Definition at line 68 of file doublematrix.h. References identity().
00068 {return identity(m,m);}
|
|
||||||||||||
|
Definition at line 709 of file doublematrix.cpp. References matvec::Matrix< double >::begin(), matvec::Matrix< double >::me, and matvec::Matrix< double >::resize(). Referenced by matvec::GLMM::AI_REML(), matvec::Model::blup_pccg(), matvec::GLMM::new_obs(), matvec::Model::setup_mme(), matvec::Model::vce_dfreml(), and matvec::Model::vce_emreml().
|
|
||||||||||||||||
|
|
|
||||||||||||||||
|
|
|
|
Definition at line 166 of file doublematrix.cpp. References matvec::Vector< T >::assign(), matvec::Vector< T >::begin(), matvec::Session::epsilon, matvec::lubksb(), matvec::ludcmp(), matvec::Matrix< double >::me, matvec::Matrix< double >::ncol, matvec::Matrix< double >::nrow, matvec::Matrix< T >::reserve(), matvec::Vector< T >::reserve(), and matvec::SESSION. Referenced by matvec::GLMM::contrast(), splines(), and matvec::Model::vce_emreml_multi_trait().
00167 {
00168 //////////////////////////////////////////////////////////
00169 // A itself turns into its inverse
00170 // returns non_singular indicator
00171 // REQUIREMENTS: A must be nonsingular
00172 /////////////////////////////////////////////////////////
00173 int nrow = Matrix<double>::nrow;
00174 int ncol = Matrix<double>::ncol;
00175 double** me = Matrix<double>::me;
00176
00177 int i,j,d;
00178 if ( nrow != ncol ) throw exception("doubleMatrix::inv(): matrix must be square");
00179 Vector<int> indx; indx.reserve(nrow);
00180 Vector<double> tempcol; tempcol.reserve(nrow);
00181 Matrix<double> temp; temp.reserve(nrow,nrow);
00182
00183 ludcmp(me,nrow,indx.begin(),d, SESSION.epsilon);
00184 for (j=0; j<nrow; j++) {
00185 tempcol.assign(0.0);
00186 tempcol[j] = 1.0;
00187 lubksb(me,nrow,indx.begin(),tempcol.begin());
00188 for (i=0; i<nrow;i++) temp[i][j] = tempcol[i];
00189 }
00190 for (i=0;i<nrow; i++) memcpy(me[i],temp[i],sizeof(double)*nrow);
00191 return *this;
00192 }
|
|
|
|
|
|
Definition at line 617 of file doublematrix.cpp. References matvec::Matrix< T >::begin(), matvec::Session::epsilon, matvec::ginverse1(), matvec::Matrix< double >::me, matvec::Matrix< double >::ncol, matvec::Matrix< double >::nrow, matvec::SESSION, and symmetric(). Referenced by matvec::normal_loginClasses(), and matvec::GLMM::setup_ww().
00618 {
00619 int nrow = Matrix<double>::nrow;
00620 int ncol = Matrix<double>::ncol;
00621 double** me = Matrix<double>::me;
00622
00623 double retval = 0.0;
00624 if (!this->symmetric()) throw exception("doubleMatrix::logdet(): matrix must be symmetric");
00625 Matrix<double> temp(*this);
00626
00627 ginverse1(temp.begin(),nrow,retval,3,SESSION.epsilon);
00628 return retval;
00629 }
|
|
|
Definition at line 416 of file doublematrix.cpp. References lu_solve(), matvec::Vector< T >::size(), and matvec::Matrix< double >::vec().
00417 {
00418 int n = rhs.size();
00419 doubleMatrix bmat(n,1);
00420 for (int i=0; i<n; i++) bmat[i][0] = rhs[i];
00421 bmat = lu_solve(bmat);
00422 return bmat.vec();
00423 }
|
|
|
Definition at line 391 of file doublematrix.cpp. References matvec::Vector< T >::begin(), matvec::Session::epsilon, matvec::lubksb(), matvec::ludcmp(), matvec::Matrix< double >::me, matvec::Matrix< double >::ncol, matvec::Matrix< double >::nrow, matvec::Matrix< double >::num_cols(), matvec::Matrix< double >::num_rows(), matvec::Vector< T >::reserve(), matvec::SESSION, and matvec::warning(). Referenced by lu_solve().
00392 {
00393 int nrow = Matrix<double>::nrow;
00394 int ncol = Matrix<double>::ncol;
00395 double** me = Matrix<double>::me;
00396
00397 if (nrow != ncol) throw exception("doubleMatrix::lu_solve(): matrix is not square ");
00398 int m = rhs.num_rows();
00399 int n = rhs.num_cols();
00400
00401 if (ncol != m ) throw exception("doubleMatrix::lu_solve(): bad arg");
00402 doubleMatrix solmat(m,n);
00403 int i,j,d;
00404 Vector<double> solvec; solvec.reserve(m);
00405 Vector<int> indx; indx.reserve(m);
00406 ludcmp(me,nrow,indx.begin(),d,SESSION.epsilon);
00407 warning("A.solve(b), A destroyed");
00408 for (j=0; j<n; j++) {
00409 for (i=0; i<m; i++) solvec[i] = rhs[i][j];
00410 lubksb(me,m,indx.begin(),solvec.begin());
00411 for (i=0; i<m; i++) solmat[i][j] = solvec[i];
00412 }
00413 return solmat;
00414 }
|
|
|
Definition at line 227 of file doublematrix.cpp. References eigen(), matvec::Session::epsilon, matvec::Vector< T >::max(), matvec::Vector< T >::min(), matvec::Matrix< double >::nrow, matvec::Matrix< double >::resize(), matvec::SESSION, and symmetric(). Referenced by matvec::normal_loginClasses().
00228 { // Matrix Exponential
00229 if(tol==0) tol=SESSION.epsilon;
00230
00231 doubleMatrix temp,P;
00232 if (!this->symmetric()) throw exception("doubleMatrix.mat_exp() must be symmetric\n");
00233 if (nrow==0) return temp;
00234 temp.resize(nrow,nrow);
00235 Vector<double> D;
00236
00237 P=*this;
00238 D=P.eigen();
00239 double dmin=D.min();
00240 double dmax=D.max();
00241 if(dmin < dmax+2.*std::log(tol)) {
00242 dmin=dmax+2.*std::log(tol);
00243 }
00244 for(int k=0;k<nrow;k++) {
00245 if(D[k]<dmin) D[k]=dmin;
00246 D[k]=std::exp(D[k]);
00247 for(int i=0;i<nrow;i++){
00248 for(int j=0;j<nrow;j++){
00249 temp[i][j]+=P[i][k]*P[j][k]*D[k];
00250 }
00251 }
00252 }
00253 return temp;
00254 }
|
|
||||||||||||
|
Definition at line 259 of file doublematrix.cpp. References eigen(), matvec::Session::epsilon, matvec::Vector< T >::max(), matvec::Vector< T >::min(), matvec::Matrix< double >::nrow, matvec::Matrix< double >::resize(), matvec::Vector< T >::resize(), matvec::SESSION, and symmetric(). Referenced by matvec::normal_loginClasses().
00260 { // Matrix exp partial derivatives
00261 if(tol==0) tol=SESSION.epsilon;
00262
00263 doubleMatrix P;
00264 if (!this->symmetric()) throw exception("doubleMatrix.mat_exp() must be symmetric\n");
00265 if (nrow==0) return;
00266 int nvc=(nrow*(nrow+1))/2;
00267 temp.resize(nvc);
00268 for(int k=0;k<nvc;k++) temp[k].resize(nrow,nrow);
00269
00270 Vector<double> D;
00271
00272 P=*this;
00273 D=P.eigen();
00274 // cout << "\nP der\n" << P;
00275 double dmin=D.min();
00276 double dmax=D.max();
00277 if(dmin < dmax*2.*tol) {
00278 dmin=dmax*2.*tol;
00279 }
00280 for(int i=0;i<nrow;i++) {
00281 if(D[i]<dmin) D[i]=dmin;
00282 }
00283 doubleMatrix Partkl;
00284 double g;
00285 for(int k=0;k<nrow;k++){
00286 double dk=D[k];
00287 for(int l=0;l<nrow;l++) {
00288 double dl=D[l];
00289 if(std::fabs(dk-dl)< tol) g=std::sqrt(dk*dl);
00290 else g=(dk-dl)/std::log(dk/dl);
00291 Partkl.resize(nrow,nrow);
00292
00293 for(int i=0;i<nrow;i++)
00294 for(int j=0;j<nrow;j++)
00295 Partkl[i][j]=g*P[i][k]*P[j][l];
00296
00297
00298 int vc=0;
00299 for(int i=0;i<nrow;i++) {
00300 for(int j=i;j<nrow;j++,vc++) {
00301 temp[vc]+=P[i][k]*P[j][l]*Partkl;
00302 if(i!=j) temp[vc]+=P[j][k]*P[i][l]*Partkl;
00303 // cout << i << ":" << j << " "<< vc << " " << temp[vc] << endl;
00304 }
00305 }
00306
00307 }
00308 }
00309 // for(int i=0;i<nvc;i++) cout << i << "->" << temp[i] << endl;
00310
00311 return;
00312 }
|
|
|
Definition at line 197 of file doublematrix.cpp. References eigen(), matvec::Session::epsilon, matvec::Vector< T >::max(), matvec::Vector< T >::min(), matvec::Matrix< double >::nrow, matvec::Matrix< double >::resize(), matvec::SESSION, and symmetric(). Referenced by matvec::normal_loginClasses().
00198 { // Matrix Log
00199 if(tol==0) tol=SESSION.epsilon;
00200 doubleMatrix temp,P;
00201 if (!this->symmetric()) throw exception("doubleMatrix.mat_exp() must be symmetric\n");
00202 if (nrow==0) return temp;
00203 temp.resize(nrow,nrow);
00204 Vector<double> D;
00205
00206 P=*this;
00207 D=P.eigen();
00208 double dmin=D.min();
00209 double dmax=D.max();
00210 if(dmin < dmax*2.*tol) {
00211 dmin=dmax*2.*tol;
00212 }
00213 for(int k=0;k<nrow;k++) {
00214 if(D[k]<dmin) D[k]=dmin;
00215 D[k]=std::log(D[k]);
00216 for(int i=0;i<nrow;i++){
00217 for(int j=0;j<nrow;j++){
00218 temp[i][j]+=P[i][k]*P[j][k]*D[k];
00219 }
00220 }
00221 }
00222 return temp;
00223 }
|
|
|
|
|
|
|
|
|
Definition at line 601 of file doublematrix.cpp. References matvec::abs(), matvec::Matrix< double >::diag(), matvec::Matrix< double >::sum(), matvec::Field::sum(), and matvec::Matrix< double >::transpose().
00602 {
00603 double retval = 0.0;
00604 doubleMatrix A = this->transpose();
00605 if (s == "inf") {
00606 retval = abs(A).sum().max();
00607 }
00608 else if (s == "fro") {
00609 retval = std::sqrt((A*(*this)).diag(0).sum());
00610 }
00611 else {
00612 throw exception("doubleMatrix::norm((): bad arg, it must be either \"inf\" or \"fro\" ");
00613 }
00614 return retval;
00615 }
|
|
|
Definition at line 572 of file doublematrix.cpp. References matvec::abs(), matvec::Vector< T >::begin(), matvec::Matrix< T >::begin(), matvec::Vector< T >::max(), matvec::Matrix< double >::me, matvec::Matrix< double >::ncol, matvec::Matrix< double >::nrow, matvec::Matrix< T >::reserve(), matvec::Vector< T >::reserve(), matvec::Field::sum(), and matvec::svdcmp().
00573 {
00574 int nrow = Matrix<double>::nrow;
00575 int ncol = Matrix<double>::ncol;
00576 double** me = Matrix<double>::me;
00577
00578 double retval = 0.0;
00579 int i;
00580 if (p == 1) {
00581 doubleMatrix A(*this);
00582 retval = abs(A).sum().max();
00583 }
00584 else if (p == 2) {
00585 Vector<double> w; w.reserve(ncol);
00586 Matrix<double> v; v.reserve(ncol,ncol);
00587 int m = nrow;
00588 if (nrow < ncol) m = ncol; //u(m,ccol) where m>=ccol
00589 Matrix<double> u(m,ncol);
00590
00591 for (i=0; i<nrow; i++) memcpy(u[i],me[i],sizeof(double)*ncol);
00592 svdcmp(u.begin(), m, ncol, w.begin(), v.begin()); // Singular value decomposition.
00593 retval = w.max();
00594 }
00595 else {
00596 throw exception("doubleMatrix::norm(): bad arg, 1 or 2 expected");
00597 }
00598 return retval;
00599 }
|
|
|
Definition at line 163 of file matrix.h. Referenced by matvec::BGMatrix::BGMatrix(), matvec::Model::contrast(), matvec::GLMM::contrast(), covariance(), matvec::Model::CovMat(), matvec::Model::estimate(), gs_solve(), identity(), lu_solve(), matvec::Plotter::plot(), matvec::Plotter::plot3D(), matvec::PoissonDist::sample(), matvec::UniformDist::sample(), and sweep().
00163 { return ncol; }
|
|
|
Definition at line 162 of file matrix.h. Referenced by matvec::BGMatrix::BGMatrix(), matvec::Model::blup_pccg(), matvec::Model::contrast(), matvec::GLMM::contrast(), covariance(), matvec::Model::CovMat(), matvec::Model::estimate(), ginv1(), gs_solve(), identity(), lu_solve(), matvec::PoissonDist::sample(), matvec::UniformDist::sample(), matvec::Model::setup_mme(), sqrtm(), matvec::GLMM::SSQCP(), sweep(), matvec::UnknownDist::UnknownDist(), matvec::Model::vce_dfreml(), and matvec::Model::vce_emreml().
00162 { return nrow; }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||||
|
|
|
|
Definition at line 117 of file matrix.h.
00117 { return *this; } //unary plus
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Definition at line 39 of file doublematrix.cpp. References matvec::Matrix< T >::num_cols(), matvec::Matrix< T >::num_rows(), and matvec::Matrix< double >::resize().
00040 {
00041 resize(a.num_rows(),a.num_cols());
00042 for (int j,i=0; i<a.num_rows(); ++i) for (j=0; j<a.num_cols(); ++j) Matrix<double>::me[i][j] = (double)a[i][j];
00043 return *this;
00044 }
|
|
|
Definition at line 37 of file doublematrix.h. References matvec::Matrix< double >::copy().
00037 { copy(a); return *this; }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Definition at line 114 of file matrix.h.
00114 { return me[i]; }
|
|
|
Definition at line 113 of file matrix.h.
00113 { return me[i]; }
|
|
|
Referenced by matvec::Model::info(), and matvec::Model::vce_emreml_multi_trait(). |
|
|
Definition at line 47 of file doublematrix.cpp. References matvec::Matrix< double >::begin(), matvec::Session::epsilon, matvec::Matrix< double >::ncol, matvec::Matrix< double >::nrow, matvec::psdefinite(), matvec::SESSION, and symmetric(). Referenced by matvec::funk(), matvec::Model::minfun_vce(), and matvec::UnknownDist::UnknownDist().
|
|
||||||||||||
|
Definition at line 631 of file doublematrix.cpp. References matvec::Matrix< double >::me, matvec::Matrix< double >::ncol, matvec::Matrix< double >::nrow, and matvec::Vector< T >::size(). Referenced by matvec::GLMM::contrast(), matvec::normal_loginClasses(), matvec::GLMM::residual(), and matvec::GLMM::SSQCP().
00632 {
00633 int nrow = Matrix<double>::nrow;
00634 int ncol = Matrix<double>::ncol;
00635 double** me = Matrix<double>::me;
00636 int m = x.size();
00637 int n = y.size();
00638
00639 double s = 0.0;
00640 if (m != nrow || n != ncol) throw exception("doubleMatrix::quadratic(): size not conformable");
00641 for (int j,i=0;i<m; i++) for (j=0; j<n; j++) s += x[i]*me[i][j]*y[j];
00642 return s;
00643 }
|
|
|
Definition at line 69 of file doublematrix.cpp. References matvec::Vector< T >::begin(), matvec::Matrix< T >::begin(), matvec::Session::epsilon, matvec::SESSION, and matvec::svdcmp().
00070 {
00071 int nrow = Matrix<double>::nrow;
00072 int ncol = Matrix<double>::ncol;
00073 double** me = Matrix<double>::me;
00074
00075 Vector<double> w(ncol);
00076 Matrix<double> v(ncol,ncol);
00077 int m = nrow;
00078 if (nrow < ncol) m = ncol; //u(m,col) where m>=col
00079 Matrix<double> u(m,ncol);
00080 int i;
00081 for (i=0; i<nrow; i++) memcpy(u[i],me[i],sizeof(double)*ncol);
00082 svdcmp(u.begin(), m, ncol, w.begin(), v.begin()); // Singular value decomposition.
00083
00084 int r = 0;
00085 for (i=0; i<ncol; i++) if (w[i] > SESSION.epsilon) r++;
00086 return (r);
00087 }
|
|
||||||||||||
|
Referenced by matvec::Model::copyfrom(), and matvec::Model::prepare_data(). |
|
||||||||||||||||
|
|
|
|
Definition at line 86 of file matrix.h.
00086 {resize(a.ne); return *this;}
|
|
||||||||||||||||
|
|
|
||||||||||||||||
|
||||||||||||
|
Referenced by matvec::Plotter::plot(), and matvec::Plotter::plot3D(). |
|
|
Definition at line 161 of file matrix.h. Referenced by matvec::Population::sum_genotype_freq1(), and matvec::Population::sum_genotype_freq2().
|
|
||||||||||||
|
Referenced by matvec::KP(). |
|
||||||||||||
|
Definition at line 93 of file doublematrix.cpp. References ginv1(), inv(), matvec::Matrix< double >::me, matvec::Matrix< double >::nrow, matvec::Matrix< double >::resize(), sqrtm(), and matvec::Matrix< double >::transpose().
00095 {
00096 int nk=knots.nrow;
00097 int n=nrow;
00098 double k;
00099 doubleMatrix R,Q,Delt,H,X;
00100 R.resize(nk,nk);
00101 Q.resize(nk,nk);
00102 H.resize(nk-1,1);
00103 k=(knots[nk-1][0]-knots[0][0])/((double)(nk-1));
00104 for(int i=0;i<nk-1;i++) H[i][0]=(knots[i+1][0]-knots[i][0]);
00105 for(int i=2;i<nk;i++){
00106 R(i,i)=(H(i,1)+H(i-1,1))/3;
00107 if (i<(nk-1)) {
00108 R(i+1,i)=H(i,1)/6;
00109 R(i,i+1)=H(i,1)/6;
00110 }
00111 Q(i-1,i)=1/H(i-1,1);
00112 Q(i,i)=-(1/H(i-1,1))-(1/H(i,1));
00113 Q(i+1,i)=1/H(i,1);
00114 }
00115 doubleMatrix Xg,Xgamma;
00116 Xg.resize(n,nk);
00117 Xgamma.resize(n,nk);
00118 for(int i=1;i<=n;i++){
00119 double t=me[i-1][0];
00120 int j;
00121 for( j=1;t>knots.me[j][0];j++);
00122 double th=knots.me[j][0];
00123 double tl=knots.me[j-1][0];
00124 Xg(i,j+1)=(t-tl)/H(j,1);
00125 Xg(i,j)=(th-t)/H(j,1);
00126 Xgamma(i,j+1)=-(t-tl)*(th-t)*(1+(t-tl)/H(j,1))/6;
00127 Xgamma(i,j)=-(t-tl)*(th-t)*(1+(th-t)/H(j,1))/6;
00128 }
00129 X=Xg+Xgamma*R.ginv1()*Q.transpose();
00130 if(type==1) {
00131 doubleMatrix Xs,Delt,G;
00132 G.resize(nk-2,nk-2);
00133 Xs.resize(nk,nk);
00134 Delt.resize(nk,nk-2);
00135 for(int i=0;i<nk-2;i++){
00136 Delt[i][i]=1./H[i][0];
00137 Delt[i+1][i]=-(1./H[i][0]+1./H[i+1][0]);
00138 Delt[i+2][i]=1./H[i+1][0];
00139 G[i][i]=(H[i][0]+H[i+1][0])/3.;
00140 if(i) {
00141 G[i-1][i]=H[i][0]/6.;
00142 G[i][i-1]=H[i][0]/6.;
00143 }
00144 }
00145 G.sqrtm();
00146 doubleMatrix Zs,dpd;
00147 dpd=Delt.transpose()*Delt;
00148 Zs=Delt*dpd.inv()*G;
00149 Zs/=pow(k,1.5);
00150 for(int i=0;i<nk;i++){
00151 Xs[i][0]=1.;
00152 Xs[i][1]=knots[i][0];
00153 for(int j=0;j<nk-2;j++){
00154 Xs[i][j+2]=Zs[i][j];
00155 }
00156
00157 }
00158 X=X*Xs;
00159 }
00160 //cout << X;
00161
00162 return(X);
00163 }
|
|
|
Definition at line 701 of file doublematrix.cpp. References matvec::Matrix< double >::begin(), matvec::Session::epsilon, matvec::ginverse1(), matvec::Matrix< double >::num_rows(), matvec::SESSION, and symmetric(). Referenced by splines(), and matvec::Model::vce_emreml_multi_trait().
|
|
|
|
|
||||||||||||||||||||
|
Referenced by matvec::GLMM::contrast(), and matvec::KP(). |
|
|
Referenced by matvec::Individual::get_m_penetrance(), and norm(). |
|
|
|
|
||||||||||||||||
|
Definition at line 53 of file doublematrix.cpp. References matvec::Vector< T >::begin(), matvec::Matrix< double >::me, matvec::Matrix< double >::ncol, matvec::Matrix< double >::nrow, matvec::Vector< T >::resize(), matvec::Matrix< double >::resize(), matvec::Vector< T >::size(), and matvec::svdcmp().
00054 {
00055 int nrow = Matrix<double>::nrow;
00056 int ncol = Matrix<double>::ncol;
00057 double** me = Matrix<double>::me;
00058
00059 if (nrow < ncol) throw exception("doubleMatrix::svd(), bad args");
00060
00061 if (u.nrow != nrow || u.ncol != ncol) u.resize(nrow,ncol);
00062 if (s.size() != ncol ) s.resize(ncol);
00063 if (v.nrow != ncol || v.ncol != ncol) v.resize(ncol,ncol);
00064
00065 for (int i=0; i<nrow; i++) memcpy(u[i],me[i],sizeof(double)*ncol);
00066 svdcmp(u.me, nrow, ncol, s.begin(), v.me); // Singular value decomposition.
00067 }
|
|
|
Definition at line 71 of file doublematrix.h. References matvec::Matrix< double >::num_cols(), and matvec::Matrix< double >::num_rows().
|
|
||||||||||||
|
|
|
|
Reimplemented from matvec::Matrix< double >. Definition at line 72 of file doublematrix.h. References matvec::Session::epsilon, matvec::Matrix< double >::me, matvec::Matrix< double >::nrow, and matvec::SESSION. Referenced by eigen(), ginv1(), logdet(), mat_exp(), mat_exp_der(), mat_log(), psd(), and sqrtm().
|
|
|
|
|
|
Referenced by matvec::GLMM::AI_REML(), matvec::bound_pd(), matvec::GLMM::Build_SinMat(), matvec::GLMM::contrast(), matvec::exp(), matvec::log(), norm(), matvec::normal_loginClasses(), matvec::operator *(), matvec::operator/(), matvec::GLMM::residual(), splines(), matvec::sqrt(), matvec::GLMM::SSQCP(), matvec::GLMM::Var2Vec(), matvec::Model::vce_emreml_multi_trait(), and matvec::GLMM::Vec2Var(). |
|
|
Definition at line 644 of file doublematrix.cpp. References matvec::Matrix< double >::COLUMN, matvec::Matrix< double >::me, matvec::Matrix< double >::ncol, matvec::Matrix< double >::nrow, matvec::Vector< T >::resize(), matvec::Matrix< double >::ROW, and matvec::warning().
00645 {
00646 int nrow = Matrix<double>::nrow;
00647 int ncol = Matrix<double>::ncol;
00648 double** me = Matrix<double>::me;
00649 int i,j;
00650 double ss,s;
00651 Vector<double> temp;
00652 if (orien == COLUMN) {
00653 assert(nrow > 1);
00654 temp.resize(ncol);
00655 for (j=0; j<ncol; j++) {
00656 for (ss=0.0,s=0.0,i=0; i<nrow; ++i) {
00657 ss += me[i][j]*me[i][j];
00658 s += me[i][j];
00659 }
00660 temp[j] = (ss-s*s/nrow)/(nrow - 1);
00661 }
00662 } else if (orien == ROW) {
00663 assert(ncol > 1);
00664 temp.resize(nrow);
00665 for (i=0; i<nrow; i++) {
00666 for (ss=0.0,s=0.0,j=0; j<ncol; ++j) {
00667 ss += me[i][j]*me[i][j];
00668 s += me[i][j];
00669 }
00670 temp[i] = (ss-s*s/ncol)/(ncol - 1);
00671 }
00672 } else {
00673 warning("unknown orientation");
00674 }
00675 return temp;
00676 }
|
|
|
Referenced by lu_solve(). |
|
|
Definition at line 29 of file doublematrix.h. |
|
|
Definition at line 170 of file matrix.h. Referenced by cond(), matvec::Model::covariance(), covariance(), eigen(), ginv0(), gs_solve(), identity(), inv(), logdet(), lu_solve(), norm(), quadratic(), splines(), svd(), symmetric(), and variance(). |
|
|
Definition at line 171 of file matrix.h. Referenced by cond(), covariance(), det(), eigen(), ginv0(), gs_solve(), inv(), logdet(), lu_solve(), norm(), psd(), quadratic(), svd(), and variance(). |
|
|
Definition at line 171 of file matrix.h. Referenced by cond(), covariance(), det(), eigen(), ginv0(), gs_solve(), inv(), logdet(), lu_solve(), mat_exp(), mat_exp_der(), mat_log(), norm(), psd(), quadratic(), splines(), svd(), symmetric(), and variance(). |
1.2.16