GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GNdarray.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GNdarray.cpp - N-dimensional array class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2016-2018 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 GNdarray.cpp
23  * @brief N-dimensional array class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GTools.hpp"
32 #include "GException.hpp"
33 #include "GNdarray.hpp"
34 
35 
36 /* __ Method name definitions ____________________________________________ */
37 #define G_OP_ADD "GNdarray::operator+=(GNdarray&)"
38 #define G_OP_SUB "GNdarray::operator-=(GNdarray&)"
39 #define G_OP_MUL "GNdarray::operator*=(GNdarray&)"
40 #define G_OP_DIV "GNdarray::operator/=(GNdarray&)"
41 #define G_SHAPE "GNdarray::shape(std::vector<int>&)"
42 #define G_AT1 "GNdarray::at(int&)"
43 #define G_AT2 "GNdarray::at(int&, int&)"
44 #define G_AT3 "GNdarray::at(int&, int&, int&)"
45 #define G_ATN "GNdarray::at(std::vector<int>&)"
46 
47 
48 /*==========================================================================
49  = =
50  = Constructors/destructors =
51  = =
52  ==========================================================================*/
53 
54 /***********************************************************************//**
55  * @brief Void constructor
56  ***************************************************************************/
58 {
59  // Initialise class members
60  init_members();
61 
62  // Return
63  return;
64 }
65 
66 
67 /***********************************************************************//**
68  * @brief 1-dimensional constructor
69  *
70  * @param[in] nx Array dimension.
71  *
72  * Constructs a 1-dimensional array.
73  ***************************************************************************/
74 GNdarray::GNdarray(const int& nx)
75 {
76  // Initialise class members
77  init_members();
78 
79  // Set shape
80  m_shape.push_back(nx);
81 
82  // Set strides
83  m_strides.push_back(1);
84 
85  // Initialise data array with zeros
86  m_data.assign(nx, 0.0);
87 
88  // Return
89  return;
90 }
91 
92 
93 /***********************************************************************//**
94  * @brief 2-dimensional constructor
95  *
96  * @param[in] nx Array's first dimension.
97  * @param[in] ny Array's second dimension.
98  *
99  * Constructs a 2-dimensional array.
100  ***************************************************************************/
101 GNdarray::GNdarray(const int& nx, const int& ny)
102 {
103  // Initialise class members
104  init_members();
105 
106  // Compute size of array
107  int size = nx * ny;
108 
109  // Set shape
110  m_shape.push_back(nx);
111  m_shape.push_back(ny);
112 
113  // Set strides
114  m_strides.push_back(1);
115  m_strides.push_back(nx);
116 
117  // Initialise data array with zeros
118  m_data.assign(size, 0.0);
119 
120  // Return
121  return;
122 }
123 
124 
125 /***********************************************************************//**
126  * @brief 3-dimensional constructor
127  *
128  * @param[in] nx Array's first dimension.
129  * @param[in] ny Array's second dimension.
130  * @param[in] nz Array's third dimension.
131  *
132  * Constructs a 3-dimensional array.
133  ***************************************************************************/
134 GNdarray::GNdarray(const int& nx, const int& ny, const int& nz)
135 {
136  // Initialise class members
137  init_members();
138 
139  // Compute size of array
140  int size = nx * ny * nz;
141 
142  // Set shape
143  m_shape.push_back(nx);
144  m_shape.push_back(ny);
145  m_shape.push_back(nz);
146 
147  // Set strides
148  m_strides.push_back(1);
149  m_strides.push_back(nx);
150  m_strides.push_back(nx*ny);
151 
152  // Initialise data array with zeros
153  m_data.assign(size, 0.0);
154 
155  // Return
156  return;
157 }
158 
159 
160 /***********************************************************************//**
161  * @brief n-dimensional constructor
162  *
163  * @param[in] n Array dimensions.
164  *
165  * Constructs a n-dimensional array. If the array dimension vector @p n is
166  * empty an empty n-dimensional array is constructed (similar to the void
167  * constructor).
168  ***************************************************************************/
169 GNdarray::GNdarray(const std::vector<int>& n)
170 {
171  // Initialise class members
172  init_members();
173 
174  // Continue only if there are dimensions in vector
175  if (n.size() > 0) {
176 
177  // Compute size of array
178  int size = n[0];
179  for (int i = 1; i < n.size(); ++i) {
180  size *= n[i];
181  }
182 
183  // Set shape
184  m_shape = n;
185 
186  // Set strides
187  int stride = 1;
188  for (int i = 0; i < n.size(); ++i) {
189  m_strides.push_back(stride);
190  stride *= n[i];
191  }
192 
193  // Initialise data array with zeros
194  m_data.assign(size, 0.0);
195 
196  } // endif: there were dimensions in vector
197 
198  // Return
199  return;
200 }
201 
202 
203 /***********************************************************************//**
204  * @brief Copy constructor
205  *
206  * @param[in] array Array.
207  ***************************************************************************/
209 {
210  // Initialise class members
211  init_members();
212 
213  // Copy members
214  copy_members(array);
215 
216  // Return
217  return;
218 }
219 
220 
221 /***********************************************************************//**
222  * @brief Destructor
223  ***************************************************************************/
225 {
226  // Free members
227  free_members();
228 
229  // Return
230  return;
231 }
232 
233 
234 /*==========================================================================
235  = =
236  = Operators =
237  = =
238  ==========================================================================*/
239 
240 /***********************************************************************//**
241  * @brief Assignment operator
242  *
243  * @param[in] array Array.
244  * @return Array.
245  ***************************************************************************/
247 {
248  // Execute only if object is not identical
249  if (this != &array) {
250 
251  // Free members
252  free_members();
253 
254  // Initialise private members
255  init_members();
256 
257  // Copy members
258  copy_members(array);
259 
260  } // endif: object was not identical
261 
262  // Return this object
263  return *this;
264 }
265 
266 
267 /***********************************************************************//**
268  * @brief Equality operator
269  *
270  * @param[in] array Array.
271  * @return True if arrays are identical.
272  *
273  * Returns true if the arrays have identical shape and content.
274  ***************************************************************************/
275 bool GNdarray::operator==(const GNdarray& array) const
276 {
277  // Check of arrays have the same shape
278  bool identity = has_same_shape(array);
279 
280  // If the arrays have identical dimensions then check if all elements
281  // are identical. Break as soon as one array element differs.
282  if (identity) {
283  for (int i = 0; i < m_data.size(); ++i) {
284  if (m_data[i] != array.m_data[i]) {
285  identity = false;
286  break;
287  }
288  }
289  }
290 
291  // Return result
292  return identity;
293 }
294 
295 
296 /***********************************************************************//**
297  * @brief Non-equality operator
298  *
299  * @param[in] array Array.
300  * @return True if arrays are different.
301  *
302  * Returns true if the arrays have different shape or content.
303  ***************************************************************************/
304 bool GNdarray::operator!=(const GNdarray& array) const
305 {
306  // Get negated result of equality operation
307  bool difference = !(this->operator==(array));
308 
309  // Return result
310  return difference;
311 }
312 
313 
314 /***********************************************************************//**
315  * @brief Unary addition operator
316  *
317  * @param[in] array Array.
318  * @return Array.
319  *
320  * Adds an array to the current array.
321  ***************************************************************************/
323 {
324  // Throw an exception if the arrays have not the same shape
326 
327  // Add elements
328  for (int i = 0; i < m_data.size(); ++i) {
329  m_data[i] += array.m_data[i];
330  }
331 
332  // Return array
333  return *this;
334 }
335 
336 
337 /***********************************************************************//**
338  * @brief Unary multiplication operator
339  *
340  * @param[in] array Array.
341  * @return Array.
342  *
343  * Multiply an array by the current array.
344  ***************************************************************************/
346 {
347  // Throw an exception if the arrays have not the same shape
349 
350  // Multiply elements
351  for (int i = 0; i < m_data.size(); ++i) {
352  m_data[i] *= array.m_data[i];
353  }
354 
355  // Return array
356  return *this;
357 }
358 
359 
360 /***********************************************************************//**
361  * @brief Unary division operator
362  *
363  * @param[in] array Array.
364  * @return Array.
365  *
366  * Divide current array by array.
367  ***************************************************************************/
369 {
370  // Throw an exception if the arrays have not the same shape
372 
373  // Divide elements
374  for (int i = 0; i < m_data.size(); ++i) {
375 
376  // Only proceed if second array is non zero
377  if (array.m_data[i] != 0.0){
378  m_data[i] /= array.m_data[i];
379  }
380 
381  // ... otherwise set ratio to 0
382  else {
383  m_data[i] = 0.0;
384  }
385 
386  } // endfor: looped over array elements
387 
388  // Return array
389  return *this;
390 }
391 
392 
393 /***********************************************************************//**
394  * @brief Unary subtraction operator
395  *
396  * @param[in] array Array.
397  * @return Array.
398  *
399  * Subtracts an array from the current array.
400  ***************************************************************************/
402 {
403  // Throw an exception if the arrays have not the same shape
405 
406  // Subtract elements
407  for (int i = 0; i < m_data.size(); ++i) {
408  m_data[i] -= array.m_data[i];
409  }
410 
411  // Return array
412  return *this;
413 }
414 
415 
416 /***********************************************************************//**
417  * @brief Unary value addition operator
418  *
419  * @param[in] value Value.
420  * @return Array.
421  *
422  * Adds a value to all array elements.
423  ***************************************************************************/
424 GNdarray& GNdarray::operator+=(const double& value)
425 {
426  // Add value to elements
427  for (int i = 0; i < m_data.size(); ++i) {
428  m_data[i] += value;
429  }
430 
431  // Return array
432  return *this;
433 }
434 
435 
436 /***********************************************************************//**
437  * @brief Unary value subtraction operator
438  *
439  * @param[in] value Value.
440  * @return Array.
441  *
442  * Subtracts a value from all array elements.
443  ***************************************************************************/
444 GNdarray& GNdarray::operator-=(const double& value)
445 {
446  // Subtract value from elements
447  for (int i = 0; i < m_data.size(); ++i) {
448  m_data[i] -= value;
449  }
450 
451  // Return array
452  return *this;
453 }
454 
455 
456 /***********************************************************************//**
457  * @brief Unary multiplication operator
458  *
459  * @param[in] value Value.
460  * @return Array.
461  *
462  * Multiply all array elements by a value.
463  ***************************************************************************/
464 GNdarray& GNdarray::operator*=(const double& value)
465 {
466  // Multiply all elements
467  for (int i = 0; i < m_data.size(); ++i) {
468  m_data[i] *= value;
469  }
470 
471  // Return array
472  return *this;
473 }
474 
475 
476 /***********************************************************************//**
477  * @brief Unary division operator
478  *
479  * @param[in] value Value.
480  * @return Array.
481  *
482  * Divide all array elements by a value.
483  ***************************************************************************/
484 GNdarray& GNdarray::operator/=(const double& value)
485 {
486  // Divide all elements
487  for (int i = 0; i < m_data.size(); ++i) {
488  m_data[i] /= value;
489  }
490 
491  // Return array
492  return *this;
493 }
494 
495 
496 /***********************************************************************//**
497  * @brief Unary minus operator
498  *
499  * @return Array.
500  *
501  * Negate all array elements.
502  ***************************************************************************/
504 {
505  // Copy array
506  GNdarray result = *this;
507 
508  // Negate all elements
509  for (int i = 0; i < m_data.size(); ++i) {
510  result.m_data[i] = -result.m_data[i];
511  }
512 
513  // Return array
514  return result;
515 }
516 
517 
518 /*==========================================================================
519  = =
520  = Public methods =
521  = =
522  ==========================================================================*/
523 
524 /***********************************************************************//**
525  * @brief Clear array
526  ***************************************************************************/
527 void GNdarray::clear(void)
528 {
529  // Free members
530  free_members();
531 
532  // Initialise private members
533  init_members();
534 
535  // Return
536  return;
537 }
538 
539 
540 /***********************************************************************//**
541  * @brief Clone array
542  *
543  * @return Pointer to deep copy of array.
544  ***************************************************************************/
546 {
547  // Clone array
548  return new GNdarray(*this);
549 }
550 
551 
552 /***********************************************************************//**
553  * @brief Set shape of array
554  *
555  * @param[in] shape Shape vector.
556  *
557  * @exception GException::invalid_argument
558  * Invalid shape factorisation specified.
559  *
560  * Set the shape of the array. The shape specifies how the information is
561  * arranged in a n-dimensional array.
562  ***************************************************************************/
563 void GNdarray::shape(const std::vector<int>& shape)
564 {
565  // Computes the number of array elements
566  int nelements = 0;
567  if (shape.size() > 0) {
568  nelements = 1;
569  for (int i = 0; i < shape.size(); ++i) {
570  nelements *= shape[i];
571  }
572  }
573 
574  // Throw an exception if resulting number of elements is not equal to
575  // the existing number of elements
576  if (nelements != size()) {
577  std::string msg = "Number of elements "+gammalib::str(nelements)+
578  " in specified shape is not identical to the "
579  "number of elements "+gammalib::str(size())+
580  " in the array.";
582  }
583 
584  // Set shape
585  m_shape = shape;
586 
587  // Return
588  return;
589 }
590 
591 
592 /***********************************************************************//**
593  * @brief 1-dimensional array element access with range checking
594  * (const version)
595  *
596  * @param[in] ix Element index [0,...,shape(0)-1].
597  * @return Const reference to array element.
598  *
599  * @exception GException::invalid_value
600  * Array is not 1-dimensional.
601  * @exception GException::out_of_range
602  * Element index is out of range.
603  ***************************************************************************/
604 const double& GNdarray::at(const int& ix) const
605 {
606  // Throw an exception if the array is not 1-dimensional
607  if (m_shape.size() != 1) {
608  std::string msg = "Invalid access of "+
609  gammalib::str(m_shape.size())+"-dimensional array "
610  "with 1-dimensional access operator.";
611  throw GException::invalid_value(G_AT1, msg);
612  }
613 
614  // Throw an exception if the index is outside the valid range
615  if (ix < 0 || ix >= m_shape[0]) {
616  throw GException::out_of_range(G_AT1, "Array index", ix, m_shape[0]);
617  }
618 
619  // Return array element
620  return (*this)(ix);
621 }
622 
623 
624 /***********************************************************************//**
625  * @brief 2-dimensional array element access with range checking
626  * (const version)
627  *
628  * @param[in] ix Index in first dimension [0,...,shape(0)-1].
629  * @param[in] iy Index in second dimension [0,...,shape(1)-1].
630  * @return Const reference to array element.
631  *
632  * @exception GException::invalid_value
633  * Array is not 2-dimensional.
634  * @exception GException::out_of_range
635  * Element index is out of range.
636  ***************************************************************************/
637 const double& GNdarray::at(const int& ix, const int& iy) const
638 {
639  // Throw an exception if the array is not 2-dimensional
640  if (m_shape.size() != 2) {
641  std::string msg = "Invalid access of "+
642  gammalib::str(m_shape.size())+"-dimensional array "
643  "with 2-dimensional access operator.";
644  throw GException::invalid_value(G_AT2, msg);
645  }
646 
647  // Throw an exception if the indices are outside the valid range
648  if (ix < 0 || ix >= m_shape[0]) {
649  throw GException::out_of_range(G_AT2, "First array index", ix,
650  m_shape[0]);
651  }
652  if (iy < 0 || iy >= m_shape[1]) {
653  throw GException::out_of_range(G_AT2, "Second array index", iy,
654  m_shape[1]);
655  }
656 
657  // Return array element
658  return (*this)(ix,iy);
659 }
660 
661 
662 /***********************************************************************//**
663  * @brief 3-dimensional array element access with range checking
664  * (const version)
665  *
666  * @param[in] ix Index in first dimension [0,...,shape(0)-1].
667  * @param[in] iy Index in second dimension [0,...,shape(1)-1].
668  * @param[in] iz Index in third dimension [0,...,shape(2)-1].
669  * @return Const reference to array element.
670  *
671  * @exception GException::invalid_value
672  * Array is not 3-dimensional.
673  * @exception GException::out_of_range
674  * Element index is out of range.
675  ***************************************************************************/
676 const double& GNdarray::at(const int& ix, const int& iy, const int& iz) const
677 {
678  // Throw an exception if the array is not 2-dimensional
679  if (m_shape.size() != 3) {
680  std::string msg = "Invalid access of "+
681  gammalib::str(m_shape.size())+"-dimensional array "
682  "with 3-dimensional access operator.";
683  throw GException::invalid_value(G_AT3, msg);
684  }
685 
686  // Throw an exception if the indices are outside the valid range
687  if (ix < 0 || ix >= m_shape[0]) {
688  throw GException::out_of_range(G_AT3, "First array index", ix,
689  m_shape[0]);
690  }
691  if (iy < 0 || iy >= m_shape[1]) {
692  throw GException::out_of_range(G_AT3, "Second array index", iy,
693  m_shape[1]);
694  }
695  if (iz < 0 || iz >= m_shape[2]) {
696  throw GException::out_of_range(G_AT3, "Third array index", iz,
697  m_shape[2]);
698  }
699 
700  // Return array element
701  return (*this)(ix,iy,iz);
702 }
703 
704 
705 /***********************************************************************//**
706  * @brief n-dimensional array element access with range checking
707  * (const version)
708  *
709  * @param[in] i Index vector.
710  * @return Const reference to array element.
711  *
712  * @exception GException::invalid_value
713  * Array is not of correct dimension.
714  * @exception GException::out_of_range
715  * Element index is out of range.
716  ***************************************************************************/
717 const double& GNdarray::at(const std::vector<int>& i) const
718 {
719  // Throw an exception if the array is not 2-dimensional
720  if (m_shape.size() != i.size()) {
721  std::string msg = "Invalid access of "+
722  gammalib::str(m_shape.size())+"-dimensional array "
723  "with "+
724  gammalib::str(i.size())+"-dimensional access "
725  "operator.";
726  throw GException::invalid_value(G_ATN, msg);
727  }
728 
729  // Throw an exception if the indices are outside the valid range
730  for (int k = 0; k < m_shape.size(); ++k) {
731  if (i[k] < 0 || i[k] >= m_shape[k]) {
732  throw GException::out_of_range(G_ATN, "Dimension "+
733  gammalib::str(k)+" index",
734  i[k], m_shape[k]);
735  }
736  }
737 
738  // Return array element
739  return (*this)(i);
740 }
741 
742 
743 /***********************************************************************//**
744  * @brief Print array information
745  *
746  * @param[in] chatter Chattiness.
747  * @return String containing array information.
748  ***************************************************************************/
749 std::string GNdarray::print(const GChatter& chatter) const
750 {
751  // Initialise result string
752  std::string result = "(";
753 
754  // Continue only if chatter is not silent
755  if (chatter != SILENT) {
756 
757  // Put all elements in string
758  for (int i = 0; i < m_data.size(); ++i) {
759  if (i > 0) {
760  result += ", ";
761  }
762  result += gammalib::str(m_data[i]);
763  }
764 
765  // Append closing parentheses
766  result += ")";
767 
768  } // endif: chatter was not silent
769 
770  // Return result
771  return result;
772 }
773 
774 
775 /*==========================================================================
776  = =
777  = Private methods =
778  = =
779  ==========================================================================*/
780 
781 /***********************************************************************//**
782  * @brief Initialise class members
783  ***************************************************************************/
785 {
786  // Initialise members
787  m_shape.clear();
788  m_strides.clear();
789  m_data.clear();
790 
791  // Actually reallocate the memory of m_data
792  std::vector<double>(m_data).swap(m_data);
793 
794  // Return
795  return;
796 }
797 
798 
799 /***********************************************************************//**
800  * @brief Copy class members
801  *
802  * @param[in] array Array.
803  ***************************************************************************/
805 {
806  // Copy members
807  m_shape = array.m_shape;
808  m_strides = array.m_strides;
809  m_data = array.m_data;
810 
811  // Return
812  return;
813 }
814 
815 
816 /***********************************************************************//**
817  * @brief Delete class members
818  ***************************************************************************/
820 {
821  // Return
822  return;
823 }
824 
825 
826 /***********************************************************************//**
827  * @brief Compute array element index
828  *
829  * @param[in] i Index vector.
830  * @return Array element index.
831  *
832  * Computes the array element index from an index vector. The method assumes
833  * that the index vector has the same dimension as the array. If this
834  * condition is not guaranteed, the condition needs to be checked before
835  * calling the method.
836  ***************************************************************************/
837 int GNdarray::index(const std::vector<int>& i) const
838 {
839  // Initialise element index
840  int index = 0;
841 
842  // Compute element index
843  for (int k = 0; k < m_shape.size(); ++k) {
844  index += m_strides[k] * i[k];
845  }
846 
847  // Return index
848  return index;
849 }
850 
851 
852 /***********************************************************************//**
853  * @brief Check if array has the same shape
854  *
855  * @param[in] array Array.
856  * @return True if the array has the same shape.
857  ***************************************************************************/
858 bool GNdarray::has_same_shape(const GNdarray& array) const
859 {
860  // Check if the array has the same dimension
861  bool identity = m_shape.size() == array.m_shape.size();
862 
863  // Check if the shape is identical. Break as soon as one shape value
864  // differs.
865  if (identity) {
866  for (int i = 0; i < m_shape.size(); ++i) {
867  if (m_shape[i] != array.m_shape[i]) {
868  identity = false;
869  break;
870  }
871  }
872  }
873 
874  // Return result
875  return identity;
876 }
877 
878 
879 /***********************************************************************//**
880  * @brief Throw exception if array shapes differ
881  *
882  * @param[in] method Method that throws exception.
883  * @param[in] array Array.
884  ***************************************************************************/
885 void GNdarray::require_same_shape(const std::string& method,
886  const GNdarray& array) const
887 {
888  // If the shape differs then throw an exception
889  if (!has_same_shape(array)) {
890 
891  // Compose exception message
892  std::string msg = "Incompatible array dimensions (";
893  for (int i = 0; i < m_shape.size(); ++i) {
894  if (i > 0) {
895  msg += ", ";
896  }
897  msg += gammalib::str(m_shape[i]);
898  }
899  msg += ") <=> (";
900  for (int i = 0; i < array.m_shape.size(); ++i) {
901  if (i > 0) {
902  msg += ", ";
903  }
904  msg += gammalib::str(array.m_shape[i]);
905  }
906  msg += ").";
907 
908  // Throw exception
909  throw GException::invalid_argument(method, msg);
910 
911  } // endif: array shapes differed
912 
913  // Return
914  return;
915 }
916 
917 
918 /*==========================================================================
919  = =
920  = Friends =
921  = =
922  ==========================================================================*/
923 
924 /***********************************************************************//**
925  * @brief Computes minimum array element
926  *
927  * @param[in] array Array.
928  * @return Minimum array element.
929  ***************************************************************************/
930 double min(const GNdarray& array)
931 {
932  // Initialises result
933  double result = 0.0;
934 
935  // Continue only if we have elements
936  if (array.m_data.size() > 0) {
937 
938  // Search for minimum
939  result = array.m_data[0];
940  for (int i = 1; i < array.m_data.size(); ++i) {
941  if (array.m_data[i] < result) {
942  result = array.m_data[i];
943  }
944  }
945 
946  } // endif: there were elements
947 
948  // Returns minimum
949  return result;
950 }
951 
952 
953 /***********************************************************************//**
954  * @brief Computes maximum array element
955  *
956  * @param[in] array Array.
957  * @return Maximum array element.
958  ***************************************************************************/
959 double max(const GNdarray& array)
960 {
961  // Initialises result
962  double result = 0.0;
963 
964  // Continue only if we have elements
965  if (array.m_data.size() > 0) {
966 
967  // Search for maximum
968  result = array.m_data[0];
969  for (int i = 1; i < array.m_data.size(); ++i) {
970  if (array.m_data[i] > result) {
971  result = array.m_data[i];
972  }
973  }
974 
975  } // endif: there were elements
976 
977  // Returns maximum
978  return result;
979 }
980 
981 
982 /***********************************************************************//**
983  * @brief Computes array sum
984  *
985  * @param[in] array Array.
986  * @return Sum of array elements.
987  ***************************************************************************/
988 double sum(const GNdarray& array)
989 {
990  // Compute sum
991  double result = 0.0;
992  for (int i = 0; i < array.m_data.size(); ++i) {
993  result += array.m_data[i];
994  }
995 
996  // Returns sum
997  return result;
998 }
999 
1000 
1001 /***********************************************************************//**
1002  * @brief Computes arccos of array elements
1003  *
1004  * @param[in] array Array.
1005  * @return Array containing the arccos of every element.
1006  ***************************************************************************/
1007 GNdarray acos(const GNdarray& array)
1008 {
1009  // Initialise result array
1010  GNdarray result = array;
1011 
1012  // Evaluate each array element
1013  for (int i = 0; i < array.m_data.size(); ++i) {
1014  result.m_data[i] = std::acos(array.m_data[i]);
1015  }
1016 
1017  // Return array
1018  return result;
1019 }
1020 
1021 
1022 /***********************************************************************//**
1023  * @brief Computes acosh of array elements
1024  *
1025  * @param[in] array Array.
1026  * @return Array containing the acosh of every element.
1027  ***************************************************************************/
1028 GNdarray acosh(const GNdarray& array)
1029 {
1030  // Initialise result array
1031  GNdarray result = array;
1032 
1033  // Evaluate each array element
1034  for (int i = 0; i < array.m_data.size(); ++i) {
1035  result.m_data[i] = acosh(array.m_data[i]);
1036  }
1037 
1038  // Return array
1039  return result;
1040 }
1041 
1042 
1043 /***********************************************************************//**
1044  * @brief Computes arcsin of array elements
1045  *
1046  * @param[in] array Array.
1047  * @return Array containing the arcsin of every element.
1048  ***************************************************************************/
1049 GNdarray asin(const GNdarray& array)
1050 {
1051  // Initialise result array
1052  GNdarray result = array;
1053 
1054  // Evaluate each array element
1055  for (int i = 0; i < array.m_data.size(); ++i) {
1056  result.m_data[i] = std::asin(array.m_data[i]);
1057  }
1058 
1059  // Return array
1060  return result;
1061 }
1062 
1063 
1064 /***********************************************************************//**
1065  * @brief Computes asinh of array elements
1066  *
1067  * @param[in] array Array.
1068  * @return Array containing the asinh of every element.
1069  ***************************************************************************/
1070 GNdarray asinh(const GNdarray& array)
1071 {
1072  // Initialise result array
1073  GNdarray result = array;
1074 
1075  // Evaluate each array element
1076  for (int i = 0; i < array.m_data.size(); ++i) {
1077  result.m_data[i] = asinh(array.m_data[i]);
1078  }
1079 
1080  // Return array
1081  return result;
1082 }
1083 
1084 
1085 /***********************************************************************//**
1086  * @brief Computes arctan of array elements
1087  *
1088  * @param[in] array Array.
1089  * @return Array containing the arctan of every element.
1090  ***************************************************************************/
1091 GNdarray atan(const GNdarray& array)
1092 {
1093  // Initialise result array
1094  GNdarray result = array;
1095 
1096  // Evaluate each array element
1097  for (int i = 0; i < array.m_data.size(); ++i) {
1098  result.m_data[i] = std::atan(array.m_data[i]);
1099  }
1100 
1101  // Return array
1102  return result;
1103 }
1104 
1105 
1106 /***********************************************************************//**
1107  * @brief Computes atanh of array elements
1108  *
1109  * @param[in] array Array.
1110  * @return Array containing the atanh of every element.
1111  ***************************************************************************/
1112 GNdarray atanh(const GNdarray& array)
1113 {
1114  // Initialise result array
1115  GNdarray result = array;
1116 
1117  // Evaluate each array element
1118  for (int i = 0; i < array.m_data.size(); ++i) {
1119  result.m_data[i] = atanh(array.m_data[i]);
1120  }
1121 
1122  // Return array
1123  return result;
1124 }
1125 
1126 
1127 /***********************************************************************//**
1128  * @brief Computes cosine of array elements
1129  *
1130  * @param[in] array Array.
1131  * @return Array containing the cosine of every element.
1132  ***************************************************************************/
1133 GNdarray cos(const GNdarray& array)
1134 {
1135  // Initialise result array
1136  GNdarray result = array;
1137 
1138  // Evaluate each array element
1139  for (int i = 0; i < array.m_data.size(); ++i) {
1140  result.m_data[i] = std::cos(array.m_data[i]);
1141  }
1142 
1143  // Return array
1144  return result;
1145 }
1146 
1147 
1148 /***********************************************************************//**
1149  * @brief Computes cosh of array elements
1150  *
1151  * @param[in] array Array.
1152  * @return Array containing the cosh of every element.
1153  ***************************************************************************/
1154 GNdarray cosh(const GNdarray& array)
1155 {
1156  // Initialise result array
1157  GNdarray result = array;
1158 
1159  // Evaluate each array element
1160  for (int i = 0; i < array.m_data.size(); ++i) {
1161  result.m_data[i] = std::cosh(array.m_data[i]);
1162  }
1163 
1164  // Return array
1165  return result;
1166 }
1167 
1168 
1169 /***********************************************************************//**
1170  * @brief Computes exponential of array elements
1171  *
1172  * @param[in] array Array.
1173  * @return Array containing the exponential of every element.
1174  ***************************************************************************/
1175 GNdarray exp(const GNdarray& array)
1176 {
1177  // Initialise result array
1178  GNdarray result = array;
1179 
1180  // Evaluate each array element
1181  for (int i = 0; i < array.m_data.size(); ++i) {
1182  result.m_data[i] = std::exp(array.m_data[i]);
1183  }
1184 
1185  // Return array
1186  return result;
1187 }
1188 
1189 
1190 /***********************************************************************//**
1191  * @brief Computes absolute of array elements
1192  *
1193  * @param[in] array Array.
1194  * @return Array containing the absolute of every element.
1195  ***************************************************************************/
1196 GNdarray abs(const GNdarray& array)
1197 {
1198  // Initialise result array
1199  GNdarray result = array;
1200 
1201  // Evaluate each array element
1202  for (int i = 0; i < array.m_data.size(); ++i) {
1203  result.m_data[i] = std::abs(array.m_data[i]);
1204  }
1205 
1206  // Return array
1207  return result;
1208 }
1209 
1210 
1211 /***********************************************************************//**
1212  * @brief Computes natural logarithm of array elements
1213  *
1214  * @param[in] array Array.
1215  * @return Array containing the natural logarithm of every element.
1216  ***************************************************************************/
1217 GNdarray log(const GNdarray& array)
1218 {
1219  // Initialise result array
1220  GNdarray result = array;
1221 
1222  // Evaluate each array element
1223  for (int i = 0; i < array.m_data.size(); ++i) {
1224  result.m_data[i] = std::log(array.m_data[i]);
1225  }
1226 
1227  // Return array
1228  return result;
1229 }
1230 
1231 
1232 /***********************************************************************//**
1233  * @brief Computes sign of array elements
1234  *
1235  * @param[in] array Array.
1236  * @return Array containing the sign of every element.
1237  ***************************************************************************/
1238 GNdarray sign(const GNdarray& array)
1239 {
1240  // Initialise result array
1241  GNdarray result = array;
1242 
1243  // Evaluate each array element
1244  for (int i = 0; i < array.m_data.size(); ++i) {
1245 
1246  // Get array content
1247  double content = array.m_data[i];
1248 
1249  // Initialise sign
1250  double sign;
1251 
1252  // Handle the 3 cases
1253  if (content < 0.0) {
1254  sign = -1.0;
1255  }
1256  else if (content > 0.0) {
1257  sign = +1.0;
1258  }
1259  else {
1260  sign = 0.0;
1261  }
1262 
1263  // Store sign
1264  result.m_data[i] = sign;
1265 
1266  } // endfor: looped over all array elements
1267 
1268  // Return array
1269  return result;
1270 }
1271 
1272 
1273 /***********************************************************************//**
1274  * @brief Computes base10 logarithm of array elements
1275  *
1276  * @param[in] array Array.
1277  * @return Array containing the base10 logarithm of every element.
1278  ***************************************************************************/
1279 GNdarray log10(const GNdarray& array)
1280 {
1281  // Initialise result array
1282  GNdarray result = array;
1283 
1284  // Evaluate each array element
1285  for (int i = 0; i < array.m_data.size(); ++i) {
1286  result.m_data[i] = std::log10(array.m_data[i]);
1287  }
1288 
1289  // Return array
1290  return result;
1291 }
1292 
1293 
1294 /***********************************************************************//**
1295  * @brief Computes sine of array elements
1296  *
1297  * @param[in] array Array.
1298  * @return Array containing the sine of every element.
1299  ***************************************************************************/
1300 GNdarray sin(const GNdarray& array)
1301 {
1302  // Initialise result array
1303  GNdarray result = array;
1304 
1305  // Evaluate each array element
1306  for (int i = 0; i < array.m_data.size(); ++i) {
1307  result.m_data[i] = std::sin(array.m_data[i]);
1308  }
1309 
1310  // Return array
1311  return result;
1312 }
1313 
1314 
1315 /***********************************************************************//**
1316  * @brief Computes sinh of array elements
1317  *
1318  * @param[in] array Array.
1319  * @return Array containing the sinh of every element.
1320  ***************************************************************************/
1321 GNdarray sinh(const GNdarray& array)
1322 {
1323  // Initialise result array
1324  GNdarray result = array;
1325 
1326  // Evaluate each array element
1327  for (int i = 0; i < array.m_data.size(); ++i) {
1328  result.m_data[i] = std::sinh(array.m_data[i]);
1329  }
1330 
1331  // Return array
1332  return result;
1333 }
1334 
1335 
1336 /***********************************************************************//**
1337  * @brief Computes square root of array elements
1338  *
1339  * @param[in] array Array.
1340  * @return Array containing the square root of every element.
1341  ***************************************************************************/
1342 GNdarray sqrt(const GNdarray& array)
1343 {
1344  // Initialise result array
1345  GNdarray result = array;
1346 
1347  // Evaluate each array element
1348  for (int i = 0; i < array.m_data.size(); ++i) {
1349  result.m_data[i] = std::sqrt(array.m_data[i]);
1350  }
1351 
1352  // Return array
1353  return result;
1354 }
1355 
1356 
1357 /***********************************************************************//**
1358  * @brief Computes tangens of array elements
1359  *
1360  * @param[in] array Array.
1361  * @return Array containing the tangens of every element.
1362  ***************************************************************************/
1363 GNdarray tan(const GNdarray& array)
1364 {
1365  // Initialise result array
1366  GNdarray result = array;
1367 
1368  // Evaluate each array element
1369  for (int i = 0; i < array.m_data.size(); ++i) {
1370  result.m_data[i] = std::tan(array.m_data[i]);
1371  }
1372 
1373  // Return array
1374  return result;
1375 }
1376 
1377 
1378 /***********************************************************************//**
1379  * @brief Computes tanh of array elements
1380  *
1381  * @param[in] array Array.
1382  * @return Array containing the tanh of every element.
1383  ***************************************************************************/
1384 GNdarray tanh(const GNdarray& array)
1385 {
1386  // Initialise result array
1387  GNdarray result = array;
1388 
1389  // Evaluate each array element
1390  for (int i = 0; i < array.m_data.size(); ++i) {
1391  result.m_data[i] = std::tanh(array.m_data[i]);
1392  }
1393 
1394  // Return array
1395  return result;
1396 }
1397 
1398 
1399 /***********************************************************************//**
1400  * @brief Computes tanh of array elements
1401  *
1402  * @param[in] array Array.
1403  * @param[in] power Power.
1404  * @return Array containing the power of every element.
1405  ***************************************************************************/
1406 GNdarray pow(const GNdarray& array, const double& power)
1407 {
1408  // Initialise result array
1409  GNdarray result = array;
1410 
1411  // Evaluate each array element
1412  for (int i = 0; i < array.m_data.size(); ++i) {
1413  result.m_data[i] = std::pow(array.m_data[i], power);
1414  }
1415 
1416  // Return array
1417  return result;
1418 }
GVector acosh(const GVector &vector)
Computes acosh of vector elements.
Definition: GVector.cpp:995
GVector sinh(const GVector &vector)
Computes sinh of vector elements.
Definition: GVector.cpp:1247
bool operator!=(const GNdarray &array) const
Non-equality operator.
Definition: GNdarray.cpp:304
void clear(void)
Clear array.
Definition: GNdarray.cpp:527
#define G_AT2
Definition: GNdarray.cpp:43
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
#define G_OP_ADD
Definition: GNdarray.cpp:37
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
void free_members(void)
Delete class members.
Definition: GNdarray.cpp:819
GNdarray * clone(void) const
Clone array.
Definition: GNdarray.cpp:545
#define G_OP_DIV
Definition: GNdarray.cpp:40
double sum(const GVector &vector)
Computes vector sum.
Definition: GVector.cpp:901
double acos(const double &arg)
Computes acos by avoiding NaN due to rounding errors.
Definition: GMath.cpp:69
Gammalib tools definition.
const std::vector< int > & shape(void) const
Return shape of array.
Definition: GNdarray.hpp:307
double min(const GVector &vector)
Computes minimum vector element.
Definition: GVector.cpp:843
#define G_ATN
Definition: GNdarray.cpp:45
std::vector< int > m_shape
Array dimensions.
Definition: GNdarray.hpp:135
void copy_members(const GNdarray &array)
Copy class members.
Definition: GNdarray.cpp:804
#define G_AT3
Definition: GNdarray.cpp:44
#define G_AT1
Definition: GNdarray.cpp:42
double & at(const int &ix)
1-dimensional array element access with range checking
Definition: GNdarray.hpp:334
GNdarray & operator-=(const GNdarray &array)
Unary subtraction operator.
Definition: GNdarray.cpp:401
std::vector< double > m_data
Array data.
Definition: GNdarray.hpp:137
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
GNdarray & operator/=(const GNdarray &array)
Unary division operator.
Definition: GNdarray.cpp:368
GNdarray & operator=(const GNdarray &array)
Assignment operator.
Definition: GNdarray.cpp:246
void require_same_shape(const std::string &method, const GNdarray &array) const
Throw exception if array shapes differ.
Definition: GNdarray.cpp:885
GVector tan(const GVector &vector)
Computes tangens of vector elements.
Definition: GVector.cpp:1289
N-dimensional array class interface definition.
#define G_OP_SUB
Definition: GNdarray.cpp:38
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1184
GNdarray sign(const GNdarray &array)
Computes sign of array elements.
Definition: GNdarray.cpp:1238
std::string print(const GChatter &chatter=NORMAL) const
Print array information.
Definition: GNdarray.cpp:749
bool operator==(const GNdarray &array) const
Equality operator.
Definition: GNdarray.cpp:275
int size(void) const
Return number of elements in array.
Definition: GNdarray.hpp:293
GChatter
Definition: GTypemaps.hpp:33
GNdarray(void)
Void constructor.
Definition: GNdarray.cpp:57
N-dimensional array class.
Definition: GNdarray.hpp:44
double max(const GVector &vector)
Computes maximum vector element.
Definition: GVector.cpp:872
GVector asin(const GVector &vector)
Computes arcsin of vector elements.
Definition: GVector.cpp:1016
GNdarray & operator*=(const GNdarray &array)
Unary multiplication operator.
Definition: GNdarray.cpp:345
GVector tanh(const GVector &vector)
Computes tanh of vector elements.
Definition: GVector.cpp:1310
GVector atanh(const GVector &vector)
Computes atanh of vector elements.
Definition: GVector.cpp:1079
GVector cosh(const GVector &vector)
Computes cosh of vector elements.
Definition: GVector.cpp:1121
GVector asinh(const GVector &vector)
Computes asinh of vector elements.
Definition: GVector.cpp:1037
int index(const std::vector< int > &i) const
Compute array element index.
Definition: GNdarray.cpp:837
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1332
bool has_same_shape(const GNdarray &array) const
Check if array has the same shape.
Definition: GNdarray.cpp:858
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
Exception handler interface definition.
void init_members(void)
Initialise class members.
Definition: GNdarray.cpp:784
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
GNdarray & operator+=(const GNdarray &array)
Unary addition operator.
Definition: GNdarray.cpp:322
virtual ~GNdarray(void)
Destructor.
Definition: GNdarray.cpp:224
GVector atan(const GVector &vector)
Computes arctan of vector elements.
Definition: GVector.cpp:1058
#define G_OP_MUL
Definition: GNdarray.cpp:39
GNdarray operator-(void) const
Unary minus operator.
Definition: GNdarray.cpp:503
#define G_SHAPE
Definition: GNdarray.cpp:41
std::vector< int > m_strides
Steps in each dimension when traversing array.
Definition: GNdarray.hpp:136
GVector log10(const GVector &vector)
Computes base10 logarithm of vector elements.
Definition: GVector.cpp:1205
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413