btl_matrix.h

利用这个模板可以分析基因表达数据

	numeric_vector<T>
	operator*(const numeric_vector<T> &v) const
	{

	#if defined(BTL_DEBUG_VERSION)
	    if (ncols != v.size()) 
	    	FATAL_ERROR("The size of the numeric_vector<T> must equal the number of "
    	            "columns in the Matrix");
	#endif

	    numeric_vector<T> result(nrows);
	    const_iterator i;
	    const_iterator j = mat.begin();
	    for (iterator k = result.begin(); k < result.end(); k++)
	        for (i = v.begin(); i < v.end(); i++, j++)
	            *k += *i * *j; 
	    return result;
	}

//.............................................................................
	    	/**#: [Description="Multiple each element by a number."] */
	matrix<T>&
	operator*=(const value_type& v)
	{
	    for (iterator i=mat.begin(); i!=mat.end(); i++)
		*i *= v;
	    return *this;
	}


//.............................................................................
	    	/**#: [Description="Multiple each element by a number"] */
	matrix<T>
	operator*(const value_type& v) const
	{
	    matrix<T> result = *this;
	    for (iterator i=result.begin(); i!=result.end(); i++)
		*i *= v;
	    return result;
	}

//.............................................................................
	    	/**#: [Description="Divide each element by a number"] */
	matrix<T>&
	operator/=(const value_type& v)
	{
	#if defined(BTL_DEBUG_VERSION)
	    if (v == 0.0) WARNING("Divide by zero");
	#endif

	    BTL_REAL y_reciprocal = 1.0 / v;
	    for (iterator i=mat.begin(); i!=mat.end(); i++)
	    	*i *= y_reciprocal;
	    return *this;
	}


//.............................................................................
	    	/**#: [Description="Divide each element by a number"] */
	matrix<T>
	operator/(const value_type& v) const
	{
	#if defined(BTL_DEBUG_VERSION)
	    if (v == 0.0) WARNING("Divide by zero");
	#endif

	    matrix<T> result= *this;
	    BTL_REAL y_reciprocal = 1.0 / v;
	    for (iterator i=result.begin(); i!=result.end(); i++)
	    	*i *= y_reciprocal;
	    return result;
	}


//.............................................................................
	    	/**#: [Description="Subtraction of a number from each element."]
	    	      */
    	matrix<T>&
    	operator-=(const value_type &v)
	{
	    for (iterator i=mat.begin(); i!=mat.end(); i++)
	    	*i -= v;
	    return *this;
	}    	

//.............................................................................
	    	/**#: [Description="Subtraction of a number from each element."]
	    	      */
	matrix<T>
	operator-(const value_type &v) const
	{
	    matrix<T> result = *this;
	    for (iterator i=result.begin(); i!=result.end(); i++)
	    	*i -= v;
	    return result;
	}    	


//.............................................................................
	    	/**#: [Description="Addition of a number to each element."] */
    	matrix<T>&
    	operator+=(const value_type &v)
	{
	    for (iterator i=mat.begin(); i!=mat.end(); i++)
	    	*i += v;
	    return *this;
	}    	


//.............................................................................
	    	/**#: [Description="Addition of a number to each element."] */
	matrix<T>
	operator+(const value_type &v) const
	{
	    matrix<T> result = *this;
	    for (iterator i=result.begin(); i!=result.end(); i++)
	    	*i += v;
	    return result;
	}    	


//.............................................................................
	    	/**#: [Description="Subtraction of a vector from the rows
	    	       of a matrix"]
	    	      [Restrictions="The the size of numeric_vector v must equal the
	    	       number of columns in this matrix."] */
    	matrix<T>&
    	operator-=(const numeric_vector<T> &v)
	{

	#if defined(BTL_DEBUG_VERSION)
	    if (v.size() != ncols)
	    	FATAL_ERROR("The the size of the input numeric_vector<T> must equal the number of "
	    	            "columns in the Matrix");
	#endif

	    iterator i = mat.begin();
	    const_iterator j;
	    for (int k=0; k<nrows; k++)
	    	for (j=v.begin(); j!=v.end(); j++, i++)
	    	     *i -= *j;
    	    
	    return *this;
	}


