Pteros  2.0
Molecular modeling library for human beings!
matrix.h
1 // Copyright Maarten L. Hekkelman, Radboud University 2008-2011.
2 // Distributed under the Boost Software License, Version 1.0.
3 // (See accompanying file LICENSE_1_0.txt or copy at
4 // http://www.boost.org/LICENSE_1_0.txt)
5 //
6 // substitution matrix for multiple sequence alignments
7 
8 #pragma once
9 
10 #include "mas.h"
11 
12 #include <string>
13 #include <istream>
14 #include <cassert>
15 #include <stdexcept>
16 #include <algorithm>
17 
18 // Some predefined matrices
19 
20 // PAM250 is used by hssp-nt in aligning the sequences
21 
22 extern const int8 kMBlosum45[], kMBlosum50[], kMBlosum62[], kMBlosum80[], kMBlosum90[],
23  kMPam250[], kMPam30[], kMPam70[];
24 extern const float kMPam250ScalingFactor, kMPam250MisMatchAverage;
25 
26 struct MMtrxStats
27 {
28  double lambda, kappa, entropy, alpha, beta;
29 };
30 
31 struct MMatrixData
32 {
33  const char* mName;
34  int8 mGapOpen, mGapExtend;
35  const int8* mMatrix;
36  MMtrxStats mGappedStats, mUngappedStats;
37 };
38 
39 extern const MMatrixData kMMatrixData[];
40 
41 // Dayhoff matrix is used for calculating similarity in HSSP
42 
43 extern const float kDayhoffData[];
44 
45 // Simple scoring function using the predefined matrices
46 template<typename T>
47 inline T score(const T inMatrix[], uint8 inAA1, uint8 inAA2)
48 {
49  T result;
50 
51  if (inAA1 >= inAA2)
52  result = inMatrix[(inAA1 * (inAA1 + 1)) / 2 + inAA2];
53  else
54  result = inMatrix[(inAA2 * (inAA2 + 1)) / 2 + inAA1];
55 
56  return result;
57 }
58 
59 // --------------------------------------------------------------------
60 // uBlas compatible matrix types
61 // matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
62 // element m i,j is mapped to [i * n + j] and thus storage is row major
63 
64 template<typename T>
65 class matrix_base
66 {
67  public:
68 
69  typedef T value_type;
70 
71  virtual ~matrix_base() {}
72 
73  virtual uint32 dim_m() const = 0;
74  virtual uint32 dim_n() const = 0;
75 
76  virtual value_type& operator()(uint32 i, uint32 j) { throw std::runtime_error("unimplemented method"); }
77  virtual value_type operator()(uint32 i, uint32 j) const = 0;
78 
79  matrix_base& operator*=(const value_type& rhs);
80 
81  matrix_base& operator-=(const value_type& rhs);
82 };
83 
84 template<typename T>
85 matrix_base<T>& matrix_base<T>::operator*=(const T& rhs)
86 {
87  for (uint32 i = 0; i < dim_m(); ++i)
88  {
89  for (uint32 j = 0; j < dim_n(); ++j)
90  {
91  operator()(i, j) *= rhs;
92  }
93  }
94 
95  return *this;
96 }
97 
98 template<typename T>
99 matrix_base<T>& matrix_base<T>::operator-=(const T& rhs)
100 {
101  for (uint32 i = 0; i < dim_m(); ++i)
102  {
103  for (uint32 j = 0; j < dim_n(); ++j)
104  {
105  operator()(i, j) -= rhs;
106  }
107  }
108 
109  return *this;
110 }
111 
112 template<typename T>
113 std::ostream& operator<<(std::ostream& lhs, const matrix_base<T>& rhs)
114 {
115  lhs << '[' << rhs.dim_m() << ',' << rhs.dim_n() << ']' << '(';
116  for (uint32 i = 0; i < rhs.dim_m(); ++i)
117  {
118  lhs << '(';
119  for (uint32 j = 0; j < rhs.dim_n(); ++j)
120  {
121  if (j > 0)
122  lhs << ',';
123  lhs << rhs(i,j);
124  }
125  lhs << ')';
126  }
127  lhs << ')';
128 
129  return lhs;
130 }
131 
132 template<typename T>
133 class matrix : public matrix_base<T>
134 {
135  public:
136  typedef T value_type;
137 
138  template<typename T2>
139  matrix(const matrix_base<T2>& m)
140  : m_m(m.dim_m())
141  , m_n(m.dim_n())
142  {
143  m_data = new value_type[m_m * m_n];
144  for (uint32 i = 0; i < m_m; ++i)
145  {
146  for (uint32 j = 0; j < m_n; ++j)
147  operator()(i, j) = m(i, j);
148  }
149  }
150 
151  matrix(const matrix& m)
152  : m_m(m.m_m)
153  , m_n(m.m_n)
154  {
155  m_data = new value_type[m_m * m_n];
156  std::copy(m.m_data, m.m_data + (m_m * m_n), m_data);
157  }
158 
159  matrix& operator=(const matrix& m)
160  {
161  value_type t = new value_type[m.m_m * m.m_n];
162  std::copy(m.m_data, m.m_data + (m_m * m_n), t);
163 
164  delete[] m_data;
165  m_data = t;
166  m_m = m.m_m;
167  m_n = m.m_n;
168 
169  return *this;
170  }
171 
172  matrix(uint32 m, uint32 n, T v = T())
173  : m_m(m)
174  , m_n(n)
175  {
176  m_data = new value_type[m_m * m_n];
177  std::fill(m_data, m_data + (m_m * m_n), v);
178  }
179 
180  virtual ~matrix()
181  {
182  delete [] m_data;
183  }
184 
185  virtual uint32 dim_m() const { return m_m; }
186  virtual uint32 dim_n() const { return m_n; }
187 
188  virtual value_type operator()(uint32 i, uint32 j) const
189  {
190  assert(i < m_m); assert(j < m_n);
191  return m_data[i * m_n + j];
192  }
193 
194  virtual value_type& operator()(uint32 i, uint32 j)
195  {
196  assert(i < m_m); assert(j < m_n);
197  return m_data[i * m_n + j];
198  }
199 
200  private:
201  value_type* m_data;
202  uint32 m_m, m_n;
203 };
204 
205 // --------------------------------------------------------------------
206 
207 template<typename T>
208 class symmetric_matrix : public matrix_base<T>
209 {
210  public:
211  typedef typename matrix_base<T>::value_type value_type;
212 
213  symmetric_matrix(uint32 n)
214  : m_owner(true)
215  , m_n(n)
216  {
217  uint32 N = (m_n * (m_n + 1)) / 2;
218  m_data = new value_type[N];
219  std::fill(m_data, m_data + N, T(0));
220  }
221 
222  symmetric_matrix(const T* data, uint32 n)
223  : m_owner(false)
224  , m_data(const_cast<T*>(data))
225  , m_n(n)
226  {
227  }
228 
229 
230  virtual ~symmetric_matrix()
231  {
232  if (m_owner)
233  delete[] m_data;
234  }
235 
236  virtual uint32 dim_m() const { return m_n; }
237  virtual uint32 dim_n() const { return m_n; }
238 
239  T operator()(uint32 i, uint32 j) const;
240  virtual T& operator()(uint32 i, uint32 j);
241 
242  // erase two rows, add one at the end (for neighbour joining)
243  void erase_2(uint32 i, uint32 j);
244 
245  private:
246  bool m_owner;
247  value_type* m_data;
248  uint32 m_n;
249 };
250 
251 template<typename T>
252 inline
253 T symmetric_matrix<T>::operator()(uint32 i, uint32 j) const
254 {
255  return i < j
256  ? m_data[(j * (j + 1)) / 2 + i]
257  : m_data[(i * (i + 1)) / 2 + j];
258 // if (i > j)
259 // std::swap(i, j);
260 // assert(j < m_n);
261 // return m_data[(j * (j + 1)) / 2 + i];
262 }
263 
264 template<typename T>
265 inline
266 T& symmetric_matrix<T>::operator()(uint32 i, uint32 j)
267 {
268  if (i > j)
269  std::swap(i, j);
270  assert(j < m_n);
271  return m_data[(j * (j + 1)) / 2 + i];
272 }
273 
274 template<typename T>
275 void symmetric_matrix<T>::erase_2(uint32 di, uint32 dj)
276 {
277  uint32 s = 0, d = 0;
278  for (uint32 i = 0; i < m_n; ++i)
279  {
280  for (uint32 j = 0; j < i; ++j)
281  {
282  if (i != di and j != dj and i != dj and j != di)
283  {
284  if (s != d)
285  m_data[d] = m_data[s];
286  ++d;
287  }
288 
289  ++s;
290  }
291  }
292 
293  --m_n;
294 }
295 
296 template<typename T>
297 class identity_matrix : public matrix_base<T>
298 {
299  public:
300  typedef typename matrix_base<T>::value_type value_type;
301 
302  identity_matrix(uint32 n)
303  : m_n(n)
304  {
305  }
306 
307  virtual uint32 dim_m() const { return m_n; }
308  virtual uint32 dim_n() const { return m_n; }
309 
310  virtual value_type operator()(uint32 i, uint32 j) const
311  {
312  value_type result = 0;
313  if (i == j)
314  result = 1;
315  return result;
316  }
317 
318  private:
319  uint32 m_n;
320 };
321 
322 // --------------------------------------------------------------------
323 // matrix functions
324 
325 template<typename T>
326 matrix<T> operator*(const matrix_base<T>& lhs, const matrix_base<T>& rhs)
327 {
328  matrix<T> result(std::min(lhs.dim_m(), rhs.dim_m()), std::min(lhs.dim_n(), rhs.dim_n()));
329 
330  for (uint32 i = 0; i < result.dim_m(); ++i)
331  {
332  for (uint32 j = 0; j < result.dim_n(); ++j)
333  {
334  for (uint32 li = 0, rj = 0; li < lhs.dim_m() and rj < rhs.dim_n(); ++li, ++rj)
335  result(i, j) += lhs(li, j) * rhs(i, rj);
336  }
337  }
338 
339  return result;
340 }
341 
342 template<typename T>
343 matrix<T> operator*(const matrix_base<T>& lhs, T rhs)
344 {
345  matrix<T> result(lhs);
346  result *= rhs;
347 
348  return result;
349 }
350 
351 template<typename T>
352 matrix<T> operator-(const matrix_base<T>& lhs, const matrix_base<T>& rhs)
353 {
354  matrix<T> result(std::min(lhs.dim_m(), rhs.dim_m()), std::min(lhs.dim_n(), rhs.dim_n()));
355 
356  for (uint32 i = 0; i < result.dim_m(); ++i)
357  {
358  for (uint32 j = 0; j < result.dim_n(); ++j)
359  {
360  result(i, j) = lhs(i, j) - rhs(i, j);
361  }
362  }
363 
364  return result;
365 }
366 
367 template<typename T>
368 matrix<T> operator-(const matrix_base<T>& lhs, T rhs)
369 {
370  matrix<T> result(lhs.dim_m(), lhs.dim_n());
371  result -= rhs;
372  return result;
373 }
374 
376 //
377 //class substitution_matrix
378 //{
379 // public:
380 // substitution_matrix(const std::string& name);
381 // substitution_matrix(
382 // const substitution_matrix& m, bool positive);
383 //
384 // virtual ~substitution_matrix() {}
385 //
386 // int8 operator()(aa a, aa b) const
387 // {
388 // return m_matrix(a, b);
389 // }
390 //
391 // int8 operator()(char a, char b) const
392 // {
393 // return m_matrix(encode(a), encode(b));
394 // }
395 //
396 // float mismatch_average() const { return m_mismatch_average; }
397 // float scale_factor() const { return m_scale_factor; }
398 //
399 // private:
400 // substitution_matrix(
401 // const substitution_matrix&);
402 // substitution_matrix&
403 // operator=(const substitution_matrix&);
404 //
405 // void read(std::istream& is);
406 //
407 // matrix<int8> m_matrix;
408 // float m_mismatch_average;
409 // float m_scale_factor;
410 //};
411 //
412 //class substitution_matrix_family
413 //{
414 // public:
415 // substitution_matrix_family(
416 // const std::string& name);
417 //
418 // ~substitution_matrix_family();
419 //
420 // const substitution_matrix&
421 // operator()(float distance, bool positive) const
422 // {
423 // const substitution_matrix* result;
424 //
425 // uint32 ix = 0;
426 // while (distance < m_cutoff[ix] and ix < 3)
427 // ++ix;
428 //
429 // if (positive)
430 // result = m_pos_smat[ix];
431 // else
432 // result = m_smat[ix];
433 //
434 // return *result;
435 // }
436 //
437 // private:
438 // substitution_matrix_family(
439 // const substitution_matrix_family&);
440 // substitution_matrix_family&
441 // operator=(const substitution_matrix_family&);
442 //
443 // float m_cutoff[4];
444 // substitution_matrix*
445 // m_smat[4];
446 // substitution_matrix*
447 // m_pos_smat[4];
448 //};
449 //