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