📄 btl_matrix_algorithms.h
字号:
//
// btl_matrix_algorithms.h
//
// This file contains the generic numeric functions for manipulation of matrices
//
// These classes are part of the Bioinformatics Template Library (BTL).
//
// Copyright (C) 1997, 1998 Birkbeck College, Malet Street, London, U.K.
// Copyright (C) 2004, 2005 University College, Gower Street, London, U.K.
//
// This library is free software; you can redistribute it and/or modify it
// under the terms of the GNU Library General Public License as published
// by the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version. This library is distributed in the
// hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
// implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the GNU Library General Public License for more details.
// You should have received a copy of the GNU Library General Public
// License along with this library; if not, write to the Free Software
// Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
///////////////////////////////////////////////////////////////////////////
#if !defined (BTL_MATRIXALGORITHMS_H)
#define BTL_MATRIXALGORITHMS_H 1
/**#: [Description ="A collection of generic numerical algorithms for
matrix algebra."]
[Summary = "generic matrix algorithms"]
[Authors = "D.S.Moss, W.R.Pitt, I.Tickle, M.A.Williams"]
[Files = "<A HREF=./btl/btl_matrix_algorithms.h>btl_matrix_algorithms.h</A>"]
[Dependencies="<A HREF=#numeric_vector>btl_numeric_vector.h</A>,
<A HREF=#matrix>btl_matrix.h</A>"]
[Prerequisites="None<P>"
<P>
*/
#include <vector>
#include <iterator>
#include "BTL.h"
#include "btl_numeric_vector.h"
#include "btl_matrix.h"
_BTL_BEGIN_NAMESPACE
using namespace std;
//..........................................................................................
// MATRIX ALGEBRA
//...............................................................................................
/**#: [Description="Matrix multiplication.
The output matrix is nrows1*ncols2 in size."]
[Restrictions="ncols1 must equal nrows2."] */
template <class InputIterator1, class InputIterator2, class OutputIterator>
OutputIterator
matrix_product(InputIterator1 first1, InputIterator1 last1, BTL_INT nrows1,
InputIterator2 first2, InputIterator2 last2, BTL_INT nrows2,
OutputIterator result)
{
BTL_INT ncols1 = (last1-first1)/nrows1;
BTL_INT ncols2 = (last2-first2)/nrows2;
#if defined(BTL_DEBUG_VERSION)
if (ncols1 != nrows2)
FATAL_ERROR("The number of rows in both matrices must be equal.");
#endif
BTL_INT i,j;
numeric_vector<typename iterator_traits<OutputIterator>::value_type> temp(nrows2);
typename iterator_traits<OutputIterator>::value_type tempresult = 0;
for (j=0;j<ncols2;j++)
{
copy_column(first2,last2,nrows2,temp.begin(),(j+1));
for (i=0;i<nrows1;i++)
{
tempresult=0;
*(result+((i*nrows1)+j)) +=
scalar_product((first1+(ncols1*i)),(first1+(ncols1*(i+1))),temp.begin(),tempresult);
}
}
return (result+(nrows1*ncols2));
}
//...............................................................................................
/**#: [Description="Matrix multiplication: Transpose(M1)*M2.
The output matrix is ncols1*ncols2 in size."]
[Restrictions="Both matrices must have the same number of rows"] */
template <class InputIterator1, class InputIterator2, class OutputIterator>
OutputIterator
matrixtranspose_matrix_product(InputIterator1 first1, InputIterator1 last1, BTL_INT nrows1,
InputIterator2 first2, InputIterator2 last2, BTL_INT nrows2,
OutputIterator result)
{
BTL_INT ncols1 = (last1-first1)/nrows1;
BTL_INT ncols2 = (last2-first2)/nrows2;
#if defined(BTL_DEBUG_VERSION)
if (nrows != nrows2)
FATAL_ERROR("The number of rows in both matrices must be equal.");
#endif
InputIterator2 j;
BTL_INT col1, col2, row = 1;
for(col1=1;first1<last1;first1++,col1++)
{
if (col1 > ncols1) col1 = 1, row = row + 1;
for (j=(first2+((row-1)*ncols2)), col2=1; j<(first2+(row*ncols2)); j++, col2++)
*(result+((col1-1)*ncols1)+(col2-1)) += *first1 * *j;
}
return (result+(ncols1*ncols2));
}
//...............................................................................................
/**#: [Description="Matrix multiplication: M1*Transpose(M2)
The output matrix is nrows1*nrows2 in size."]
[Restrictions="Both matrices must have the same number of columns"] */
template <class InputIterator1, class InputIterator2, class OutputIterator>
OutputIterator
matrix_matrixtranspose_product(InputIterator1 first1, InputIterator1 last1, BTL_INT nrows1,
InputIterator2 first2, InputIterator2 last2, BTL_INT nrows2,
OutputIterator result)
{
BTL_INT ncols = (last1-first1)/nrows1;
BTL_INT ncols2 = (last2-first2)/nrows2;
typename iterator_traits<OutputIterator>::value_type tempresult=0;
#if defined(BTL_DEBUG_VERSION)
if ( ncols != ncols2)
FATAL_ERROR("The number of columns in both matrices must be equal.");
#endif
BTL_INT i,j;
for(i=0;i < nrows1; i++)
{
for(j=0;j < nrows2; j++)
{
tempresult = 0;
*(result+((i*nrows2)+j))
+= scalar_product((first1+(i*ncols)),(first1+((i+1)*ncols)),
(first2+(j*ncols)),tempresult);
}
}
return (result+(nrows1*nrows2));
}
//..........................................................................................
/**#: [Description="Matrix transpose. Can be used in-place."] */
template <class InputIterator, class OutputIterator>
OutputIterator transpose(InputIterator first, InputIterator last, BTL_INT nrows,
OutputIterator result)
{
BTL_INT ncols= (last-first)/nrows;
BTL_INT i,j;
typename iterator_traits<InputIterator>::value_type temp;
if (nrows == ncols)
{
for (i=0;i<nrows;i++)
{
for (j=i;j<ncols;j++)
{
temp = *(first+(ncols*i)+j);
*(result+(ncols*i)+j) = *(first+i+(ncols*j));
*(result+i+(ncols*j)) = temp;
}
}
}
else
{
matrix<typename iterator_traits<InputIterator>::value_type> temp_matrix(first,last,nrows,ncols);
for (i=1;i<=nrows;i++)
{
for (j=1;j<=ncols;j++)
{
*(result+(nrows*(j-1))+(i-1)) = temp_matrix(i,j);
}
}
}
return (result+(nrows*ncols));
}
// AS ALL SWAPS OCCUR WITHIN PAIRS OF SYMMETRICALLY ARRANGED DIAGONALS THERE IS CERTAINLY A
// MORE MEMORY EFFICIENT WAY OF DOING THIS TRANSPOSE IN THE NON-SQUARE CASE
//..........................................................................................
/**#: [Description="Calculate the mean of each column of a
matrix and store the answers in a container of size >= ncols."] */
template<class InputIterator, class OutputIterator>
OutputIterator column_means(InputIterator first, InputIterator last,
BTL_INT nrows, OutputIterator result)
{
BTL_INT ncols = (last-first)/nrows;
BTL_REAL reciprocal_nrows = 1.0/nrows;
OutputIterator begin = result;
OutputIterator end = result+ncols;
while (first != last)
{
for ( result=begin; result!=end; result++, first++)
*result += *first;
}
for ( result=begin; result!=end; result++, first++)
*result *= reciprocal_nrows;
return ++result;
}
//...........................................................................................
/**#: [Description="Copy a column of this matrix into another container.
Columns are numbered from 1. in the FORTRAN manner."] */
template<class InputIterator, class OutputIterator>
OutputIterator copy_column(InputIterator first, InputIterator last,
BTL_INT nrows, OutputIterator result, BTL_INT column_index)
{
BTL_INT ncols = (last-first)/nrows;
if (column_index > ncols || column_index < 1)
FATAL_ERROR("matrix index out of range");
for ( first+=(column_index-1) ; (first - last) < 0; result++, first+=ncols)
*result = *first;
return ++result;
}
//..........................................................................................
/**#: [Description="Calculates matrix determinant"]
[Restrictions="Square matrices only."] */
template<class InputIterator, class T>
T determinant(InputIterator first, InputIterator last, BTL_INT nrows, T det)
{
if ( (nrows*nrows) != (last-first))
FATAL_ERROR("Cannot calculate determinant of non-square matrix");
// For matrices with up to 4 rows use explicit Laplace expansion
if (nrows == 1) det += *first;
if (nrows == 2)
{
det += (*first * *(first+3)) - (*(first+1) * *(first+2));
}
if (nrows == 3)
{
det += *first * ((*(first+4) * *(first+8)) - (*(first+7) * *(first+5)))
- *(first+3) * ((*(first+1) * *(first+8)) - (*(first+7) * *(first+2)))
+ *(first+6) * ((*(first+1) * *(first+5)) - (*(first+4) * *(first+2)));
}
if (nrows == 4)
{
det += ( ((*first * *(first+5)) - (*(first+4) * *(first+1)))
* ((*(first+10) * *(first+15)) - (*(first+14) * *(first+11))) )
+ ( ((*(first+1) * *(first+8)) - (*first * *(first+9)))
* ((*(first+6) * *(first+15)) - (*(first+14) * *(first+7))) )
+ ( ((*first * *(first+13)) - (*(first+1) * *(first+12)))
* ((*(first+6) * *(first+11)) - (*(first+10) * *(first+7))) )
+ ( ((*(first+4) * *(first+9)) - (*(first+5) * *(first+8)))
* ((*(first+2) * *(first+15)) - (*(first+14) * *(first+3))) )
+ ( ((*(first+5) * *(first+12)) - (*(first+4) * *(first+13)))
* ((*(first+2) * *(first+11)) - (*(first+10) * *(first+3))) )
+ ( ((*(first+8) * *(first+13)) - (*(first+9) * *(first+12)))
* ((*(first+2) * *(first+7)) - (*(first+6) * *(first+3))) );
}
// For matrices with greater than 4 rows use iterative Laplace expansion
if (nrows > 4)
{
int sign;
matrix<T> Minor((nrows-1)); // creates appropriate amount of space
typename matrix<T>::iterator end, begin = Minor.begin();
for (BTL_INT i = 1; i <= nrows; i++)
{
end = matrix_minor(first,last,nrows,begin,i,1);
sign = (2*(i%2)) - 1;
det += sign * *(first+(nrows*(i-1)))
* determinant(begin,end,(nrows-1),T(0));
}
}
return det;
}
// CAN POSSIBLY BE MORE EFFICIENTLY WRITTEN USING TEMPLATE METAPROGRAMS.
// ALSO AS THE LAPLACE EXPANSION GROWS AS N! - FOR MATRICES WITH MORE THAN 6 OR 7 ROWS
// WE SHOULD USE LU DECOMPOSITION.
//..........................................................................................
/**#: [Description="Find the minor of a matrix.
The minor denoted by minormatrix(.,.,,.,i,j) is found by
eliminating the row i and column j from the parent matrix."] */
template<class InputIterator, class OutputIterator>
OutputIterator matrix_minor(InputIterator first, InputIterator last, BTL_INT nrows,
OutputIterator result, BTL_INT i, BTL_INT j)
{
BTL_INT ncols= (last-first)/nrows;
if ( i<1 || i>nrows || j<1 || j>ncols)
FATAL_ERROR("One or both indexes are out of range");
BTL_INT row,col;
BTL_INT nless = (ncols-1);
i-=1;
j-=1;
for (row = 0; row < nrows; row++)
for (col = 0; col < ncols; col++)
{
if (row < i)
{
if (col < j)
*(result + (nless*row + col))
= *(first + (ncols*row + col));
else if (col > j)
*(result + (nless*row + col - 1))
= *(first + (ncols*row + col));
}
else if (row > i)
{
if (col < j)
*(result + (nless*(row-1) + col))
= *(first + (ncols*row + col));
else if (col > j)
*(result + (nless*(row-1) + col - 1))
= *(first + (ncols*row + col));
}
}
return (result+((nrows-1)*(ncols-1)));
}
//..........................................................................................
/**#: [Description="Matrix adjoint"]
[Restrictions="Square matrices only"].*/
template<class InputIterator, class OutputIterator>
OutputIterator adjoint(InputIterator first, InputIterator last,
BTL_INT nrows, OutputIterator result)
{
if ((nrows*nrows) != (last-first))
FATAL_ERROR("Cannot calculate adjoint of non-square matrix");
if (nrows < 3)
FATAL_ERROR("Cannot calculate adjoint of matrix with less than 3 rows");
matrix<BTL_REAL> Minor((nrows-1));
typename matrix<BTL_REAL>::iterator end, begin = Minor.begin();
BTL_INT col,row;
int sign; // must declare as signed integer can't rely on automatic type conversion
for (col = 1; col <= nrows; col++)
for (row = 1; row <= nrows; row++)
{
end=matrix_minor(first,last,nrows,begin,row,col);
sign = (1-(2*((row+col)%2)));
*(result+(nrows*(col-1)+row-1))
= sign * determinant(begin,end,(nrows-1),BTL_REAL(0)) ;
}
return (result+(nrows*nrows));
}
//..........................................................................................
/**#: [Description="Calculate inverse of a small square matrix
(inverse = adjoint/determinant). The result matrix should
be of a floating point type."] */
template<class InputIterator, class OutputIterator>
OutputIterator inverse_square(InputIterator first, InputIterator last,
BTL_INT nrows, OutputIterator result)
{
// Work out determinant
BTL_REAL det = BTL_REAL(0);
det = determinant(first,last,nrows,det);
BTL_REAL reciprocal_det = 1.0/det;
BTL_REAL small = numeric_limits<BTL_REAL>::epsilon();
if ( det*det < small*small)
FATAL_ERROR("Matrix singular to machine precision: try SV_decomposition ");
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -