C++ Matrix::Column方法代码示例

/** 
 * Compute a distance metric between two columns of a
 * matrix. <b>Note that the indexes are *1* based (not 0) as that is
 * Newmat's convention</b>. Note that dist(M,i,j) must equal dist(M,j,i);
 * 
 * @param M - Matrix whose columns represent individual items to be clustered.
 * @param col1Ix - Column index to be compared (1 based).
 * @param col2Ix - Column index to be compared (1 based).
 * 
 * @return - "Distance" or "dissimilarity" metric between two columns of matrix.
 */
double GuassianRadial::dist(const Matrix &M, int col1Ix, int col2Ix) const {
  double dist = 0;
  if(col1Ix == col2Ix) 
    return 0;
  ColumnVector V = M.Column(col1Ix) - M.Column(col2Ix);
  dist = V.SumSquare() / (2 * m_Sigma * m_Sigma);
  dist = exp(-1 * dist);
  return dist;
}

ReturnMatrix Helmert(const Matrix& X, bool full)
{
   REPORT
   Tracer et("Helmert * Matrix");
   int m = X.nrows(); int n = X.ncols();
   if (m == 0) Throw(ProgramException("Matrix has 0 rows ", X));
   Matrix Y;
   if (full) Y.resize(m,n); else Y.resize(m-1, n);
   for (int j = 1; j <= n; ++j)
   {
      ColumnVector CV = X.Column(j);
      Y.Column(j) = Helmert(CV, full);
   }
   Y.release(); return Y.for_return();
}

std::vector<int> fiedler_reorder(const SymmetricMatrix& m)
{
  SymmetricMatrix absm=m;
  const int nrows=m.Nrows();
  for (int i=0;i<nrows;++i) {
    for (int j=0;j<=i;++j){
       //absolute value
       absm.element(i,j)=std::fabs(absm.element(i,j));
    }
  }

  //laplacian
  SymmetricMatrix lap(nrows);
  lap=0.;
  for (int i=0;i<nrows;++i)
    lap.element(i,i)=absm.Row(i+1).Sum();
  lap-=absm;

  DiagonalMatrix eigs; 
  Matrix vecs;
  
  EigenValues(lap,eigs,vecs);

  ColumnVector fvec=vecs.Column(2);
  std::vector<double> fvec_stl(nrows);
  //copies over fvec to fvec_stl
  std::copy(&fvec.element(0),&fvec.element(0)+nrows,fvec_stl.begin());
  std::vector<int> findices;
  //sorts the data by eigenvalue in ascending order
  sort_data_to_indices(fvec_stl,findices);
  
  return findices;
  /* BLOCK works with findices*/

}

queue<vector<bool>> Filter_Indep_Col(queue<vector<bool>> q, Matrix A)
{
	queue<vector<bool>> indep_q;
	int nrow_A = A.Nrows();
	int ncol_A = A.Ncols();
	int rank_A = Rank(A);
	int size = q.size();

	while (size > 0)
	{
		Matrix B(nrow_A, rank_A);
		int k = 1;
		for (int j = 0; j < ncol_A; j++)
		{
			if (q.front().at(j) == true)
			{
				B.Column(k) = A.Column(j+1);
				k++;
			}
		
		}

		// see if columns of B are indepenent
		if (Rank(B) == min(B.Ncols(), B.Nrows()))
		{
			indep_q.push(q.front());
		}
		q.pop();
		size--;
	}
	return indep_q;
}

Matrix Find_T_BFS(queue<vector<bool>> q, Matrix S_BFS)
{
	int rank_S = Rank(S_BFS);
	int nrow_S = S_BFS.Nrows();
	int ncol_S = S_BFS.Ncols();
	int size = q.size();
	Matrix T(nrow_S, rank_S);
	T = 0;
	while (size > 0)
	{
		int r = 1;
		for (int i = 1; i < ncol_S; i++)
		{
			if (q.front().at(i) == true)
			{
				T.Column(r) = S_BFS.Column(i);
				r++;
			}

		}
		q.pop();
		size--;
		if (Rank(T) == rank_S)
			break;
	}
	return T;
}

// Matrix A's first n columns are orthonormal
// so A.Columns(1,n).t() * A.Columns(1,n) is the identity matrix.
// Fill out the remaining columns of A to make them orthonormal
// so A.t() * A is the identity matrix 
void extend_orthonormal(Matrix& A, int n)
{
   REPORT
   Tracer et("extend_orthonormal");
   int nr = A.nrows(); int nc = A.ncols();
   if (nc > nr) Throw(IncompatibleDimensionsException(A));
   if (n > nc) Throw(IncompatibleDimensionsException(A));
   ColumnVector SSR;
   { Matrix A1 = A.Columns(1,n); SSR = A1.sum_square_rows(); }
   for (int i = n; i < nc; ++i)
   {
      // pick row with smallest SSQ
      int k; SSR.minimum1(k);
      // orthogonalise column with 1 at element k, 0 elsewhere
      // next line is rather inefficient
      ColumnVector X = - A.Columns(1, i) * A.SubMatrix(k, k, 1, i).t();
      X(k) += 1.0;
      // normalise
      X /= sqrt(X.SumSquare());
      // update row sums of squares
      for (k = 1; k <= nr; ++k) SSR(k) += square(X(k));
      // load new column into matrix
      A.Column(i+1) = X;
   }
}