//.............................................................................
	    	/**#: [Description="Subtraction of a vector from the rows
	    	       of a matrix"]
	    	      [Restrictions="The the size of numeric_vector v must equal the
	    	       number of columns in this matrix."] */
	matrix<T>
	operator-(const numeric_vector<T> &v) const
	{
	    matrix<T> result = *this;
	    result -= v;
	    return result;
	}


//.............................................................................
	    	/**#: [Description="Addition of a vector to the rows
	    	       of a matrix"]
	    	      [Restrictions="The the size of numeric_vector v must equal the
	    	       number of columns in this matrix."] */
    	matrix<T>&
    	operator+=(const numeric_vector<T> &v)
	{

	#if defined(BTL_DEBUG_VERSION)
	    if (v.size() != ncols)
	    	FATAL_ERROR("The the size of the input numeric_vector<T> must equal the number of "
	    	            "columns in this Matrix");
	#endif

	    iterator i = mat.begin();
	    const_iterator j;
	    for (int k=0; k<nrows; k++)
	    	for (j=v.begin(); j!=v.end(); j++, i++)
	    	     *i += *j;
	    return *this;
	}


//.............................................................................
	    	/**#: [Description="Addition of a vector to the rows
	    	       of a matrix"]
	    	      [Restrictions="The the size of numeric_vector v must equal the
	    	       number of columns in this matrix."] */
	matrix<T>
	operator+(const numeric_vector<T> &v) const
	{
	    matrix<T> result = *this;
	    result += v;
	    return result;
	}


//.............................................................................
	    	/**#: [Description="matrix subtraction"]
	    	      [Restrictions="The the size of matrix m must equal the
	    	       size of this matrix."] */
    	matrix<T>&
    	operator-=(const matrix<T> &m)
	{
	#if defined(BTL_DEBUG_VERSION)
	    if (nrows != m.num_rows() || ncols != m.num_cols())
	    	FATAL_ERROR("Both matrices must be the same size");
	#endif

	    iterator i;
	    const_iterator j;
	    for (i=mat.begin(), j=m.begin(); i!=mat.end(); i++, j++)
	    	*i -= *j;
	    return *this;
	}


//.............................................................................
	    	/**#: [Description="Matrix subtraction"]
	    	      [Restrictions="The the size of Matrix m must equal the
	    	       size of this Matrix."] */
	matrix<T>
	operator-(const matrix<T> &m) const
	{

	#if defined(BTL_DEBUG_VERSION)
	    if (nrows != m.num_rows() || ncols != m.num_cols())
	    	FATAL_ERROR("Both matrices must be the same size");
	#endif

	    matrix<T>  result = *this;
	    const_iterator i;
	    iterator j;
	    for (i=m.begin(), j=result.begin(); i!=m.end(); i++, j++)
	    	*j -= *i;
	    return result;
	}


//.............................................................................
	    	/**#: [Description="Matrix addition"]
	    	      [Restrictions="The the size of Matrix m must equal the
	    	       size of this Matrix."] */
    	matrix<T>&
    	operator+=(const matrix<T> &m)
	{

	#if defined(BTL_DEBUG_VERSION)
	    if (nrows != m.num_rows() || ncols != m.num_cols())
	   	FATAL_ERROR("Both matrices must be the same size");
	#endif

	    iterator i;
	    const_iterator j;
	    for (i=mat.begin(), j=m.begin(); i!=mat.end(); i++, j++)
	    	*i += *j;
	    return *this;
	}


//.............................................................................
	    	/**#: [Description="Matrix addition"]
	    	      [Restrictions="The the size of Matrix m must equal the
	    	       size of this Matrix."] */
	matrix<T>
	operator+(const matrix<T> &m) const
	{

	#if defined(BTL_DEBUG_VERSION)
	    if (nrows != m.num_rows() || ncols != m.num_cols())
	    	FATAL_ERROR("Both matrices must be the same size");
	#endif

	    matrix<T>  result = *this;
	    const_iterator i;
	    iterator j;
	    for (i=m.begin(), j=result.begin(); i!=m.end(); i++, j++)
	    	*j += *i;
	    return result;
	}


