Muster
 All Classes Namespaces Files Functions Variables Typedefs Macros
matrix_utils.h
Go to the documentation of this file.
1 #ifndef MATRIX_UTILS_H
2 #define MATRIX_UTILS_H
3 
4 #include <iostream>
5 #include <iomanip>
6 #include <sstream>
7 #include <string>
8 #include <vector>
9 #include <cfloat>
10 #include <algorithm>
11 #include <numeric>
12 #include <stdint.h>
13 #include <boost/numeric/ublas/matrix.hpp>
14 
15 
16 /// True if and only if n is divisible by 2 <level> times.
17 bool isDivisibleBy2(size_t n, int level);
18 
19 
20 template <class Matrix>
21 bool in_bounds(const Matrix& mat, size_t row, size_t col) {
22  return (row < mat.size1()) && (col < mat.size2());
23 }
24 
25 
26 bool read_matrix(const char *filename, boost::numeric::ublas::matrix<double>& mat);
27 
28 
29 template <class Matrix>
30 void output(const Matrix& mat, std::ostream& out = std::cout) {
31  const size_t width = 12;
32 
33  for (size_t i=0; i < mat.size1(); i++) {
34  if (mat.size2()) {
35  out << std::setw(width) << mat(i,0);
36  }
37 
38  for (size_t j=1; j < mat.size2(); j++) {
39  out << " " << std::setw(width) << mat(i,j);
40  }
41  out << std::endl;
42  }
43 }
44 
45 
46 struct ms_summary {
47  double sum_squares;
48  double orig_max, orig_min;
50  double both_max, both_min;
51 
52  ms_summary(double ss, double ma, double mi)
53  : sum_squares(ss),
54  orig_max(ma), orig_min(mi),
55  repro_max(ma), repro_min(mi),
56  both_max(ma), both_min(mi)
57  { }
58 };
59 
60 
61 template <class Matrix>
62 ms_summary get_summary(const Matrix& orig, const Matrix& repro,
63  size_t row_start = 0, size_t row_end = std::numeric_limits<size_t>::max(),
64  size_t col_start = 0, size_t col_end = std::numeric_limits<size_t>::max()) {
65  assert(orig.size1() == repro.size1());
66  assert(orig.size2() == repro.size2());
67 
68  if (row_end > orig.size1()) row_end = orig.size1();
69  if (col_end > orig.size2()) col_end = orig.size2();
70 
71  ms_summary summary(0, -DBL_MAX, DBL_MAX);
72 
73  for (size_t i=row_start; i < row_end; i++) {
74  for (size_t j=col_start; j < col_end; j++) {
75  double diff = repro(i,j) - orig(i,j);
76  summary.sum_squares += diff*diff;
77 
78  summary.orig_max = std::max(summary.orig_max, orig(i,j));
79  summary.orig_min = std::min(summary.orig_min, orig(i,j));
80  summary.repro_max = std::max(summary.repro_max, repro(i,j));
81  summary.repro_min = std::min(summary.repro_min, repro(i,j));
82  summary.both_max = std::max(summary.orig_max, summary.repro_max);
83  summary.both_min = std::min(summary.orig_min, summary.repro_min);
84  }
85  }
86 
87  return summary;
88 }
89 
90 
91 template <class Matrix>
92 double rmse(const Matrix& orig, const Matrix& repro,
93  size_t row_start = 0, size_t row_end = std::numeric_limits<size_t>::max(),
94  size_t col_start = 0, size_t col_end = std::numeric_limits<size_t>::max()) {
95  assert(orig.size1() == repro.size1());
96  assert(orig.size2() == repro.size2());
97 
98  if (row_end > orig.size1()) row_end = orig.size1();
99  if (col_end > orig.size2()) col_end = orig.size2();
100 
101  ms_summary summary = get_summary(orig, repro, row_start, row_end, col_start, col_end);
102  size_t rows = (row_end - row_start);
103  size_t cols = (col_end - col_start);
104  return sqrt(summary.sum_squares / (rows * cols));
105 }
106 
107 
108 /// Normalized rms error
109 template <class Matrix>
110 double nrmse(const Matrix& orig, const Matrix& repro,
111  size_t row_start = 0, size_t row_end = std::numeric_limits<size_t>::max(),
112  size_t col_start = 0, size_t col_end = std::numeric_limits<size_t>::max()) {
113  assert(orig.size1() == repro.size1());
114  assert(orig.size2() == repro.size2());
115 
116  if (row_end > orig.size1()) row_end = orig.size1();
117  if (col_end > orig.size2()) col_end = orig.size2();
118 
119  ms_summary summary = get_summary(orig, repro, row_start, row_end, col_start, col_end);
120  double range = summary.orig_max - summary.orig_min;
121  size_t rows = (row_end - row_start);
122  size_t cols = (col_end - col_start);
123  return sqrt(summary.sum_squares / (rows * cols)) / range;
124 }
125 
126 
127 /// Peak Signal to Noise Ratio
128 template <class Matrix>
129 double psnr(const Matrix& orig, const Matrix& repro) {
130  assert(orig.size1() == repro.size1());
131  assert(orig.size2() == repro.size2());
132 
133  ms_summary summary = get_summary(orig, repro);
134  double range = summary.orig_max - summary.orig_min;
135  return 20 * log10(range / sqrt(summary.sum_squares / (orig.size1() * orig.size2())));
136 }
137 
138 
139 /// Similarity (NRMSE defined to be symmetric, with max and min taken from both matrices)
140 template <class Matrix>
141 double similarity(const Matrix& orig, const Matrix& repro) {
142  assert(orig.size1() == repro.size1());
143  assert(orig.size2() == repro.size2());
144 
145  ms_summary summary = get_summary(orig, repro);
146  double range = summary.both_max - summary.both_min;
147  return sqrt(summary.sum_squares / (orig.size1() * orig.size2())) / range;
148 }
149 
150 
151 template <class Matrix>
152 void standardize(Matrix& mat) {
153  size_t n = mat.size1() * mat.size2();
154  double sum = 0;
155  double sum2 = 0;
156  for (size_t i=0; i < mat.size1(); i++) {
157  for (size_t j=0; j < mat.size2(); j++) {
158  sum += mat(i,j);
159  sum2 += mat(i,j) * mat(i,j);
160  }
161  }
162  double mean = sum / n;
163  double stdDevInv = 1/sqrt(sum2/n - mean*mean);
164 
165  for (size_t i=0; i < mat.size1(); i++) {
166  for (size_t j=0; j < mat.size2(); j++) {
167  mat(i,j) = (mat(i,j) - mean) * stdDevInv;
168  }
169  }
170 }
171 
172 
173 /// Template class for generic, type-inferred absolute value
174 template <typename T>
175 struct Abs {
176  T value;
177  Abs(double num) : value(::fabs(num)) { }
178  Abs(long num) : value(::labs(num)) { }
179  Abs(unsigned long num) : value(::labs(num)) { }
180  Abs(unsigned int num) : value(::labs(num)) { }
181  Abs(long long num) : value(::llabs(num)) { }
182  Abs(int num) : value(::abs(num)) { }
183 };
184 
185 
186 /// Generic, type-inferred absolute value function.
187 /// Applies fabs to doubles, labs to longs, etc.
188 template<typename T> T abs_val(T num) {
189  return Abs<T>(num).value;
190 }
191 
192 
193 
194 template <class Matrix>
195 typename Matrix::value_type sum(const Matrix& mat,
196  size_t row_start = 0, size_t row_end = std::numeric_limits<size_t>::max(),
197  size_t col_start = 0, size_t col_end = std::numeric_limits<size_t>::max())
198 {
199  if (row_end > mat.size1()) row_end = mat.size1();
200  if (col_end > mat.size2()) col_end = mat.size2();
201 
202  typename Matrix::value_type total = 0;
203  for (size_t i = row_start; i < row_end; i++) {
204  for (size_t j = col_start; j < col_end; j++) {
205  total += mat(i,j);
206  }
207  }
208  return total;
209 }
210 
211 
212 template <class Matrix>
213 double mean_val(const Matrix& mat,
214  size_t row_start = 0, size_t row_end = std::numeric_limits<size_t>::max(),
215  size_t col_start = 0, size_t col_end = std::numeric_limits<size_t>::max())
216 {
217  if (row_end > mat.size1()) row_end = mat.size1();
218  if (col_end > mat.size2()) col_end = mat.size2();
219 
220  size_t row_len = (row_end - row_start);
221  size_t col_len = (col_end - col_start);
222  return sum(mat, row_start, row_end, col_start, col_end)
223  / ((double)row_len * col_len);
224 }
225 
226 
227 template <class Matrix>
228 typename Matrix::value_type max_val(const Matrix& mat,
229  size_t row_start = 0, size_t row_end = std::numeric_limits<size_t>::max(),
230  size_t col_start = 0, size_t col_end = std::numeric_limits<size_t>::max())
231 {
232  if (row_end > mat.size1()) row_end = mat.size1();
233  if (col_end > mat.size2()) col_end = mat.size2();
234  if (row_end <= row_start || col_end <= col_start) return 0;
235 
236  typename Matrix::value_type max = mat(row_start, col_start);
237  for (size_t i = row_start; i < row_end; i++) {
238  for (size_t j = col_start; j < col_end; j++) {
239  max = mat(i,j) > max ? mat(i,j) : max;
240  }
241  }
242  return max;
243 }
244 
245 
246 template <class Matrix>
247 typename Matrix::value_type min_val(const Matrix& mat,
248  size_t row_start = 0, size_t row_end = std::numeric_limits<size_t>::max(),
249  size_t col_start = 0, size_t col_end = std::numeric_limits<size_t>::max())
250 {
251  if (row_end > mat.size1()) row_end = mat.size1();
252  if (col_end > mat.size2()) col_end = mat.size2();
253  if (row_end <= row_start || col_end <= col_start) return 0;
254 
255  typename Matrix::value_type min = mat(row_start, col_start);
256  for (size_t i = row_start; i < row_end; i++) {
257  for (size_t j = col_start; j < col_end; j++) {
258  min = mat(i,j) < min ? mat(i,j) : min;
259  }
260  }
261  return min;
262 }
263 
264 template <class Matrix>
265 typename Matrix::value_type abs_max_val(const Matrix& mat,
266  size_t row_start = 0, size_t row_end = std::numeric_limits<size_t>::max(),
267  size_t col_start = 0, size_t col_end = std::numeric_limits<size_t>::max())
268 {
269  if (row_end > mat.size1()) row_end = mat.size1();
270  if (col_end > mat.size2()) col_end = mat.size2();
271 
272  typename Matrix::value_type amax = 0;
273  for (size_t i = row_start; i < row_end; i++) {
274  for (size_t j = col_start; j < col_end; j++) {
275  typename Matrix::value_type tmp = abs_val(mat(i,j));
276  amax = (tmp > amax) ? tmp : amax;
277  }
278  }
279  return amax;
280 }
281 
282 template <class Matrix, typename V>
283 void set_all(const Matrix& mat, V value,
284  size_t row_start = 0, size_t row_end = std::numeric_limits<size_t>::max(),
285  size_t col_start = 0, size_t col_end = std::numeric_limits<size_t>::max())
286 {
287  if (row_end > mat.size1()) row_end = mat.size1();
288  if (col_end > mat.size2()) col_end = mat.size2();
289 
290  for (size_t i = row_start; i < row_end; i++) {
291  for (size_t j = col_start; j < col_end; j++) {
292  mat(i,j) = value;
293  }
294  }
295 }
296 
297 
298 template <typename Iterator1, typename Iterator2>
299 double manhattan_distance(Iterator1 first1, Iterator1 last1, Iterator2 first2) {
300  double sum = 0.0;
301  Iterator1 i = first1;
302  Iterator2 j = first2;
303  while (i != last1) {
304  sum += abs_val(*i - *j);
305  i++;
306  j++;
307  }
308  return sum;
309 }
310 
311 
312 template <typename Iterator1, typename Iterator2>
313 double euclidean_distance(Iterator1 first1, Iterator1 last1, Iterator2 first2) {
314  double sum = 0.0;
315  Iterator1 i = first1;
316  Iterator2 j = first2;
317  while (i != last1) {
318  double diff = (*i - *j);
319  sum += diff * diff;
320  i++;
321  j++;
322  }
323  return sqrt(sum);
324 }
325 
326 
327 /// Interpolates a value for point (x,y) based on values in the Matrix.
328 /// Finds the nearest values by taking floor and ceil of x and y, then
329 /// uses bilinear interpolation to estimate the value at (x,y).
330 template <class Matrix>
331 double interp_bilinear(const Matrix& mat, double x, double y) {
332  if (x < 0) {
333  x = 0;
334  } else if (x > (mat.size1() - 1)) {
335  x = mat.size1() - 1;
336  }
337 
338  if (y < 0) {
339  y = 0;
340  } else if (y > (mat.size2() - 1)) {
341  y = mat.size2() - 1;
342  }
343 
344  // get corners of interpolation
345  size_t x1 = (size_t)floor(x);
346  size_t x2 = (size_t)ceil(x);
347 
348  size_t y1 = (size_t)floor(y);
349  size_t y2 = (size_t)ceil(y);
350 
351  size_t xd = (x2-x1);
352  size_t yd = (y2-y1);
353 
354  double value;
355  if (xd == 0 && yd == 0) {
356  value = mat(x1,y1);
357 
358  } else if (yd == 0) {
359  value = (x2 - x) * mat(x1,y1) + (x - x1) * mat(x2, y2);
360 
361  } else if (xd == 0) {
362  value = (y2 - y) * mat(x1,y1) + (y - y1) * mat(x2, y2);
363 
364  } else {
365  double d = xd * yd;
366  value =
367  mat(x1,y1)/d * ((x2-x) * (y2-y)) +
368  mat(x2,y1)/d * ((x-x1) * (y2-y)) +
369  mat(x1,y2)/d * ((x2-x) * (y-y1)) +
370  mat(x2,y2)/d * ((x-x1) * (y-y1));
371  }
372 
373  return value;
374 }
375 
376 #endif // MATRIX_UTILS_H
double repro_min
Definition: matrix_utils.h:49
Abs(unsigned long num)
Definition: matrix_utils.h:179
void output(const Matrix &mat, std::ostream &out=std::cout)
Definition: matrix_utils.h:30
Matrix::value_type max_val(const Matrix &mat, size_t row_start=0, size_t row_end=std::numeric_limits< size_t >::max(), size_t col_start=0, size_t col_end=std::numeric_limits< size_t >::max())
Definition: matrix_utils.h:228
T value
Definition: matrix_utils.h:176
double sum_squares
Definition: matrix_utils.h:47
double orig_max
Definition: matrix_utils.h:48
double orig_min
Definition: matrix_utils.h:48
double euclidean_distance(Iterator1 first1, Iterator1 last1, Iterator2 first2)
Definition: matrix_utils.h:313
Abs(unsigned int num)
Definition: matrix_utils.h:180
double psnr(const Matrix &orig, const Matrix &repro)
Peak Signal to Noise Ratio.
Definition: matrix_utils.h:129
double interp_bilinear(const Matrix &mat, double x, double y)
Interpolates a value for point (x,y) based on values in the Matrix.
Definition: matrix_utils.h:331
Matrix::value_type sum(const Matrix &mat, size_t row_start=0, size_t row_end=std::numeric_limits< size_t >::max(), size_t col_start=0, size_t col_end=std::numeric_limits< size_t >::max())
Definition: matrix_utils.h:195
Matrix::value_type min_val(const Matrix &mat, size_t row_start=0, size_t row_end=std::numeric_limits< size_t >::max(), size_t col_start=0, size_t col_end=std::numeric_limits< size_t >::max())
Definition: matrix_utils.h:247
double rmse(const Matrix &orig, const Matrix &repro, size_t row_start=0, size_t row_end=std::numeric_limits< size_t >::max(), size_t col_start=0, size_t col_end=std::numeric_limits< size_t >::max())
Definition: matrix_utils.h:92
ms_summary get_summary(const Matrix &orig, const Matrix &repro, size_t row_start=0, size_t row_end=std::numeric_limits< size_t >::max(), size_t col_start=0, size_t col_end=std::numeric_limits< size_t >::max())
Definition: matrix_utils.h:62
bool read_matrix(const char *filename, boost::numeric::ublas::matrix< double > &mat)
bool in_bounds(const Matrix &mat, size_t row, size_t col)
Definition: matrix_utils.h:21
Abs(long num)
Definition: matrix_utils.h:178
double both_min
Definition: matrix_utils.h:50
double repro_max
Definition: matrix_utils.h:49
void set_all(const Matrix &mat, V value, size_t row_start=0, size_t row_end=std::numeric_limits< size_t >::max(), size_t col_start=0, size_t col_end=std::numeric_limits< size_t >::max())
Definition: matrix_utils.h:283
void standardize(Matrix &mat)
Definition: matrix_utils.h:152
double mean_val(const Matrix &mat, size_t row_start=0, size_t row_end=std::numeric_limits< size_t >::max(), size_t col_start=0, size_t col_end=std::numeric_limits< size_t >::max())
Definition: matrix_utils.h:213
Template class for generic, type-inferred absolute value.
Definition: matrix_utils.h:175
Abs(int num)
Definition: matrix_utils.h:182
double both_max
Definition: matrix_utils.h:50
double similarity(const Matrix &orig, const Matrix &repro)
Similarity (NRMSE defined to be symmetric, with max and min taken from both matrices) ...
Definition: matrix_utils.h:141
Abs(long long num)
Definition: matrix_utils.h:181
T abs_val(T num)
Generic, type-inferred absolute value function.
Definition: matrix_utils.h:188
double manhattan_distance(Iterator1 first1, Iterator1 last1, Iterator2 first2)
Definition: matrix_utils.h:299
Abs(double num)
Definition: matrix_utils.h:177
ms_summary(double ss, double ma, double mi)
Definition: matrix_utils.h:52
Matrix::value_type abs_max_val(const Matrix &mat, size_t row_start=0, size_t row_end=std::numeric_limits< size_t >::max(), size_t col_start=0, size_t col_end=std::numeric_limits< size_t >::max())
Definition: matrix_utils.h:265
double nrmse(const Matrix &orig, const Matrix &repro, size_t row_start=0, size_t row_end=std::numeric_limits< size_t >::max(), size_t col_start=0, size_t col_end=std::numeric_limits< size_t >::max())
Normalized rms error.
Definition: matrix_utils.h:110
bool isDivisibleBy2(size_t n, int level)
True if and only if n is divisible by 2 &lt;level&gt; times.
Muster. Copyright © 2010, Lawrence Livermore National Laboratory, LLNL-CODE-433662.
Distribution of Muster and its documentation is subject to terms of the Muster LICENSE.
Generated on Thu Sep 1 2016 using Doxygen 1.8.5