//Predict on a chunk of data.
ReturnMatrix SOGP::predictM(const Matrix& in, ColumnVector &sigconf,bool conf){
  //printf("SOGP::Predicting on %d points\n",in.Ncols());
  Matrix out(alpha.Ncols(),in.Ncols());
  sigconf.ReSize(in.Ncols());
  for(int c=1;c<=in.Ncols();c++)
    out.Column(c) = predict(in.Column(c),sigconf(c),conf);
  out.Release();
  return out;
}

void CircularShift(const Matrix& X1, int first, int last)
{
      Matrix X; UpperTriangularMatrix U1, U2;
      int n = X1.Ncols();

      // Try right circular shift of columns
      X = X1; QRZ(X, U1);
      RightCircularUpdateCholesky(U1, first, last);
      X = X1.Columns(1,first-1) | X1.Column(last)
         | X1.Columns(first,last-1) | X1.Columns(last+1,n);
      QRZ(X, U2);
      X = U1 - U2; Clean(X, 0.000000001); Print(X);

      // Try left circular shift of columns
      X = X1; QRZ(X, U1);
      LeftCircularUpdateCholesky(U1, first, last);
      X = X1.Columns(1,first-1) | X1.Columns(first+1,last)
         | X1.Column(first) | X1.Columns(last+1,n);
      QRZ(X, U2);
      X = U1 - U2; Clean(X, 0.000000001); Print(X);
}

ReturnMatrix Helmert_transpose(const Matrix& Y, bool full)
{
   REPORT
   Tracer et("Helmert_transpose * Matrix ");
   int m = Y.nrows(); int n = Y.ncols(); if (!full) ++m;
   Matrix X(m, n);
   for (int j = 1; j <= n; ++j)
   {
      ColumnVector CV = Y.Column(j);
      X.Column(j) = Helmert_transpose(CV, full);
   }
   X.release(); return X.for_return();
}

Matrix Find_S_BFS(queue<vector<bool>> q, Matrix A, ColumnVector y)
{
	int nrow_A = A.Nrows();
	int ncol_A = A.Ncols();
	int rank_A = Rank(A);
	int size = q.size();

	Matrix S_BFS(ncol_A, size);
	S_BFS = 0;
	int S_j = 1; // column index for S
	while (size > 0)
	{
		int k = 1;
		Matrix B(nrow_A, rank_A);
		for (int j = 0; j < ncol_A; j++)
		{
			if (q.front().at(j) == true)
			{
				for (int i = 1; i <= nrow_A; i++)
				{
					B.Column(k) = A.Column(j + 1);
				}
				k++;
			}
		}
		// Find B already, and then find solution
		ColumnVector s(rank_A);
        s = Find_Solution(B, y);

	    // Insert zero into appropriate positions
	    // Find S_BFS matrix
		int kk = 1;
		ColumnVector soln(ncol_A);
		for (int i = 0; i < ncol_A; i++)
		{
			if (q.front().at(i) == true)
			{
				S_BFS(i + 1, S_j) = s(kk);
				kk++;
			}
			else
			{
				S_BFS(i + 1, S_j) = 0;
			}
		}
		q.pop();
		size--;
		S_j++;
	}
	return S_BFS;
}

bool SpectClust::findNLargestSymEvals(const SymmetricMatrix &W, int numLamda, std::vector<Numeric> &eVals, Matrix &EVec) {
  bool converged = false;
  eVals.clear();
  eVals.reserve(numLamda);
  DiagonalMatrix D(W.Ncols());
  Matrix E;
  try {
    EigenValues(W, D, E);
    converged = true;
    EVec.ReSize(W.Ncols(), numLamda);
    int count = 0;
    for(int i = W.Ncols(); i > W.Ncols() - numLamda; i--) {
      eVals.push_back(D(i));
      EVec.Column(++count) << E.Column(i);
    }
  }
  catch(const Exception &e) {
    Err::errAbort("Exception: " + ToStr(e.what()));
  }
  catch(...) {
    Err::errAbort("Yikes couldn't calculate eigen vectors.");
  }
  return converged;
}

Matrix& operator *(const Matrix& a, const Matrix& b) {
	if (a.Column() != b.Row()) {
		throw "Exception : Invailed matrix size.\n";
		exit(EXIT_FAILURE);
	}

	Matrix *re = new Matrix(a.Row(), b.Column());
	for (int i = 0; i < re->Row(); i++) {
		for (int j = 0; j < re->Column(); j++) {
			register double f = 0.0;
			for (int k = 0; k < a.Column(); k++) {
				f += a(i, k) * b(k, j);
			}
			(*re)(i, j) = f;
		}
	}
	return *re;
}

