← Back to branch summary

~chris-rogers/maus/emr_mc_digitization

« back to all changes in this revision

Viewing changes to src/common_cpp/Maths/Matrix.cc

  • Committer: Chris Rogers
  • Date: 2014-04-16 11:48:45 UTC
  • mfrom: (707 merge)
  • mto: This revision was merged to the branch mainline in revision 711.
  • Revision ID: chris.rogers@stfc.ac.uk-20140416114845-h3u3q7pdcxkxvovs
Update to trunk

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/common_cpp/Maths/Matrix.cc
38
38
#include "gsl/gsl_eigen.h"
39
39
#include "json/json.h"
40
40
 
41
 
#include "Interface/Squeal.hh"
 
41
#include "Utils/Exception.hh"
42
42
#include "Maths/Vector.hh"
43
43
#include "Utils/JsonWrapper.hh"
44
44
 
167
167
double& MatrixBase<double, gsl_matrix>::operator()(
168
168
  const size_t row, const size_t column) {
169
169
  if (matrix_ == NULL) {
170
 
    throw(Squeal(Squeal::recoverable,
 
170
    throw(Exception(Exception::recoverable,
171
171
                 "Attempting to index an empty matrix.",
172
172
                 "MatrixBase<double, gsl_matrix>::operator()"));
173
173
  } else if ((row < 1) || (row > matrix_->size1)) {
174
 
    throw(Squeal(Squeal::recoverable,
 
174
    throw(Exception(Exception::recoverable,
175
175
                 "Row index out of bounds.",
176
176
                 "MatrixBase<double, gsl_matrix>::operator()"));
177
177
  } else if ((column < 1) || (column > matrix_->size2)) {
178
 
    throw(Squeal(Squeal::recoverable,
 
178
    throw(Exception(Exception::recoverable,
179
179
                 "Column index out of bounds.",
180
180
                 "MatrixBase<double, gsl_matrix>::operator()"));
181
181
  }
186
186
complex& MatrixBase<complex, gsl_matrix_complex>::operator()(
187
187
    const size_t row, const size_t column) {
188
188
  if (matrix_ == NULL) {
189
 
    throw(Squeal(Squeal::recoverable,
 
189
    throw(Exception(Exception::recoverable,
190
190
                 "Attempting to index an empty matrix.",
191
191
                 "MatrixBase<complex, gsl_matrix_complex>::operator()"));
192
192
  } else if ((row < 1) || (row > matrix_->size1)) {
193
 
    throw(Squeal(Squeal::recoverable,
 
193
    throw(Exception(Exception::recoverable,
194
194
                 "Row index out of bounds.",
195
195
                 "MatrixBase<complex, gsl_matrix_complex>::operator()"));
196
196
  } else if ((column < 1) || (column > matrix_->size2)) {
197
 
    throw(Squeal(Squeal::recoverable,
 
197
    throw(Exception(Exception::recoverable,
198
198
                 "Column index out of bounds.",
199
199
                 "MatrixBase<complex, gsl_matrix_complex>::operator()"));
200
200
  }
205
205
const double& MatrixBase<double, gsl_matrix>::operator()(
206
206
  const size_t row, const size_t column) const {
207
207
  if (matrix_ == NULL) {
208
 
    throw(Squeal(Squeal::recoverable,
 
208
    throw(Exception(Exception::recoverable,
209
209
                 "Attempting to index an empty matrix.",
210
210
                 "MatrixBase<double, gsl_matrix>::operator()"));
211
211
  } else if ((row < 1) || (row > matrix_->size1)) {
212
 
    throw(Squeal(Squeal::recoverable,
 
212
    throw(Exception(Exception::recoverable,
213
213
                 "Row index out of bounds.",
214
214
                 "MatrixBase<double, gsl_matrix>::operator()"));
215
215
  } else if ((column < 1) || (column > matrix_->size2)) {
216
 
    throw(Squeal(Squeal::recoverable,
 
216
    throw(Exception(Exception::recoverable,
217
217
                 "Column index out of bounds.",
218
218
                 "MatrixBase<double, gsl_matrix>::operator()"));
219
219
  }
224
224
const complex& MatrixBase<complex, gsl_matrix_complex>::operator()(
225
225
  const size_t row, const size_t column) const {
226
226
  if (matrix_ == NULL) {
227
 
    throw(Squeal(Squeal::recoverable,
 
227
    throw(Exception(Exception::recoverable,
228
228
                 "Attempting to index an empty matrix.",
229
229
                 "MatrixBase<complex, gsl_matrix_complex>::operator()"));
230
230
  } else if ((row < 1) || (row > matrix_->size1)) {
231
 
    throw(Squeal(Squeal::recoverable,
 
231
    throw(Exception(Exception::recoverable,
232
232
                 "Row index out of bounds.",
233
233
                 "MatrixBase<complex, gsl_matrix_complex>::operator()"));
234
234
  } else if ((column < 1) || (column > matrix_->size2)) {
235
 
    throw(Squeal(Squeal::recoverable,
 
235
    throw(Exception(Exception::recoverable,
236
236
                 "Column index out of bounds.",
237
237
                 "MatrixBase<complex, gsl_matrix_complex>::operator()"));
238
238
  }
242
242
template <typename StdType, typename GslType> Vector<StdType>
243
243
MAUS::MatrixBase<StdType, GslType>::row(const size_t row) const {
244
244
  if (matrix_ == NULL) {
245
 
    throw(Squeal(Squeal::recoverable,
 
245
    throw(Exception(Exception::recoverable,
246
246
                 "Attempting to index an empty matrix.",
247
247
                 "MatrixBase<double, gsl_matrix>::row()"));
248
248
  } else if ((row < 1) || (row > matrix_->size1)) {
249
 
    throw(Squeal(Squeal::recoverable,
 
249
    throw(Exception(Exception::recoverable,
250
250
                 "Row index out of bounds.",
251
251
                 "MatrixBase<double, gsl_matrix>::operator()"));
252
252
  }
266
266
template <typename StdType, typename GslType> Vector<StdType>
267
267
MatrixBase<StdType, GslType>::column(const size_t column) const {
268
268
  if (matrix_ == NULL) {
269
 
    throw(Squeal(Squeal::recoverable,
 
269
    throw(Exception(Exception::recoverable,
270
270
                 "Attempting to index and empty matrix.",
271
271
                 "MatrixBase<double, gsl_matrix>::operator()"));
272
272
  } else if ((column < 1) || (column > matrix_->size2)) {
273
 
    throw(Squeal(Squeal::recoverable,
 
273
    throw(Exception(Exception::recoverable,
274
274
                 "Column index out of bounds.",
275
275
                 "MatrixBase<double, gsl_matrix>::operator()"));
276
276
  }
355
355
    } else if (   (rhs.matrix_ == NULL)
356
356
               || (matrix_->size1 != rhs.matrix_->size1)
357
357
               || (matrix_->size2 != rhs.matrix_->size2)) {
358
 
      throw(Squeal(Squeal::recoverable,
 
358
      throw(Exception(Exception::recoverable,
359
359
                   "Attempted to assign a matrix of a different size.",
360
360
                   "MatrixBase<double>::operator=()"));
361
361
    }
377
377
    } else if (   (rhs.matrix_ == NULL)
378
378
               || (matrix_->size1 != rhs.matrix_->size1)
379
379
               || (matrix_->size2 != rhs.matrix_->size2)) {
380
 
      throw(Squeal(Squeal::recoverable,
 
380
      throw(Exception(Exception::recoverable,
381
381
                   "Attempted to assign a matrix of a different size.",
382
382
                   "MatrixBase<complex>::operator=()"));
383
383
    }
392
392
  const MatrixBase<double, gsl_matrix>& rhs) {
393
393
  if (   (number_of_rows() != rhs.number_of_rows())
394
394
      || (number_of_columns() != rhs.number_of_columns())) {
395
 
    throw(Squeal(Squeal::recoverable,
 
395
    throw(Exception(Exception::recoverable,
396
396
                 "Attempted to add two matrices of different sizes",
397
397
                 "MatrixBase<double>::operator+=()"));
398
398
  }
407
407
  const MatrixBase<complex, gsl_matrix_complex>& rhs) {
408
408
  if (   (number_of_rows() != rhs.number_of_rows())
409
409
      || (number_of_columns() != rhs.number_of_columns())) {
410
 
    throw(Squeal(Squeal::recoverable,
 
410
    throw(Exception(Exception::recoverable,
411
411
                 "Attempted to add two matrices of different sizes",
412
412
                 "MatrixBase<complex>::operator+=()"));
413
413
  }
422
422
  const MatrixBase<double, gsl_matrix>& rhs) {
423
423
  if (   (number_of_rows() != rhs.number_of_rows())
424
424
      || (number_of_columns() != rhs.number_of_columns())) {
425
 
    throw(Squeal(Squeal::recoverable,
 
425
    throw(Exception(Exception::recoverable,
426
426
                 "Attempted to subtract two matrices of different sizes",
427
427
                 "MatrixBase<double,gsl_matrix>::operator-=()"));
428
428
  }
437
437
  const MatrixBase<complex, gsl_matrix_complex>& rhs) {
438
438
  if (   (number_of_rows() != rhs.number_of_rows())
439
439
      || (number_of_columns() != rhs.number_of_columns())) {
440
 
    throw(Squeal(Squeal::recoverable,
 
440
    throw(Exception(Exception::recoverable,
441
441
                 "Attempted to subtract two matrices of different sizes",
442
442
                 "MatrixBase<complex,gsl_matrix_complex>::operator-=()"));
443
443
  }
451
451
MatrixBase<double, gsl_matrix>::operator*=(
452
452
  const MatrixBase<double, gsl_matrix>& rhs) {
453
453
  if (number_of_columns() != rhs.number_of_rows()) {
454
 
    throw(Squeal(Squeal::recoverable,
 
454
    throw(Exception(Exception::recoverable,
455
455
                 "Attempted to multiply two matrices of incompatible sizes",
456
456
                 "MatrixBase<complex,gsl_matrix_complex>::operator*=()"));
457
457
  }
465
465
MatrixBase<complex, gsl_matrix_complex>::operator*=(
466
466
  const MatrixBase<complex, gsl_matrix_complex>& rhs) {
467
467
  if (number_of_columns() != rhs.number_of_rows()) {
468
 
    throw(Squeal(Squeal::recoverable,
 
468
    throw(Exception(Exception::recoverable,
469
469
                 "Attempted to multiply two matrices of incompatible sizes",
470
470
                 "MatrixBase<complex,gsl_matrix_complex>::operator*=()"));
471
471
  }
683
683
                                                     const char *   file,
684
684
                                                     int            line,
685
685
                                                     int            gsl_errno) {
686
 
  throw(Squeal(
687
 
      Squeal::recoverable,
 
686
  throw(Exception(
 
687
      Exception::recoverable,
688
688
      reason,
689
689
      "MatrixBase<StdType, GslType>::gsl_error_handler()"));
690
690
}
725
725
  size_t columns = real_matrix.number_of_columns();
726
726
  if (   (rows != imaginary_matrix.number_of_rows())
727
727
      || (columns != imaginary_matrix.number_of_columns())) {
728
 
    throw(Squeal(
729
 
      Squeal::recoverable,
 
728
    throw(Exception(
 
729
      Exception::recoverable,
730
730
      "Attempted to build a complex matrix using "
731
731
      "real and imaginary matrices of different sizes",
732
732
      "Matrix<complex>::Matrix()"));
969
969
  size_t columns = real_matrix.number_of_columns();
970
970
  if (   (rows != imaginary_matrix.number_of_rows())
971
971
      || (columns != imaginary_matrix.number_of_columns())) {
972
 
    throw(Squeal(
973
 
      Squeal::recoverable,
 
972
    throw(Exception(
 
973
      Exception::recoverable,
974
974
      "Attempted to build a complex matrix using "
975
975
      "real and imaginary matrices of different sizes",
976
976
      "MAUS::MAUS::Complex::complex"));
995
995
  size_t rows = matrix.number_of_rows();
996
996
  size_t columns = matrix.number_of_columns();
997
997
  if (rows != columns) {
998
 
    throw(Squeal(Squeal::recoverable,
 
998
    throw(Exception(Exception::recoverable,
999
999
                 "Attempt to get determinant of non-square matrix",
1000
1000
                 "MAUS::determinant()"));
1001
1001
  }
1014
1014
  size_t rows = matrix.number_of_rows();
1015
1015
  size_t columns = matrix.number_of_columns();
1016
1016
  if (rows != columns) {
1017
 
    throw(Squeal(Squeal::recoverable,
 
1017
    throw(Exception(Exception::recoverable,
1018
1018
                 "Attempt to get determinant of non-square matrix",
1019
1019
                 "MAUS::determinant()"));
1020
1020
  }
1031
1031
  size_t rows = matrix.number_of_rows();
1032
1032
  size_t columns = matrix.number_of_columns();
1033
1033
  if (rows != columns) {
1034
 
    throw(Squeal(Squeal::recoverable,
 
1034
    throw(Exception(Exception::recoverable,
1035
1035
                 "Attempted to get inverse of non-square matrix",
1036
1036
                 "MAUS::inverse()"));
1037
1037
  }
1038
1038
  if (determinant(matrix) == 0.) {
1039
 
    throw(Squeal(Squeal::recoverable,
 
1039
    throw(Exception(Exception::recoverable,
1040
1040
                 "Attempted to get inverse of singular matrix",
1041
1041
                 "MAUS::inverse()"));
1042
1042
  }
1053
1053
    for (size_t column = 1; column <= columns; ++column) {
1054
1054
      test_value = matrix(row, column);
1055
1055
      if (test_value != test_value) {
1056
 
        throw(Squeal(Squeal::recoverable,
 
1056
        throw(Exception(Exception::recoverable,
1057
1057
          "Failed to invert matrix: result contained NaN elements - singular?",
1058
1058
          "MAUS::inverse()"));
1059
1059
      } else if (test_value == std::numeric_limits<double>::infinity()) {
1060
 
        throw(Squeal(Squeal::recoverable,
 
1060
        throw(Exception(Exception::recoverable,
1061
1061
          "Failed to invert matrix: "
1062
1062
          "result contained infinite elements - singular?",
1063
1063
          "MAUS::inverse()"));
1073
1073
  size_t rows = matrix.number_of_rows();
1074
1074
  size_t columns = matrix.number_of_columns();
1075
1075
  if (rows != columns) {
1076
 
    throw(Squeal(Squeal::recoverable,
 
1076
    throw(Exception(Exception::recoverable,
1077
1077
                 "Attempt to get inverse of non-square matrix",
1078
1078
                 "MAUS::inverse()"));
1079
1079
  }
1090
1090
    for (size_t column = 1; column <= columns; ++column) {
1091
1091
      test_value = matrix(row, column);
1092
1092
      if (test_value != test_value) {
1093
 
        throw(Squeal(Squeal::recoverable,
 
1093
        throw(Exception(Exception::recoverable,
1094
1094
          "Failed to invert matrix: result contained NaN elements - singular?",
1095
1095
          "MAUS::inverse()"));
1096
1096
      } else if (   real(test_value) == infinity
1097
1097
                 || imag(test_value) == infinity) {
1098
 
        throw(Squeal(Squeal::recoverable,
 
1098
        throw(Exception(Exception::recoverable,
1099
1099
          "Failed to invert matrix: "
1100
1100
          "result contained infinite elements - singular?",
1101
1101
          "MAUS::inverse()"));
1161
1161
  size_t rows = matrix.number_of_rows();
1162
1162
  size_t columns = matrix.number_of_columns();
1163
1163
  if (rows != columns) {
1164
 
    throw(Squeal(Squeal::recoverable,
 
1164
    throw(Exception(Exception::recoverable,
1165
1165
                 "Attempt to get eigenvalues of non-square matrix",
1166
1166
                 "MAUS::eigenvalues"));
1167
1167
  }
1173
1173
  gsl_eigen_nonsymm_free(workspace);
1174
1174
  if (ierr != 0) {
1175
1175
    gsl_vector_complex_free(eigenvalues);
1176
 
    throw(Squeal(Squeal::recoverable,
 
1176
    throw(Exception(Exception::recoverable,
1177
1177
                 "Failed to calculate eigenvalue",
1178
1178
                 "MAUS::eigenvalues"));
1179
1179
  }
1221
1221
  size_t rows = matrix.number_of_rows();
1222
1222
  size_t columns = matrix.number_of_columns();
1223
1223
  if (rows != columns) {
1224
 
    throw(Squeal(Squeal::recoverable,
 
1224
    throw(Exception(Exception::recoverable,
1225
1225
                 "Attempt to get eigensystem of non-square matrix",
1226
1226
                 "MAUS::eigensystem"));
1227
1227
  }
1236
1236
  if (ierr != 0) {
1237
1237
    gsl_vector_complex_free(eigenvalues);
1238
1238
    gsl_matrix_complex_free(eigenvectors);
1239
 
    throw(Squeal(Squeal::recoverable,
 
1239
    throw(Exception(Exception::recoverable,
1240
1240
                 "Failed to calculate eigenvalue",
1241
1241
                 "MAUS::eigenvectors"));
1242
1242
  }