//.............................................................................
	    	/**#: [Description="Matrix multiplication Transpose(m1) * m2
	    	       N.B. this is less efficient than (but not equivalent to)
	    	       operator%"]
	    	      [Restrictions="Both *this and the input Matrix must have
	    	       the same number of rows"] */
	matrix<T>
	operator&(const matrix<T>& m) const
	{

	#if defined(BTL_DEBUG_VERSION)
	    if (m.num_rows() != nrows)
	    	FATAL_ERROR("The input Matrix must have the same number of rows as this Matrix");
	#endif
    
	    matrix<T> result(ncols,m.num_cols());
	    numeric_vector<T> temp(nrows);
	    size_type i,j,k,iarr;
	    for (i=0, iarr=0; i<ncols; i++, iarr=result.ncols*i)
	    {
	    	for (j=0; j<nrows; j++) temp[j] = mat[ncols*j + i];
	    	for (j=0; j<m.ncols; j++)
	    	    for (k=0; k<nrows; k++)    	    
	    	    	result.mat[iarr + j] += temp[k] * m.mat[m.ncols*k + j]; // m1(i,j)+= *this(k,i) * m(k,j)
	    }   	
	    return result;
	}

//.............................................................................
	    	/**#: [Description="Matrix multiplication m1 * Transpose(m2)
	    	       N.B. this is more efficient than (but not equivalent to)
	    	       operator&"]
	    	      [Restrictions="Both *this and the input Matrix must have
	    	       the same number of columns"] */
	matrix<T> 
	operator%(const matrix<T>& m) const
	{

	#if defined(BTL_DEBUG_VERSION)
	    if (m.num_cols() != ncols)
	    	FATAL_ERROR("Both this and the input Matrix must have the same number "
	    	            "of columns");
	#endif
    
	    matrix<T> result(nrows,m.num_rows());
	    for (size_type i=1; i<=nrows; i++)
	    	for (size_type j=1; j<=m.nrows; j++)
	    	    for (size_type k=1; k<=ncols; k++)
	    	    	result(i,j) += (*this)(i,k) * m(j,k);

	    return result;
	}


//.............................................................................
    	    	/**#: [Description="Equality operator"]*/
	bool 
	operator==(const matrix<T>& m) const
	{ 
	    if ( nrows != m.num_rows() || ncols != m.num_cols() )
	    	return false;
    	
	    iterator i;
	    const_iterator j;
	    for (i=mat.begin(), j=m.begin(); i!=mat.end(); i++, j++)
            {
	    	if (*i != *j) return false;
            }
            return true;
	}

//.............................................................................
// I/O FUNCTIONS
//.............................................................................
// ostream operator

	friend ostream&
	operator<<(ostream &os, const matrix<T> &m)
	{
	        typedef typename matrix<T>::size_type size_type;

		os.setf(ios::showpoint);
		os.setf(ios::fixed, ios::floatfield);

		for ( size_type i=1; i<=m.nrows; i++)
		{
			os << '\n';
			for (size_type j=1; j<=m.ncols; j++)
				os << m(i,j) << " ";
		}
		os << '\n';
		return os;
	}

//.............................................................................
// istream operator

    	friend istream&
	operator>>(istream &is, matrix<T> &m)
	{
    
	    typedef typename matrix<T>::size_type size_type;
	    typedef typename matrix<T>::value_type value_type;

	    size_type M, N;

	    is >> M >> N;

	    if (M<1 || N<1) 
	    {
	    	WARNING("Both dimensions must be greater than 0");
		return is;
	    }

//	    if ( !(M == m.nrows && N == m.ncols) )
//	    {
//	        ~m();
//		matrix<T> m(M,N);
//	    	m.nrows = M;
//		m.ncols = N;
//	    }

	    for (size_type i=0; i<M; i++)
	        for (size_type j=0; j<N; j++)
	        {
	            is >>  m(i,j);
	        }

	    return is;
	}

//.............................................................................
// OTHER FUNCTIONS.
//.............................................................................

    	
}; // class matrix

_BTL_END_NAMESPACE

#endif