// produces the Cholesky decomposition of EAE where A = chol.t() * chol
// and E produces a LEFT circular shift of the rows and columns from
// 1,...,k-1,k,k+1,...l,l+1,...,p to
// 1,...,k-1,k+1,...l,k,l+1,...,p to
void left_circular_update_Cholesky(UpperTriangularMatrix &chol, int k, int l)
{
   int nRC = chol.Nrows();
   int i, j;

   // I. compute shift of column k to the lth position
   Matrix cholCopy = chol;
   // a. grab column k
   ColumnVector columnK = cholCopy.Column(k);
   // b. shift columns k+1,...l to the LEFT
   for(j = k+1; j <= l; ++j)
      cholCopy.Column(j-1) = cholCopy.Column(j);
   // c. copy the elements of columnK into the lth column of cholCopy
   cholCopy.Column(l) = 0.0;
   for(i = 1; i <= k; ++i)
      cholCopy(i,l) = columnK(i);

   // II. apply and compute Given's rotations
   int nGivens = l-k;
   ColumnVector cGivens(nGivens); cGivens = 0.0;
   ColumnVector sGivens(nGivens); sGivens = 0.0;
   for(j = k; j <= nRC; ++j)
   {
      ColumnVector columnJ = cholCopy.Column(j);

      // apply the previous Givens rotations to columnJ
      int imax = j - k; if (imax > nGivens) imax = nGivens;
      for(int i = 1; i <= imax; ++i)
      {
         int gIndex = i;
         int topRowIndex = k + i - 1;
         GivensRotationR(cGivens(gIndex), sGivens(gIndex),
            columnJ(topRowIndex), columnJ(topRowIndex+1));
      }

      // compute a new Given's rotation when j < l
      if(j < l)
      {
         int gIndex = j-k+1;
         columnJ(j) = pythag(columnJ(j), columnJ(j+1), cGivens(gIndex),
            sGivens(gIndex));
         columnJ(j+1) = 0.0;
      }

      cholCopy.Column(j) = columnJ;
   }

   chol << cholCopy;
	
}

// produces the Cholesky decomposition of EAE where A = chol.t() * chol
// and E produces a RIGHT circular shift of the rows and columns from
// 1,...,k-1,k,k+1,...l,l+1,...,p to
// 1,...,k-1,l,k,k+1,...l-1,l+1,...p
void right_circular_update_Cholesky(UpperTriangularMatrix &chol, int k, int l)
{
   int nRC = chol.Nrows();
   int i, j;
	
   // I. compute shift of column l to the kth position
   Matrix cholCopy = chol;
   // a. grab column l
   ColumnVector columnL = cholCopy.Column(l);
   // b. shift columns k,...l-1 to the RIGHT
   for(j = l-1; j >= k; --j)
      cholCopy.Column(j+1) = cholCopy.Column(j);
   // c. copy the top k-1 elements of columnL into the kth column of cholCopy
   cholCopy.Column(k) = 0.0;
   for(i = 1; i < k; ++i) cholCopy(i,k) = columnL(i);

    // II. determine the l-k Given's rotations
   int nGivens = l-k;
   ColumnVector cGivens(nGivens); cGivens = 0.0;
   ColumnVector sGivens(nGivens); sGivens = 0.0;
   for(i = l; i > k; i--)
   {
      int givensIndex = l-i+1;
      columnL(i-1) = pythag(columnL(i-1), columnL(i),
         cGivens(givensIndex), sGivens(givensIndex));
      columnL(i) = 0.0;
   }
   // the kth entry of columnL is the new diagonal element in column k of cholCopy
   cholCopy(k,k) = columnL(k);
	
   // III. apply these Given's rotations to subsequent columns
   // for columns k+1,...,l-1 we only need to apply the last nGivens-(j-k) rotations
   for(j = k+1; j <= nRC; ++j)
   {
      ColumnVector columnJ = cholCopy.Column(j);
      int imin = nGivens - (j-k) + 1; if (imin < 1) imin = 1;
      for(int gIndex = imin; gIndex <= nGivens; ++gIndex)
      {
         // apply gIndex Given's rotation
         int topRowIndex = k + nGivens - gIndex;
         GivensRotationR(cGivens(gIndex), sGivens(gIndex),
            columnJ(topRowIndex), columnJ(topRowIndex+1));
      }
      cholCopy.Column(j) = columnJ;
   }

   chol << cholCopy;
}

bool GenSetBase::generateAll(Matrix& M, ColumnVector& X, double D){ 
  if (Size<=0 || Vdim<=0) {
    cerr << "***ERROR: GenSetBase::generateAll(Matrix,...) "
	 << "called with size=" << Size << ", vdim=" << Vdim << endl;
    return false;
  }
  if (M.Ncols() != Size || M.Nrows() != Vdim) {
    cerr << "***ERROR: GenSetBase::generateAll(Matrix,...) "
	 << "dimesion of M expected to be "
	 << Vdim << "-by-" << Size 
	 << " but is " << M.Nrows() << "-by-" << M.Ncols()
	 << endl;
    return false;
  }
  ColumnVector xi(Vdim);
  for (int i=1; i<=Size; i++) {
    generate(i, D, X, xi);
    M.Column(i) = xi;
  }
  return true;
}