LOOS  v2.3.2
MatrixOps.hpp
1 /*
2  MatrixOps.hpp
3 */
4 
5 /*
6  This file is part of LOOS.
7 
8  LOOS (Lightweight Object-Oriented Structure library)
9  Copyright (c) 2009, Tod D. Romo, Alan Grossfield
10  Department of Biochemistry and Biophysics
11  School of Medicine & Dentistry, University of Rochester
12 
13  This package (LOOS) is free software: you can redistribute it and/or modify
14  it under the terms of the GNU General Public License as published by
15  the Free Software Foundation under version 3 of the License.
16 
17  This package is distributed in the hope that it will be useful,
18  but WITHOUT ANY WARRANTY; without even the implied warranty of
19  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  GNU General Public License for more details.
21 
22  You should have received a copy of the GNU General Public License
23  along with this program. If not, see <http://www.gnu.org/licenses/>.
24 */
25 
26 #if !defined LOOS_MATRIXOPS_HPP
27 #define LOOS_MATRIXOPS_HPP
28 
29 #include <loos_defs.hpp>
30 #include <MatrixImpl.hpp>
31 #include <exceptions.hpp>
32 #include <sorting.hpp>
33 #include <utils_random.hpp>
34 #include <TimeSeries.hpp>
35 #include <stdexcept>
36 #include <cmath>
37 
38 
39 
40 
41 
42 namespace loos {
43 
44 
45  typedef Math::Matrix<float, Math::ColMajor> RealMatrix;
46  typedef Math::Matrix<double, Math::ColMajor> DoubleMatrix;
47 
55  namespace Math {
56 
57 
59  boost::tuple<RealMatrix, RealMatrix, RealMatrix> svd(RealMatrix& M);
60 
62 
65  boost::tuple<DoubleMatrix, DoubleMatrix, DoubleMatrix> svd(DoubleMatrix& M);
66 
68 
73  DoubleMatrix eigenDecomp(DoubleMatrix& M);
74 
75  RealMatrix eigenDecomp(RealMatrix& M);
76 
78  RealMatrix MMMultiply(const RealMatrix& A, const RealMatrix& B, const bool transa = false, const bool transb = false);
79  DoubleMatrix MMMultiply(const DoubleMatrix& A, const DoubleMatrix& B, const bool transa = false, const bool transb = false);
80 
82  RealMatrix invert(RealMatrix& A, const float eps = 1e-5);
83  DoubleMatrix invert(DoubleMatrix& A, const double eps = 1e-5);
84 
86  template<typename T>
87  T eye(const uint n) {
88  T I(n, n);
89  for (uint i=0; i<n; ++i)
90  I(i, i) = 1.0;
91  return(I);
92  }
93 
94  // ***DEPRECATED***
95  DoubleMatrix deye(const uint n);
96 
97 
99  template<typename T>
100  T submatrix(const T& M, const Range& rows, const Range& cols) {
101  uint m = rows.second - rows.first;
102  uint n = cols.second - cols.first;
103 
104  T A(m,n);
105  for (uint i=0; i < n; ++i)
106  for (uint j=0; j < m; ++j)
107  A(j,i) = M(j+rows.first, i+cols.first);
108 
109  return(A);
110  }
111 
114 
116  void operator+=(RealMatrix& A, const RealMatrix& B);
117  RealMatrix operator+(const RealMatrix& A, const RealMatrix& B);
118  void operator-=(RealMatrix& A, const RealMatrix& B);
119  RealMatrix operator-(const RealMatrix& A, const RealMatrix& B);
120  void operator*=(RealMatrix& A, const float d);
121  RealMatrix operator*(const RealMatrix& A, const float d);
122  void operator*=(RealMatrix& A, const RealMatrix& B);
123  RealMatrix operator*(const RealMatrix& A, const RealMatrix& B);
124  RealMatrix operator-(RealMatrix& A);
125 
126  void operator+=(DoubleMatrix& A, const DoubleMatrix& B);
127  DoubleMatrix operator+(const DoubleMatrix& A, const DoubleMatrix& B);
128  void operator-=(DoubleMatrix& A, const DoubleMatrix& B);
129  DoubleMatrix operator-(const DoubleMatrix& A, const DoubleMatrix& B);
130  void operator*=(DoubleMatrix& A, const double d);
131  DoubleMatrix operator*(const DoubleMatrix& A, const double d);
132  void operator*=(DoubleMatrix& A, const DoubleMatrix& B);
133  DoubleMatrix operator*(const DoubleMatrix& A, const DoubleMatrix& B);
134  DoubleMatrix operator-(DoubleMatrix& A);
135 
137  template<typename T>
138  T permuteColumns(const T& A, const std::vector<uint>& indices) {
139  if (indices.size() != A.cols())
140  throw(std::logic_error("indices to permuteColumns must match the size of the matrix"));
141 
142  T B(A.rows(), A.cols());
143 
144  for (uint i=0; i<A.cols(); ++i) {
145  if (indices[i] > A.cols())
146  throw(std::out_of_range("Permutation index is out of bounds"));
147  for (uint j=0; j<A.rows(); ++j)
148  B(j, i) = A(j, indices[i]);
149  }
150 
151  return(B);
152  }
153 
155  template<typename T>
156  T permuteRows(const T& A, const std::vector<uint>& indices) {
157  if (indices.size() != A.rows())
158  throw(std::logic_error("indices to permuteRows must match the size of the matrix"));
159 
160  T B(A.rows(), A.cols());
161 
162  for (uint j=0; j<A.rows(); ++j) {
163  if (indices[j] > A.rows())
164  throw(std::out_of_range("Permutation index is out of bounds"));
165  for (uint i=0; i<A.cols(); ++i)
166  B(j, i) = A(indices[j], i);
167  }
168 
169  return(B);
170  }
171 
173  template<typename T>
174  T transpose(const T& A) {
175  T B(A.cols(), A.rows());
176 
177  for (uint j=0; j<A.rows(); ++j)
178  for (uint i=0; i<A.cols(); ++i)
179  B(i, j) = A(j, i);
180 
181  return(B);
182  }
183 
184 
186  template<typename T>
187  T shuffleColumns(const T& A) {
188  std::vector<float> random_numbers(A.cols());
189  base_generator_type& rng = rng_singleton();
190  boost::uniform_real<> rngmap(0.0, 1.0);
191  boost::variate_generator<base_generator_type&, boost::uniform_real<> > rnd(rng, rngmap);
192 
193  for (uint i=0; i<A.cols(); ++i)
194  random_numbers[i] = rnd();
195 
196  std::vector<uint> indices = sortedIndex(random_numbers);
197  return(permuteColumns(A, indices));
198  }
199 
200 
202  template<typename T>
203  T shuffleColumnVector(const T& v) {
204  std::vector<float> random_numbers(v.size());
205  base_generator_type& rng = rng_singleton();
206  boost::uniform_real<> rngmap(0.0, 1.0);
207  boost::variate_generator<base_generator_type&, boost::uniform_real<> > rnd(rng, rngmap);
208 
209  for (uint i=0; i<v.size(); ++i)
210  random_numbers[i] = rnd();
211 
212  std::vector<uint> indices = sortedIndex(random_numbers);
213 
214  // This forces the result to -always- be a column vector...
215  T u(v.size(), 1);
216  for (uint i=0; i<v.size(); ++i)
217  u[i] = v[indices[i]];
218 
219  return(u);
220  }
221 
222 
223  template<typename T>
224  void reverseColumns(T& A) {
225  uint m = A.rows();
226  uint n = A.cols();
227  uint k = n / 2;
228 
229  for (uint i=0; i<k; ++i)
230  for (uint j=0; j<m; ++j)
231  std::swap(A(j,i), A(j,n-i-1));
232  }
233 
234  template<typename T>
235  void reverseRows(T& A) {
236  uint m = A.rows();
237  uint n = A.cols();
238  uint k = m / 2;
239 
240  for (uint j=0; j<k; ++j)
241  for (uint i=0; i<n; ++i)
242  std::swap(A(j, i), A(m-j-1, i));
243  }
244 
245 
246 
247  namespace {
248  template<typename T>
249  double colDotProd(const T& A, const uint i, const T& B, const uint j) {
250  double sum = 0.0;
251  for (uint k=0; k<A.rows(); ++k)
252  sum += A(k,i) * B(k,j);
253  return(sum);
254  }
255  }
256 
257  template<typename T>
258  double subspaceOverlap(const T& A, const T& B, uint nmodes = 0) {
259  if (A.rows() != B.rows())
260  throw(NumericalError("subspaceOverlap: Matrices have different dimensions"));
261 
262  if (nmodes == 0)
263  nmodes = A.cols();
264  if (nmodes > A.cols() || nmodes > B.cols())
265  throw(NumericalError("Requested number of modes exceeds matrix dimensions"));
266 
267  double sum = 0.0;
268  for (uint i=0; i<nmodes; ++i)
269  for (uint j=0; j<nmodes; ++j) {
270  double d = colDotProd(A, i, B, j);
271  sum += d*d;
272  }
273 
274  return(sum / nmodes);
275  }
276 
277 
278 
280 
302  template<typename T>
303  double covarianceOverlap(const T& lamA, const T& UA, const T& lamB, const T& UB) {
304  if (!(UA.rows() == UB.rows() && lamA.rows() <= UA.cols() && lamB.rows() <= UB.cols()))
305  throw(NumericalError("covarianceOverlap: Matrices have incorrect dimensions"));
306 
307  // X = abs(UB'*UA)
308  T X = MMMultiply(UB,UA,true,false);
309  for (ulong i = 0; i < X.size(); ++i)
310  X[i] = fabs(X[i]);
311 
312  // L = abs(lamB*lamA')
313  T L = MMMultiply(lamB, lamA, false, true);
314 
315  // y = sum(sum(sqrt(L).*X.*X));
316  double y = 0;
317  for (ulong i = 0; i<X.size(); ++i)
318  y += sqrt(L[i]) * X[i] * X[i];
319 
320  // e = sum(s.*t);
321  double e =0;
322  for (ulong i=0; i<lamA.size(); ++i)
323  e += lamA[i] + lamB[i];
324 
325  double num = e - 2.0 * y;
326  double co = 1.0 - sqrt( fabs(num) / e );
327 
328  return(co);
329  }
330 
331 
332  // Returns: z-score, raw covariance overlap, and stddev used in the z-score
333  template<typename T>
334  boost::tuple<double, double, double> zCovarianceOverlap(const T& lamA, const T& UA, const T& lamB, const T& UB, const uint tries) {
335  double coverlap = covarianceOverlap(lamA, UA, lamB, UB);
336  std::vector<double> random_coverlaps(tries);
337 
338  for (uint i=0; i<tries; ++i) {
339  T shuffled_lamA = shuffleColumnVector(lamA);
340  T shuffled_lamB = shuffleColumnVector(lamB);
341  random_coverlaps[i] = covarianceOverlap(shuffled_lamA, UA, shuffled_lamB, UB);
342  }
343 
344  TimeSeries<double> ts(random_coverlaps);
345  double score = (coverlap - ts.average()) / ts.stdev();
346 
347  boost::tuple<double, double, double> result(score, coverlap, ts.stdev());
348  return(result);
349  }
350 
351 
352  };
353 
354 
355 };
356 
357 
358 #endif
T transpose(const T &A)
Transposes a matrix.
Definition: MatrixOps.hpp:174
Simple matrix template class using policy classes to determine behavior.
Definition: MatrixImpl.hpp:53
boost::tuple< RealMatrix, RealMatrix, RealMatrix > svd(RealMatrix &M)
Compute the SVD of a single precision matrix.
Definition: MatrixOps.cpp:94
DoubleMatrix eigenDecomp(DoubleMatrix &M)
Compute eigendecomposition of M.
Definition: MatrixOps.cpp:36
T permuteRows(const T &A, const std::vector< uint > &indices)
Returns a copy of the matrix with the rows permuted by the indices.
Definition: MatrixOps.hpp:156
double covarianceOverlap(const T &lamA, const T &UA, const T &lamB, const T &UB)
Computes the covariance overlap between two subspaces.
Definition: MatrixOps.hpp:303
T shuffleColumns(const T &A)
Randomly shuffle the columns of a matrix.
Definition: MatrixOps.hpp:187
RealMatrix invert(RealMatrix &A, const float eps)
Pseudo-inverse of a matrix using the SVD.
Definition: MatrixOps.cpp:230
void operator+=(RealMatrix &A, const RealMatrix &B)
Overloaded operators for RealMatrix matrices (see important note below)
Definition: MatrixOps.cpp:286
T shuffleColumnVector(const T &v)
! Randomly shuffle the rows of a single column vector
Definition: MatrixOps.hpp:203
def subspaceOverlap(A, B, nmodes=0)
Definition: subspace.py:10
Coord< T > operator*(const Matrix44< T > &, const Coord< T > &)
Definition: Matrix44.hpp:244
std::pair< uint, uint > Range
Specify a range for columns/rows [first,second)
Definition: MatrixImpl.hpp:49
base_generator_type & rng_singleton(void)
Suite-wide random number generator singleton.
Exception caused by a blas/atlas error.
Definition: exceptions.hpp:91
T submatrix(const T &M, const Range &rows, const Range &cols)
Extracts a submatrix.
Definition: MatrixOps.hpp:100
T eye(const uint n)
An identity matrix of size n.
Definition: MatrixOps.hpp:87
T permuteColumns(const T &A, const std::vector< uint > &indices)
Returns a copy of the matrix with the columns permuted by the indices.
Definition: MatrixOps.hpp:138
std::vector< uint > sortedIndex(const T &A)
Definition: sorting.hpp:71
Namespace for most things not already encapsulated within a class.
RealMatrix MMMultiply(const RealMatrix &A, const RealMatrix &B, const bool transa, const bool transb)
Matrix-matrix multiply (using BLAS)
Definition: MatrixOps.cpp:170
void normalizeColumns(DoubleMatrix &A)
Normalizes each column as a column-vector.
Definition: MatrixOps.cpp:418