GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GMatrixBase.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GMatrixBase.cpp - Abstract matrix base class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2006-2013 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GMatrixBase.cpp
23  * @brief Abstract matrix base class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #include <cmath>
29 #include "GTools.hpp"
30 #include "GException.hpp"
31 #include "GVector.hpp"
32 #include "GMatrixBase.hpp"
33 
34 /* __ Method name definitions ____________________________________________ */
35 #define G_COPY_MEMBERS "GMatrixBase::copy_members(GMatrixBase)"
36 #define G_SELECT_NON_ZERO "GMatrixBase::select_non_zero()"
37 
38 /* __ Debugging definitions ______________________________________________ */
39 //#define G_DEBUG_COPY_MEMBERS // Dump copy_members info
40 
41 
42 /*==========================================================================
43  = =
44  = Constructors/destructors =
45  = =
46  ==========================================================================*/
47 
48 /***********************************************************************//**
49  * @brief Void constructor
50  ***************************************************************************/
52 {
53  // Initialise members
54  init_members();
55 
56  // Return
57  return;
58 }
59 
60 
61 /***********************************************************************//**
62  * @brief Copy constructor
63  *
64  * @param[in] matrix Matrix.
65  *
66  * Constructs matrix by copying information from another matrix. The
67  * constructor is sufficiently generic to provide the base constructor for
68  * all derived classes, including sparse matrices.
69  ***************************************************************************/
71 {
72  // Initialise class members for clean destruction
73  init_members();
74 
75  // Copy members
76  copy_members(matrix);
77 
78  // Return
79  return;
80 }
81 
82 
83 /***********************************************************************//**
84  * @brief Destructor
85  ***************************************************************************/
87 {
88  // Free class members
89  free_members();
90 
91  // Return
92  return;
93 }
94 
95 
96 /*==========================================================================
97  = =
98  = Operators =
99  = =
100  ==========================================================================*/
101 
102 /***********************************************************************//**
103  * @brief Assignment operator
104  *
105  * @param[in] matrix Matrix.
106  * @return Matrix.
107  ***************************************************************************/
109 {
110  // Execute only if object is not identical
111  if (this != &matrix) {
112 
113  // Free members
114  free_members();
115 
116  // Initialise members
117  init_members();
118 
119  // Copy members
120  copy_members(matrix);
121 
122  } // endif: object was not identical
123 
124  // Return this object
125  return *this;
126 }
127 
128 
129 /***********************************************************************//**
130  * @brief Equalty operator
131  *
132  * @param[in] matrix Matrix.
133  * @return True if both matrices are identical.
134  *
135  * Checks if two matrices are identical. Two matrices are considered
136  * identical if they have the same dimensions and identical elements.
137  ***************************************************************************/
138 bool GMatrixBase::operator==(const GMatrixBase& matrix) const
139 {
140  // Initalise result to true (are identical)
141  bool result = true;
142 
143  // Test for difference. The loop will stop once then first difference
144  // is encountered.
145  if (m_rows == matrix.m_rows &&
146  m_cols == matrix.m_cols &&
147  m_elements == matrix.m_elements) {
148  for (int i = 0; i < m_elements; ++i) {
149  if (m_data[i] != matrix.m_data[i]) {
150  result = false;
151  break;
152  }
153  }
154  }
155  else {
156  result = false;
157  }
158 
159  // Return result
160  return result;
161 }
162 
163 
164 /***********************************************************************//**
165  * @brief Non-equality operator
166  *
167  * @param[in] matrix Matrix.
168  * @return True if both matrices are different.
169  *
170  * Checks if two matrices are different. Two matrices are considered
171  * different if they either have different dimensions or at least one element
172  * that differs.
173  ***************************************************************************/
174 bool GMatrixBase::operator!=(const GMatrixBase& matrix) const
175 {
176  // Get negated result of equality operation
177  bool result = !(this->operator==(matrix));
178 
179  // Return result
180  return result;
181 }
182 
183 
184 /*==========================================================================
185  = =
186  = Protected methods =
187  = =
188  ==========================================================================*/
189 
190 /***********************************************************************//**
191  * @brief Initialise class members
192  ***************************************************************************/
194 {
195  // Initialise members
196  m_rows = 0; // Number of rows
197  m_cols = 0; // Number of columns
198  m_elements = 0; // Logically used number of elements
199  m_alloc = 0; // Allocated # of elements (>= m_elements)
200  m_num_rowsel = 0; // Number of selected rows (for comp. decomp.)
201  m_num_colsel = 0; // Number of selected columns (for comp. decomp.)
202  m_data = NULL; // Matrix data
203  m_colstart = NULL; // Column start indices (m_cols+1)
204  m_rowsel = NULL; // Row selection (for compressed decomposition)
205  m_colsel = NULL; // Column selection (for compressed decomposition)
206 
207  // Return
208  return;
209 }
210 
211 
212 /***********************************************************************//**
213  * @brief Copy class members
214  *
215  * @param[in] matrix Matrix.
216  ***************************************************************************/
218 {
219  // Copy matrix attributes
220  m_rows = matrix.m_rows;
221  m_cols = matrix.m_cols;
222  m_elements = matrix.m_elements;
223  m_alloc = matrix.m_alloc;
224  m_num_rowsel = matrix.m_num_rowsel;
225  m_num_colsel = matrix.m_num_colsel;
226 
227  // Allocate only memory if we have rows and columns and data to copy
228  if (m_rows > 0 && m_cols > 0) {
229 
230  // Allocate memory for column start array and copy content
231  if (matrix.m_colstart != NULL) {
232  m_colstart = new int[m_cols+1];
233  for (int i = 0; i <= m_cols; ++i) {
234  m_colstart[i] = matrix.m_colstart[i];
235  }
236  }
237 
238  // Allocate memory for elements and copy them
239  if (matrix.m_data != NULL && m_alloc > 0) {
240  m_data = new double[m_alloc];
241  for (int i = 0; i < m_elements; ++i) {
242  m_data[i] = matrix.m_data[i];
243  }
244  }
245 
246  // If there is a row selection then copy it
247  if (matrix.m_rowsel != NULL && m_num_rowsel > 0) {
248  m_rowsel = new int[m_num_rowsel];
249  for (int i = 0; i < m_num_rowsel; ++i) {
250  m_rowsel[i] = matrix.m_rowsel[i];
251  }
252  }
253 
254  // If there is a column selection then copy it
255  if (matrix.m_colsel != NULL && m_num_colsel > 0) {
256  m_colsel = new int[m_num_colsel];
257  for (int i = 0; i < m_num_colsel; ++i) {
258  m_colsel[i] = matrix.m_colsel[i];
259  }
260  }
261 
262  } // endif: there were data to copy
263 
264  // Optionally show debug information
265  #if defined(G_DEBUG_COPY_MEMBERS)
266  std::cout << "GMatrixBase::copy_members:"
267  << " m_colstart=" << matrix.m_colstart << "->" << m_colstart
268  << " m_data=" << matrix.m_data << "->" << m_data
269  << " m_rowsel=" << matrix.m_rowsel << "->" << m_rowsel
270  << " m_colsel=" << matrix.m_colsel << "->" << m_colsel
271  << std::endl;
272  #endif
273 
274  // Return
275  return;
276 }
277 
278 
279 /***********************************************************************//**
280  * @brief Delete class members
281  ***************************************************************************/
283 {
284  // De-allocate only if memory has indeed been allocated
285  if (m_colsel != NULL) delete [] m_colsel;
286  if (m_rowsel != NULL) delete [] m_rowsel;
287  if (m_data != NULL) delete [] m_data;
288  if (m_colstart != NULL) delete [] m_colstart;
289 
290  // Properly mark members as free
291  m_colstart = NULL;
292  m_data = NULL;
293  m_rowsel = NULL;
294  m_colsel = NULL;
295 
296  // Return
297  return;
298 }
299 
300 
301 /***********************************************************************//**
302  * @brief Setup compressed matrix factorisation
303  *
304  * Determines the non-zero rows and columns in matrix and set up index
305  * arrays that point to these rows/columns. These arrays are used for
306  * compressed matrix factorisations.
307  ***************************************************************************/
309 {
310  // Free existing selection arrays
311  if (m_rowsel != NULL) delete [] m_rowsel;
312  if (m_colsel != NULL) delete [] m_colsel;
313 
314  // Allocate selection arrays
315  m_rowsel = new int[m_rows];
316  m_colsel = new int[m_cols];
317 
318  // Initialise non-zero row and column counters
319  m_num_rowsel = 0;
320  m_num_colsel = 0;
321 
322  // Declare loop variables
323  int row;
324  int col;
325 
326  // Find all non-zero rows
327  for (row = 0; row < m_rows; ++row) {
328  for (col = 0; col < m_cols; ++col) {
329  if ((*this)(row,col) != 0.0) {
330  break;
331  }
332  }
333  // Found a non-zero element in row
334  if (col < m_cols) {
336  }
337  }
338 
339  // Find all non-zero columns
340  for (col = 0; col < m_cols; ++col) {
341  for (row = 0; row < m_rows; ++row) {
342  if ((*this)(row,col) != 0.0)
343  break;
344  }
345  // Found a non-zero element in column
346  if (row < m_rows) {
347  m_colsel[m_num_colsel++] = col;
348  }
349  }
350 
351  // Return
352  return;
353 }
354 
355 
356 
357 /***********************************************************************//**
358  * @brief Scale all matrix elements with a scalar
359  *
360  * @param[in] scalar Scalar.
361  *
362  * Multiply all matrix elements with a scalar. There are three cases:
363  * - the multiplier is 0, then reset all elements to 0,
364  * - the multiplier is +/-1, then do nothing or negate,
365  * - in any other case, multiply by multiplier.
366  ***************************************************************************/
367 void GMatrixBase::scale_elements(const double& scalar)
368 {
369  // Case 1: If multiplicator is 0 then set entire matrix to 0
370  if (scalar == 0.0) {
371  for (int i = 0; i < m_elements; ++i) {
372  m_data[i] = 0.0;
373  }
374  }
375 
376  // Case 3: If multiplicator is not +/- 1 then perform multiplication
377  else if (std::abs(scalar) != 1.0) {
378  for (int i = 0; i < m_elements; ++i) {
379  m_data[i] *= scalar;
380  }
381  }
382 
383  // Case 2: If multiplication is -1 then negate
384  else if (scalar == -1.0) {
385  for (int i = 0; i < m_elements; ++i) {
386  m_data[i] = -m_data[i];
387  }
388  }
389 
390  // Return
391  return;
392 }
393 
394 
395 /***********************************************************************//**
396  * @brief Set all elements to a specific value
397  *
398  * @param[in] value Value.
399  *
400  * Sets all matrix elements to a specific value.
401  ***************************************************************************/
402 void GMatrixBase::set_all_elements(const double& value)
403 {
404  // Set all matrix elements
405  for (int i = 0; i < m_elements; ++i) {
406  m_data[i] = value;
407  }
408 
409  // Return
410  return;
411 }
412 
413 
414 /***********************************************************************//**
415  * @brief Return minimum matrix element
416  *
417  * @return Minimum element in matrix.
418  *
419  * Returns the minimum matrix element. If the matrix is empty, returns 0.
420  ***************************************************************************/
422 {
423  // Initialise result
424  double result = 0.0;
425 
426  // Continue only if there are elements
427  if (m_elements > 0) {
428 
429  // Set actual minimum to first elements
430  result = m_data[0];
431 
432  // Search all elements for the smallest one
433  for (int i = 1; i < m_elements; ++i) {
434  if (m_data[i] < result) {
435  result = m_data[i];
436  }
437  }
438 
439  } // endif: there were elements
440 
441  // Return result
442  return result;
443 }
444 
445 
446 /***********************************************************************//**
447  * @brief Returns maximum matrix element
448  *
449  * @return Maximum element in matrix.
450  *
451  * Returns the maximum matrix element. If the matrix is empty, returns 0.
452  ***************************************************************************/
454 {
455  // Initialise result
456  double result = 0.0;
457 
458  // Continue only if there are elements
459  if (m_elements > 0) {
460 
461  // Set actual maximum to first elements
462  result = m_data[0];
463 
464  // Search all elements for the largest one
465  for (int i = 1; i < m_elements; ++i) {
466  if (m_data[i] > result) {
467  result = m_data[i];
468  }
469  }
470 
471  } // endif: there were elements
472 
473  // Return result
474  return result;
475 }
476 
477 
478 /***********************************************************************//**
479  * @brief Returns sum over all matrix elements
480  ***************************************************************************/
482 {
483  // Initialise matrix sum
484  double result = 0.0;
485 
486  // Add all elements
487  for (int i = 0; i < m_elements; ++i) {
488  result += m_data[i];
489  }
490 
491  // Return result
492  return result;
493 }
494 
495 
496 /***********************************************************************//**
497  * @brief Print all matrix elements
498  *
499  * @param[in] chatter Chattiness (defaults to NORMAL).
500  * @param[in] num Maximum number of rows and columns to print (default: 10)
501  *
502  * Prints all matrix elements into a string. The parameter max_elements
503  * allows to control the maximum number of matrix elements that should be
504  * printed. If set to 0, all elements will be printed. Otherwise, the number
505  * of rows and columns will be limited by ommitting the central values.
506  ***************************************************************************/
507 std::string GMatrixBase::print_elements(const GChatter& chatter,
508  const int& num) const
509 {
510  // Initialise result string
511  std::string result;
512 
513  // Set row and column limits
514  int row_stop = 0;
515  int row_start = 0;
516  int col_stop = 0;
517  int col_start = 0;
518  if (num > 0) {
519  if (m_rows > num) {
520  row_stop = num/2;
521  row_start = m_rows - row_stop;
522  }
523  if (m_cols > num) {
524  col_stop = num/2;
525  col_start = m_cols - col_stop;
526  }
527  }
528 
529  // Print matrix elements row by row using the access function
530  for (int row = 0; row < row_stop; ++row) {
531  result += "\n ";
532  for (int col = 0; col < col_stop; ++col) {
533  result += gammalib::str((*this)(row,col));
534  if (col != m_cols-1) {
535  result += ", ";
536  }
537  }
538  if (col_start > col_stop) {
539  result += "... ";
540  }
541  for (int col = col_start; col < m_cols; ++col) {
542  result += gammalib::str((*this)(row,col));
543  if (col != m_cols-1) {
544  result += ", ";
545  }
546  }
547  }
548  if (row_start > row_stop) {
549  result += "\n ";
550  for (int col = 0; col < col_stop; ++col) {
551  result += "... ";
552  }
553  if (col_start > col_stop) {
554  result += "... ";
555  }
556  for (int col = col_start; col < m_cols; ++col) {
557  result += "... ";
558  }
559  }
560  for (int row = row_start; row < m_rows; ++row) {
561  result += "\n ";
562  for (int col = 0; col < col_stop; ++col) {
563  result += gammalib::str((*this)(row,col));
564  if (col != m_cols-1) {
565  result += ", ";
566  }
567  }
568  if (col_start > col_stop) {
569  result += "... ";
570  }
571  for (int col = col_start; col < m_cols; ++col) {
572  result += gammalib::str((*this)(row,col));
573  if (col != m_cols-1) {
574  result += ", ";
575  }
576  }
577  }
578 
579  // Return result
580  return result;
581 }
582 
583 
584 /***********************************************************************//**
585  * @brief Print row compression scheme if it exists
586  *
587  * @param[in] chatter Chattiness (defaults to NORMAL).
588  ***************************************************************************/
589 std::string GMatrixBase::print_row_compression(const GChatter& chatter) const
590 {
591  // Initialise result string
592  std::string result;
593 
594  // If there is a row compression the print the scheme
595  if (m_rowsel != NULL) {
596  result.append("\n"+gammalib::parformat("Row selection"));
597  for (int row = 0; row < m_num_rowsel; ++row) {
598  result.append(" ");
599  result.append(gammalib::str(m_rowsel[row]));
600  }
601  }
602 
603  // Return result
604  return result;
605 }
606 
607 
608 /***********************************************************************//**
609  * @brief Print column compression scheme if it exists
610  *
611  * @param[in] chatter Chattiness (defaults to NORMAL).
612  ***************************************************************************/
613 std::string GMatrixBase::print_col_compression(const GChatter& chatter) const
614 {
615  // Initialise result string
616  std::string result;
617 
618  // If there is a row compression the print the scheme
619  if (m_colsel != NULL) {
620  result.append("\n"+gammalib::parformat("Column selection"));
621  for (int col = 0; col < m_num_colsel; ++col) {
622  result.append(" ");
623  result.append(gammalib::str(m_colsel[col]));
624  }
625  }
626 
627  // Return result
628  return result;
629 }
double get_min_element(void) const
Return minimum matrix element.
void scale_elements(const double &scalar)
Scale all matrix elements with a scalar.
Abstract matrix base class interface definition.
std::string print_row_compression(const GChatter &chatter=NORMAL) const
Print row compression scheme if it exists.
void init_members(void)
Initialise class members.
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
double get_element_sum(void) const
Returns sum over all matrix elements.
virtual GVector row(const int &row) const =0
int m_rows
Number of rows.
Gammalib tools definition.
int m_num_colsel
Number of selected columns (for compressed decomposition)
void set_all_elements(const double &value)
Set all elements to a specific value.
std::string print_elements(const GChatter &chatter=NORMAL, const int &num=15) const
Print all matrix elements.
double * m_data
Matrix data.
Abstract matrix base class definition.
void copy_members(const GMatrixBase &matrix)
Copy class members.
int m_cols
Number of columns.
double get_max_element(void) const
Returns maximum matrix element.
GChatter
Definition: GTypemaps.hpp:33
int m_alloc
Size of allocated matrix memory.
virtual bool operator==(const GMatrixBase &matrix) const
Equalty operator.
std::string print_col_compression(const GChatter &chatter=NORMAL) const
Print column compression scheme if it exists.
Vector class interface definition.
virtual ~GMatrixBase(void)
Destructor.
Definition: GMatrixBase.cpp:86
int m_elements
Number of elements stored in matrix.
GMatrixBase(void)
Void constructor.
Definition: GMatrixBase.cpp:51
virtual bool operator!=(const GMatrixBase &matrix) const
Non-equality operator.
Exception handler interface definition.
virtual GMatrixBase & operator=(const GMatrixBase &matrix)
Assignment operator.
void free_members(void)
Delete class members.
int m_num_rowsel
Number of selected rows (for compressed decomposition)
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
int * m_colstart
Column start indices (m_cols+1)
int * m_rowsel
Row selection (for compressed decomposition)
void select_non_zero(void)
Setup compressed matrix factorisation.
int * m_colsel
Column selection (for compressed decomposition)
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413