From 381fb883a159164f722de042e4771b9fa6b4f08e Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Tue, 28 May 2019 10:55:39 +0200 Subject: [PATCH] Clean up of old vector and matrix documentation --- doc/source/users/tutorial/index.rst | 2 +- doc/source/users/tutorial/intro.rst | 25 +++++++++++++++++++++++++ doc/source/users/tutorial/linalg.rst | 763 ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- doc/source/users/tutorial/sky/convert_coordinates.rst | 24 +++++++++++------------- doc/source/users/tutorial/sky/convert_projections.rst | 21 +++++++++------------ doc/source/users/user_manual/modules/linalg/index.rst | 4 +++- doc/source/users/user_manual/modules/linalg/matrix.rst | 127 ------------------------------------------------------------------------------------------------------------------------------- doc/source/users/user_manual/modules/linalg/matrix_arithmetics.rst | 102 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/users/user_manual/modules/linalg/matrix_storage.rst | 248 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/users/user_manual/modules/linalg/vector.rst | 264 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 10 files changed, 663 insertions(+), 917 deletions(-) create mode 100644 doc/source/users/tutorial/intro.rst delete mode 100644 doc/source/users/tutorial/linalg.rst delete mode 100644 doc/source/users/user_manual/modules/linalg/matrix.rst create mode 100644 doc/source/users/user_manual/modules/linalg/matrix_arithmetics.rst create mode 100644 doc/source/users/user_manual/modules/linalg/matrix_storage.rst create mode 100644 doc/source/users/user_manual/modules/linalg/vector.rst diff --git a/doc/source/users/tutorial/index.rst b/doc/source/users/tutorial/index.rst index 1798311..ed0c4b0 100644 --- a/doc/source/users/tutorial/index.rst +++ b/doc/source/users/tutorial/index.rst @@ -4,6 +4,6 @@ Tutorials .. toctree:: :maxdepth: 1 + intro sky/index - linalg opt/index diff --git a/doc/source/users/tutorial/intro.rst b/doc/source/users/tutorial/intro.rst new file mode 100644 index 0000000..45a0019 --- /dev/null +++ b/doc/source/users/tutorial/intro.rst @@ -0,0 +1,25 @@ +Introduction +============ + +These tutorials provide some usage examples of the GammaLib C++ library and +its Python module. + +To test the C++ examples they need to be enclosed in the following code +structure + +.. code-block:: cpp + :linenos: + + #include "GammaLib.hpp" + int main(void) { + // example code goes here + return 0; + } + +For the Python examples the gammalib module needs to be imported + +.. code-block:: python + :linenos: + + import gammalib + # example code goes here diff --git a/doc/source/users/tutorial/linalg.rst b/doc/source/users/tutorial/linalg.rst deleted file mode 100644 index 0bd621d..0000000 --- a/doc/source/users/tutorial/linalg.rst +++ /dev/null @@ -1,763 +0,0 @@ -Linear algebra --------------- - -Vectors -~~~~~~~ - -General -^^^^^^^ - -The :doxy:`GVector` class implements a vector as a one-dimensional array of -successive double precision floating point values. You may allocate -an empty vector using - -.. code-block:: cpp - - GVector vector; - -and a vector of specific dimension using - -.. code-block:: cpp - - GVector vector(10); // Allocates a vector with 10 elements - -On allocation, all elements of a vector are set to 0. - -You may also allocate a vector by copying from another vector - -.. code-block:: cpp - - GVector vector(10); // Allocates a vector with 10 elements - GVector another = vector; // Allocates vector by copying - -or by using - -.. code-block:: cpp - - GVector vector = GVector(10); // Allocates a vector with 10 elements - -Note that there are 3 special constructors which you can use to allocate -vectors with 1, 2 or 3 elements: - -.. code-block:: cpp - - GVector vector(1.0); // Allocates vector [1.0] - GVector vector(1.0, 2.0); // Allocates vector [1.0, 2.0] - GVector vector(1.0, 2.0, 3.0); // Allocates vector [1.0, 2.0, 3.0] - -You can access vector elements using the :doxy:`GVector::operator[]` operator: - -.. code-block:: cpp - - GVector vector(10); // Allocates a vector with 10 elements - for (int i = 0; i < 10; ++i) { - vector[i] = (i+1)*10.0; // Set elements 10, 20, ..., 100 - } - for (int i = 0; i < 10; ++i) { - std::cout << vector[i] << std::endl; // Dump all elements, one by row - } - -Note that the :doxy:`GVector::operator[]` operator does not check the validity of the element -index. If index checking is needed, use the :doxy:`GVector::at` method: - - -.. code-block:: cpp - - GVector vector(10); // Allocates a vector with 10 elements - for (int i = 0; i < 10; ++i) { - vector.at(i) = (i+1)*10.0; // Set elements 10, 20, ..., 100 - } - for (int i = 0; i < 10; ++i) { - std::cout << vector.at(i) << std::endl; // Dump all elements, one by row - } - -The :doxy:`GVector::at` method will throw an :doxy:`GException::out_of_range` exception -in case that the specified index is not valid. - -You may also dump the content of a vector using - -.. code-block:: cpp - - std::cout << vector << std::endl; // Dump entire vector - -which in the above example will put the sequence :: - - (10, 20, 30, 40, 50, 60, 70, 80, 90, 100) - -on the screen. - - -Vector arithmetics -^^^^^^^^^^^^^^^^^^ - -You can handle vectors pretty much the same way you handle floating point -variables. :doxy:`GVector` supports various arithmetic operations: - -.. code-block:: cpp - - GVector a; // A vector - GVector b; // Another vector - GVector c; // Yet another vector - double s; // A double precision value - ... - c = a + b; // Vector + Vector addition - c = a + s; // Vector + Scalar addition - c = s + b; // Scalar + Vector addition - c = a - b; // Vector - Vector subtraction - c = a - s; // Vector - Scalar subtraction - c = s - b; // Scalar - Vector subtraction - s = a * b; // Vector * Vector multiplication (dot product) - c = a * s; // Vector * Scalar multiplication - c = s * b; // Scalar * Vector multiplication - c = a / s; // Vector * Scalar division - -Most of these operations operate element-wise. For example, scalar -additions or subtractions add or subtract a given scalar value from every -vector element. And scalar multiplications and divisions multiply or -divide every vector element by a given value. The dot product implements -the usual formula - -.. math:: - s = \sum_{i=0}^{N-1} a_i b_i - -(where :math:`N` is the number of vector elements). -It is obvious that the dot product, as well as vector addition and -subtraction, require vectors of identical dimensions. If vectors are not -identical, an :doxy:`GException::vector_mismatch` exception will be thrown: - -.. code-block:: cpp - - try { - GVector a(10); - GVector b(11); - GVector c = a + b; // WRONG: Vectors have incompatible dimensions - } - catch (GVector::vector_mismatch &e) { - std::cout << e.what() << std::endl; // Dimension exception is catched here - throw; - } - -Further vector operations are - -.. code-block:: cpp - - c = a; // Vector assignment - c = s; // Assigns scalar to all vector elements - s = c[index]; // Vector element access - c += a; // c = c + a; - c -= a; // c = c - a; - c += s; // c = c + s; - c -= s; // c = c - s; - c *= s; // c = c * s; - c /= s; // c = c / s; - c = -a; // Vector negation - -Finally, you can use the comparison operators - -.. code-block:: cpp - - int equal = (a == b); // True if all elements equal - int unequal = (a != b); // True if at least one elements unequal - -to compare all elements of a vector. If all elements are -identical, the :doxy:`GVector::operator==` operator returns ``true``, otherwise ``false``. -If at least one element differs, the :doxy:`GVector::operator!=` operator returns ``true``, -if all elements are identical it returns ``false``. - -In addition to the operators, you can apply the following mathematical -functions to vectors:: - - acos atan exp sin tanh - acosh atanh abs sinh pow - asin cos log sqrt - asinh cosh log10 tan - -Again, these functions should be understood to be applied element-wise. -They all take a vector as argument and produce a vector as result. For -example - -.. code-block:: cpp - - c = sin(a); - -attributes the sine of each element of vector ``a`` to vector ``c``. -Additional implemented functions are - -.. code-block:: cpp - - c = cross(a, b); // Vector cross product (for 3d only) - s = norm(a); // Vector norm |a| - s = min(a); // Minimum element of vector - s = max(a); // Maximum element of vector - s = sum(a); // Sum of vector elements - -Finally, the following methods exist: - -.. code-block:: cpp - - int n = a.size(); // Returns number of vector elements - int n = a.non_zeros(); // Returns number of non-zero vector elements - - -Matrixes -~~~~~~~~ - -General -^^^^^^^ - -A matrix is a two-dimensional array of double precision floating point -values, arranged in rows and columns. Three matrix storage classes are -implemented in GammaLib: - -- :doxy:`GMatrix` which explicitly stores all elements of a matrix - -- :doxy:`GMatrixSymmetric` which implements a symmetric matrix and only stores - the lower-left triangle of matrix elements - -- :doxy:`GMatrixSparse` which implements a sparse matrix and only stores the - non-zero elements of a matrix - -All matrix classes derive from the abstract :doxy:`GMatrixBase` class. - - -.. _sec_matrix_storage: - -Matrix storage classes -^^^^^^^^^^^^^^^^^^^^^^ - -In the most general case, the rows and columns of a matrix are stored in -a continuous array of :math:`{\tt rows} \times {\tt columns}` memory -locations. This storage type is referred to as a *full matrix*, and is -implemented by the class :doxy:`GMatrix`. Operations on full matrixes are in -general relatively fast, but memory requirements may be important to -hold all the elements. In general matrixes are stored column-wise -(or in column-major format). For example, the matrix :: - - 1 2 3 4 5 - 6 7 8 9 10 - 11 12 13 14 15 - -is stored in memory as :: - - | 1 6 11 | 2 7 12 | 3 8 13 | 4 9 14 | 5 10 15 | - -Many physical or mathematical problems treat with a subclass of matrixes -that is symmetric, i.e. for which the element :math:`(row,col)` is identical to -the element :math:`(col,row)`. In this case, the duplicated elements need not to -be stored. The class :doxy:`GMatrixSymmetric` implements such a storage type. -:doxy:`GMatrixSymmetric` stores the lower-left triangle of the matrix in -column-major format. For illustration, the matrix :: - - 1 2 3 4 - 2 5 6 7 - 3 6 8 9 - 4 7 9 10 - -is stored in memory as :: - - | 1 2 3 4 | 5 6 7 | 8 9 | 10 | - -This divides the storage requirements to hold the matrix elements by -almost a factor of two. - -Finally, quite often one has to deal with matrixes that contain a large -number of zeros. Such matrixes are called *sparse matrixes*. If only the -non-zero elements of a sparse matrix are stored the memory requirements -are considerably reduced. This goes however at the expense of matrix -element access, which has become now more complex. In particular, -filling efficiently a sparse matrix is a non-trivial problem (see -:ref:`sec_matrix_filling`). Sparse matrix storage is implemented by -the :doxy:`GMatrixSparse` class. A :doxy:`GMatrixSparse` object contains -three one-dimensional arrays to store the matrix elements: a double type -array that contains in continuous column-major order all non-zero -elements, an int type array that contains for each non-zero element the -row number of its location, and an int type array that contains the -storage location of the first non-zero element for each matrix column. -To illustrate this storage format, the matrix :: - - 1 0 0 7 - 2 5 0 0 - 3 0 6 0 - 4 0 0 8 - -is stored in memory as :: - - | 1 2 3 4 | 5 | 6 | 7 8 | Matrix elements - | 0 1 2 3 | 1 | 2 | 0 3 | Row indices for all elements - | 0 | 4 | 5 | 6 | Storage location of first element of each column - -This example is of course not very economic, since the total number of -Bytes used to store the matrix is -:math:`8 \times 8 + (8 + 4) \times 4 = 112` Bytes, while a full -:math:`4 \times 4` matrix is stored in -:math:`(4 \times 4) \times 8 = 128` Bytes (recall: a double type values -takes 8 Bytes, an int type value takes 4 Bytes). For realistic large -systems, however, the gain in memory space can be dramatical. - -The usage of the :doxy:`GMatrix`, :doxy:`GMatrixSymmetric` and :doxy:`GMatrixSparse` -classes is analoguous in that they implement basically all functions and -methods in an identical way. So from the semantics the user has not to worry -about the storage class. However, matrix element access speeds are not -identical for all storage types, and if performance is an issue (as it -certainly always will be), the user has to consider matrix access more -carefully (see :ref:`sec_matrix_filling`). - -You allocate a matrix using the constructors: - -.. code-block:: cpp - - GMatrix A(10,20); // Full 10 x 20 matrix - GMatrixSymmetric B(10,10); // Symmetric 10 x 10 matrix - GMatrixSparse C(1000,10000); // Sparse 1000 x 10000 matrix - - GMatrix A(0,0); // WRONG: empty matrix not allowed - GMatrixSymmetric B(20,22); // WRONG: symmetric matrix requested - -In the constructor, the first argument specifies the number of rows, the -second the number of columns: ``A(row,column)``. A symmetric matrix needs of -course an equal number of rows and columns. And an empty matrix is not -allowed. All matrix elements are initialised to 0 by the matrix -allocation. - -You can access matrix elements using the ``()`` operator, with the first -argument specifying the row and the second argument the column to be -accessed (row and column indices run from 0 to the number of elements -minus one): - -.. code-block:: cpp - - for (int row = 0; row < n_rows; ++row) { - for (int col = 0; col < n_cols; ++col) { - A(row,col) = (row+col)/2.0; // Set value of matrix element - } - } - ... - double sum2 = 0.0; - for (int row = 0; row < n_rows; ++row) { - for (int col = 0; col < n_cols; ++col) { - sum2 *= A(row,col) * A(row,col); // Get value of matrix element - } - } - -You can dump the content of a matrix to the console using - -.. code-block:: cpp - - std::cout << A << std::endl; // Dump matrix - - -Matrix arithmetics -^^^^^^^^^^^^^^^^^^ - -The following description of matrix arithmetics applies to all storage -classes (see :ref:`sec_matrix_storage`). The following matrix -operators have been implemented: - -.. code-block:: cpp - - GMatrix A; - GMatrix B; - GMatrix C; - ... - C = A + B; // Matrix Matrix addition - C = A - B; // Matrix Matrix subtraction - C = A * B; // Matrix Matrix multiplication - C = A * v; // Matrix Vector multiplication - C = A * s; // Matrix Scalar multiplication - C = s * A; // Scalar Matrix multiplication - C = A / s; // Matrix Scalar division - C = -A; // Negation - A += B; // Matrix inplace addition - A -= B; // Matrix inplace subtraction - A *= B; // Matrix inplace multiplications - A *= s; // Matrix inplace scalar multiplication - A /= s; // Matrix inplace scalar division - -The comparison operators - -.. code-block:: cpp - - int equal = (A == B); // True if all elements equal - int unequal = (A != B); // True if at least one elements unequal - -allow to compare all elements of a matrix. If all elements are -identical, the ``==`` operator returns ``true``, otherwise ``false``. -If at least one element differs, the ``!=`` operator returns true, -if all elements are identical it returns false. - - -General matrix methods and functions -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -A number of methods have been implemented to manipulate matrixes. The -methods described in this section are available for all storage classes. - -The methods - -.. code-block:: cpp - - int rows = A.rows(); // Returns number of rows in matrix - int cols = A.columns(); // Returns number of columns in matrix - int rows = A.size(); // Returns number of elements in matrix - -provide access to the matrix dimensions, the methods - -.. code-block:: cpp - - double sum = A.sum(); // Sum of all elements in matrix - double min = A.min(); // Returns minimum element of matrix - double max = A.max(); // Returns maximum element of matrix - double fill = A.fill(); // Returns fraction of non-zero elements - -inform about some matrix properties. The methods - -.. code-block:: cpp - - GVector row_vector = A.row(row); // Extract matrix row into vector - GVector column_vector = A.col(column); // Extract matrix column into vector - -extract entire rows and columns from a matrix into a vector. Conversely, -you may set entire rows or columns of a matrix using the methods - - A.row(row, row_vector); // Add vector to column - A.column(column, column_vector); // Puts vector in column - - -The methods - -.. code-block:: cpp - - A.insert_col(v_col,col); // Puts vector in column - A.add_col(v_col,col); // Add vector to column - -inserts or adds the elements of a vector into a matrix column. Note that -no row insertion routines have been implemented (so far) since they -would be less efficient (recall that all matrix types are stored in -column-major format). - -Conversion from one storage type to another is performed using - -.. code-block:: cpp - - B = A.convert_to_full(); // Converts A -> GMatrix - B = A.convert_to_sym(); // Converts A -> GMatrixSymmetric - B = A.convert_to_sparse(); // Converts A -> GMatrixSparse - -Note that ``convert_to_sym()`` can only be applied to a matrix that is -indeed symmetric. - -The transpose of a matrix can be obtained by using one of - -.. code-block:: cpp - - A.transpose(); // Transpose method - B = transpose(A); // Transpose function - -The absolute value of a matrix is provided by - -.. code-block:: cpp - - B = fabs(A); // B = |A| - - -Specific matrix methods and functions -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -Some methods do only exist for specific matrix storage classes. Although -these methods could in principle apply to all matrix classes, they have -for practical reasons not yet been implemented for all storage classes. -This may change in future versions of GammaLib. - - -GMatrix methods -~~~~~~~~~~~~~~~ - - -GMatrixSymmetric methods -~~~~~~~~~~~~~~~~~~~~~~~~ - - -GMatrixSparse methods -~~~~~~~~~~~~~~~~~~~~~ - -Other methods -~~~~~~~~~~~~~ - -You can extract the lower or upper triangle of a matrix into another -matrix using - -.. code-block:: cpp - - GMatrix B = A.extract_lower_triangle(); // B holds lower triangle - GMatrix B = A.extract_upper_triangle(); // B holds upper triangle - -This method is implemented for storage classes :doxy:`GMatrix` and -:doxy:`GMatrixSymmetric`. - - - -Matrix factorisations -^^^^^^^^^^^^^^^^^^^^^ - -A general tool of numeric matrix calculs is factorisation. - -Solve linear equation ``Ax = b``. Inverse a matrix (by solving successively -``Ax = e``, where ``e`` are the unit vectors for all dimensions). - -For symmetric and positive definite matrices the most efficient -factorisation is the Cholesky decomposition. The following code fragment -illustrates the usage: - -.. code-block:: cpp - - GMatrix A(n_rows, n_cols}; - GVector x(n_rows}; - GVector b(n_rows}; - ... - A.cholesky_decompose(); // Perform Cholesky factorisation - x = A.cholesky_solver(b); // Solve Ax=b for x - -Note that once the function ``A.cholesky_decompose()`` has been applied, -the original matrix content has been replaced by its Cholesky -decomposition. Since the Cholesky decomposition can be performed inplace -(i.e. without the allocation of additional memory to hold the result), -the matrix replacement is most memory economic. In case that the -original matrix should be kept, one may either copy it before into -another :doxy:`GMatrix` object or use the function - -.. code-block:: cpp - - GMatrix L = cholesky_decompose(A); - x = L.cholesky_solver(b); - -A symmetric and positive definite matrix can be inverted using the -Cholesky decomposition using - -.. code-block:: cpp - - A.cholesky_invert(); // Inverse matrix using Cholesky fact. - -Alternatively, the function - -.. code-block:: cpp - - GMatrix A_inv = cholesky_invert(A); - -may be used. - -The Cholesky decomposition, solver and inversion routines may also be -applied to matrices that contain rows or columns that are filled by -zeros. In this case the functions provide the option to (logically) -compress the matrices by skipping the zero rows and columns during the -calculation. - -For compressed matrix Cholesky factorisation, only the non-zero rows and -columns have to be symmetric and positive definite. In particular, the -full matrix may even be non-symmetric. - -Sparse matrixes -^^^^^^^^^^^^^^^ - -The only exception that does not work is - -.. code-block:: cpp - - GMatrixSparse A(10,10); - A(0,0) = A(1,1) = A(2,2) = 1.0; // WRONG: Cannot assign multiple at once - -In this case the value ``1.0`` is only assigned to the last element, i.e. -``A(2,2)``, the other elements will remain ``0``. This feature has to do with -the way how the compiler translates the code and how  implements sparse -matrix filling. :doxy:`GMatrixSparse` provides a pointer for a new element to be -filled. Since there is only one such *fill pointer*, only one element -can be filled at once in a statement. **So it is strongly advised to -avoid multiple matrix element assignment in a single row.** Better write -the above code like - -.. code-block:: cpp - - GMatrixSparse A; - A(0,0) = 1.0; - A(1,1) = 1.0; - A(2,2) = 1.0; - -This way, element assignment works fine. - -Inverting a sparse matrix produces in general a full matrix, so the -inversion function should be used with caution. Note that a full matrix -that is stored in sparse format takes roughly twice the memory than a -normal :doxy:`GMatrix` object. If nevertheless the inverse of a sparse matrix -should be examined, it is recommended to perform the analysis -column-wise - -.. code-block:: cpp - - GMatrixSparse A(rows,cols); // Allocate sparse matrix - GVector unit(rows); // Allocate vector - ... - A.cholesky_decompose(); // Factorise matrix - - // Column-wise solving the matrix equation - for (int col = 0; col < cols; ++col) { - unit(col) = 1.0; // Set unit vector - GVector x = cholesky_solver(unit); // Get column x of inverse - ... - unit(col) = 0.0; // Clear unit vector for next round - } - -.. _sec_matrix_filling: - -Filling sparse matrixes -^^^^^^^^^^^^^^^^^^^^^^^ - -The filling of a sparse matrix is a tricky issue since the storage of -the elements depends on their distribution in the matrix. If one would -know beforehand this distribution, sparse matrix filling would be easy -and fast. In general, however, the distribution is not known a priori, -and matrix filling may become a quite time consuming task. - -If a matrix has to be filled element by element, the access through the -operator - -.. code-block:: cpp - - m(row,col) = value; - -may be mandatory. In principle, if a new element is inserted into a -matrix a new memory cell has to be allocated for this element, and other -elements may be moved. Memory allocation is quite time consuming, and to -reduce the overhead, :doxy:`GMatrixSparse` can be configured to allocate memory -in bunches. By default, each time more matrix memory is needed, -:doxy:`GMatrixSparse` allocates 512 cells at once (or 6144 Bytes since each -element requires a double and a int storage location). If this amount of -memory is not adequat one may change this value by using - -.. code-block:: cpp - - m.set_mem_block(size); - -where size is the number of matrix elements that should be allocated at -once (corresponding to a total memory of :math:`12 \times {\tt size}` -Bytes). - -Alternatively, a matrix may be filled column-wise using the functions - -.. code-block:: cpp - - m.insert_col(vector,col); // Insert a vector in column - m.add_col(vector,col); // Add content of a vector to column - -While ``insert_col`` sets the values of column ``col`` (deleting thus any -previously existing entries), ``add_col`` adds the content of vector to all -elements of column col. Using these functions is considerably more rapid -than filling individual values. - -Still, if the matrix is big (i.e. several thousands of rows and -columns), filling individual columns may still be slow. To speed-up -dynamical matrix filling, an internal fill-stack has been implemented in -:doxy:`GMatrixSparse`. Instead of inserting values column-by-column, the columns -are stored in a stack and filled into the matrix once the stack is full. -This reduces the number of dynamic memory allocations to let the matrix -grow as it is built. By default, the internal stack is disabled. The -stack can be enabled and used as follows: - -.. code-block:: cpp - - m.stack_init(size, entries); // Initialise stack - ... - m.add_col(vector,col); // Add columns - ... - m.stack_destroy(); // Flush and destory stack - -The method ``stack_init`` initialises a stack with a number of size -elements and a maximum of entries columns. The larger the values size -and entries are chosen, the more efficient the stack works. The total -amount of memory of the stack can be estimated as -:math:`12 \times {\tt size} + 8 \times {\tt entries}` Bytes. If a rough -estimate of the total number of non-zero elements is available it is -recommended to set size to this value. As a rule of thumb, size should -be at least of the dimension of either the number of rows or the number -of columns of the matrix (take the maximum of both). entries is best set -to the number of columns of the matrix. If memory limits are an issue -smaller values may be set, but if the values are too small, the speed -increase may become negligible (or stack-filling may even become slower -than normal filling). - -Stack-filling only works with the method ``add_col``. Note also that -filling sub-sequently the same column leads to stack flushing. In the -code - -.. code-block:: cpp - - for (int col = 0; col < 100; ++col) { - column = 0.0; // Reset column - column(col) = col; // Set column - m.add_col(column,col); // Add column - } - -stack flushing occurs in each loop, and consequently, the stack-filling -approach will be not very efficient (it would probably be even slover -than normal filling). If successive operations are to be performed on -columns, it is better to perform them before adding. The code - -.. code-block:: cpp - - column = 0.0; // Reset column - for (int col = 0; col < 100; ++col) - column(col) = col; // Set column - m.add_col(column,col); // Add column - -would be far more efficient. - -A avoidable overhead occurs for the case that the column to be added is -sparse. The vector may contain many zeros, and :doxy:`GMatrixSparse` has to -filter them out. If the sparsity of the column is known, this overhead -can be avoided by directly passing a compressed array to ``add_col``: - -.. code-block:: cpp - - int number = 5; // 5 elements in array - double* values = new double[number]; // Allocate values - int* rows = new int[number]; // Allocate row index - ... - m.stack_init(size, entries); // Initialise stack - ... - for (int i = 0; i < number; ++i) { // Initialise array - values[i] = ... // ... set values - rows[i] = ... // ... set row indices - } - ... - m.add_col(values,rows,number,col); // Add array - ... - m.stack_destroy(); // Flush and destory stack - ... - delete [] values; // Free array - delete [] rows; - -The method ``add_col`` calls the method ``stack_push_column`` for stack -filling. ``add_col`` is more general than ``stack_push_column`` in that it -decides which of stack- or direct filling is more adequate. In -particular, ``stack_push_column`` may refuse pushing a column onto the -stack if there is not enough space. In that case, ``stack_push_column`` -returns a non-zero value that corresponds to the number of non-zero -elements in the vector that should be added. However, it is recommended -to not use ``stack_push_column`` and call instead ``add_col``. - -The method ``stack_destroy`` is used to flush and destroy the stack. After -this call the stack memory is liberated. If the stack should be flushed -without destroying it, the method ``stack_flush`` may be used: - -.. code-block:: cpp - - m.stack_init(size, entries); // Initialise stack - ... - m.add_col(vector,col); // Add columns - ... - m.stack_flush(); // Simply flush stack - -Once flushed, the stack can be filled anew. - -Note that stack flushing is not automatic! This means, if one trys to -use a matrix for calculs without flushing, the calculs may be wrong. -**If a stack is used for filling, always flush the stack before using -the matrix.** diff --git a/doc/source/users/tutorial/sky/convert_coordinates.rst b/doc/source/users/tutorial/sky/convert_coordinates.rst index 452dca5..a28b32d 100644 --- a/doc/source/users/tutorial/sky/convert_coordinates.rst +++ b/doc/source/users/tutorial/sky/convert_coordinates.rst @@ -5,27 +5,25 @@ To convert from celestial to Galactic coordinates you can set a sky direction in celestial coordinates and read it back in Galactic coordinates. +**C++** + +.. code-block:: cpp + :linenos: + + GSkyDir dir; + dir.radec_deg(83.63,22.01); + double l = dir.l_deg(); + double b = dir.b_deg(); + std::cout << l << ", " << b << std::endl; + **Python** .. code-block:: python :linenos: - import gammalib dir=gammalib.GSkyDir() dir.radec_deg(83.63,22.01) l=dir.l_deg() b=dir.b_deg() print(l,b) (184.55973405309402, -5.7891829467816827) - -**C++** - -.. code-block:: cpp - :linenos: - - #include "GammaLib.hpp" - GSkyDir dir; - dir.radec_deg(83.63,22.01); - double l = dir.l_deg(); - double b = dir.b_deg(); - std::cout << l << ", " << b << std::endl; diff --git a/doc/source/users/tutorial/sky/convert_projections.rst b/doc/source/users/tutorial/sky/convert_projections.rst index c3cf7b9..b8f6c77 100644 --- a/doc/source/users/tutorial/sky/convert_projections.rst +++ b/doc/source/users/tutorial/sky/convert_projections.rst @@ -6,25 +6,22 @@ into a map in cartesian projection. Map conversion is performed using the += operator that adds the bilinearly interpolated intensity values from one map to another. The example code applies to any kind of map projections. -**Python** - -.. code-block:: python - :linenos: - - import gammalib - healpix = gammalib.GSkyMap("healpix.fits") - map = gammalib.GSkyMap("CAR","GAL",0.0,0.0,0.5,0.5,100,100) - map += healpix - map.save("carmap.fits") - **C++** .. code-block:: cpp :linenos: - #include "GammaLib.hpp" GSkyMap healpix("healpix.fits"); GSkyMap map("CAR","GAL",0.0,0.0,0.5,0.5,100,100); map += healpix; map.save("carmap.fits"); +**Python** + +.. code-block:: python + :linenos: + + healpix = gammalib.GSkyMap('healpix.fits') + map = gammalib.GSkyMap('CAR','GAL',0.0,0.0,0.5,0.5,100,100) + map += healpix + map.save('carmap.fits') diff --git a/doc/source/users/user_manual/modules/linalg/index.rst b/doc/source/users/user_manual/modules/linalg/index.rst index 3232688..47b9efc 100644 --- a/doc/source/users/user_manual/modules/linalg/index.rst +++ b/doc/source/users/user_manual/modules/linalg/index.rst @@ -7,4 +7,6 @@ Linear algebra :maxdepth: 1 overview - matrix + vector + matrix_storage + matrix_arithmetics diff --git a/doc/source/users/user_manual/modules/linalg/matrix.rst b/doc/source/users/user_manual/modules/linalg/matrix.rst deleted file mode 100644 index ed0d9d0..0000000 --- a/doc/source/users/user_manual/modules/linalg/matrix.rst +++ /dev/null @@ -1,127 +0,0 @@ -Matrix Storage classes -====================== - -Three matrix storage classes are implemented: - -* :doxy:`GMatrix` implements a plain matrix storage class which stores - all elements of the matrix in memory; -* :doxy:`GMatrixSymmetric` implements a symmetric matrix storage class - which stores only a triangle of matrix elements, imposing thus - strict matrix symmetry; -* :doxy:`GMatrixSparse` implements a sparse matrix storage class which - stores only non-zero matrix elements in a column-wise organisation. - -Storage class conversion constructors exist for all three classes -to transform one storage class into another: - -**Python** - -.. code-block:: python - :linenos: - - plain = gammalib.GMatrix(10,10) - symmetric = gammalib.GMatrixSymmetric(plain) - sparse = gammalib.GMatrixSparse(symmetric) - -**C++** - -.. code-block:: cpp - :linenos: - - GMatrix plain(10,10); - GMatrixSymmetric symmetric(plain); - GMatrixSparse sparse(symmetric); - -In the above example, a plain matrix is instantiated in line 1, the -plain matrix is converted into a symmetric matrix in line 2, and the -symmetric matrix is converted into a sparse matrix in line 3. -Additional complementary storage class conversion constructors exist, -but conversion to a symmetric matrix is of course only possible if the -matrix is indeed symmetric. - -Matrix elements are accessed using the ``operator()``. Filling of -sparse matrix elements using this operator is possible, although this -can be time consuming due to internal memory management. In general, -each fill of a new non-zero element needs to shift all elements -that are located after that element in memory. - -To reduce the memory management overhead in the filling of a sparse -matrix, methods have been implemented that allow to fill a -matrix column wise: - -**Python** - -.. code-block:: cpp - :linenos: - - sparse = gammalib.GMatrixSparse(10,5) - column = gammalib.GVector(10) - column[0] = 1.0 - column[1] = 2.0 - column[5] = 8.0 - sparse.column(0, column) - sparse.add_to_column(0, column) - -**C++** - -.. code-block:: cpp - :linenos: - - GMatrixSparse sparse(10,5); - GVector column(10); - column[0] = 1.0; - column[1] = 2.0; - column[5] = 8.0; - sparse.column(0, column); - sparse.add_to_column(0, column); - -Line 1 allocates a sparse matrix with 10 rows and 5 columns, line 2 -instantiates a vector with 10 elements. In lines 3-5, 3 elements of -the vector are set to specific values, all other elements will default -to 0. In line 6, the elements of the vector are set as the elements -of the first matrix column (column 0). Line 7 differs from line 6 in -that the elements are now not set but added to the existing matrix -elements. - -To further reduce the memory management overhead for the column-wise -fill of a sparse matrix, a "fill-stack" has been implemented. The -"fill-stack" is a buffer that implements a queue for columns that are -to be set or added to the matrix. The columns will be stored in this -"fill-stack" in the order they are provided, and only once the -"fill-stack" is full, or upon request, the "fill-stack" will be flushed -into memory. The "fill-stack" is used as follows: - -**C++** - -.. code-block:: cpp - :linenos: - - sparse.stack_init(size, entries); - ... - sparse.column(0, column); - ... - sparse.stack_flush(); - ... - sparse.stack_destroy(); - -The ``stack_init(size, entries)`` method initialises the "fill-stack", -where ``size`` is the size of the allocated memory buffer and ``entries`` -is the maximum number of columns that will be held by the buffer. -If ``size`` is set to 0 (the default value), a default ``size`` value of -512 is used. If ``entries`` is set to 0 (the default value), the number of -matrix columns is taken as default ``entries`` value. Note that a too large -number of elements will produce some overhead due to "fill-stack" -management, hence ``entries`` should not exceeed a value of the order of -10-100. - -The ``stack_flush()`` method flushes the stack, which is mandatory -before any usage of the matrix. Note that the "fill-stack" **is not -inserted automatically** before any matrix operation, hence manual stack -flushing is needed to make all filled matrix elements available for usage. -The ``stack_destroy()`` method will flush the stack and free all stack -elements. This method should be called once no filling is required anymore. -If ``stack_destroy()`` is called immediately after filling, no call to -``stack_flush()`` is needed as the ``stack_destroy()`` method flushes the -stack before destroying it. The matrix stack is also destroyed by the -sparse matrix destructor, hence manual stack destruction is not -mandatory. diff --git a/doc/source/users/user_manual/modules/linalg/matrix_arithmetics.rst b/doc/source/users/user_manual/modules/linalg/matrix_arithmetics.rst new file mode 100644 index 0000000..3afa283 --- /dev/null +++ b/doc/source/users/user_manual/modules/linalg/matrix_arithmetics.rst @@ -0,0 +1,102 @@ +Matrix arithmetics +================== + +The following description of matrix arithmetics applies to all storage classes. +The following matrix operators have been implemented: + +**C++** + +.. code-block:: cpp + :linenos: + + GMatrix A(10,10); // A 10 x 10 matrix + GMatrix B(10,10); // Another 10 x 10 matrix + GMatrix C; // Result matrix + GVector v(10); // Vector with 10 elements + double s = 2.0; // Floating point value + C = A + B; // Matrix addition + C = A - B; // Matrix subtraction + C = A * B; // Matrix multiplication + C = A * v; // Vector multiplication + C = A * s; // Right-handed scalar multiplication + C = s * A; // Left-handed scalar Matrix multiplication (only C++) + C = A / s; // Scalar division + C = -A; // Negation + A += B; // Unary matrix addition + A -= B; // Unary matrix subtraction + A *= B; // Unary matrix multiplications + A *= s; // Unary matrix scalar multiplication + A /= s; // Unary matrix scalar division + +**Python** + +.. code-block:: python + :linenos: + + A = gammalib.GMatrix(10,10) # A 10 x 10 matrix + B = gammalib.GMatrix(10,10) # Another 10 x 10 matrix + v = gammalib.GVector(10) # Vector with 10 elements + s = 2.0 # Floating point value + C = A + B # Matrix addition + C = A - B # Matrix subtraction + C = A * B # Matrix multiplication + C = A * v # Vector multiplication + C = A * s # Scalar multiplication + C = A / s # Scalar division + C = -A # Negation + A += B # Unary matrix addition + A -= B # Unary matrix subtraction + A *= B # Unary matrix multiplications + A *= s # Unary matrix scalar multiplication + A /= s # Unary matrix scalar division + +You can use the comparison operators + +**C++** + +.. code-block:: cpp + :linenos: + + int equal = (A == B); // True if all elements equal + int unequal = (A != B); // True if at least one elements unequal + +**Python** + +.. code-block:: python + :linenos: + + equal = (A == B) # True if all elements equal + unequal = (A != B) # True if at least one elements unequal + +In addition to the operators, you can apply the following mathematical functions +to a matrix: + +**C++** + +.. code-block:: cpp + :linenos: + + C = A.abs(); // Matrix with absolute values of all matrix elements + C = A.transpose(); // Transpose matrix + C = A.invert(); // Invert matrix + v = A.solve(v); // Solve matrix equation x = M x v + s = A.min(); // Minimum element of matrix + s = A.max(); // Maximum element of matrix + s = A.sum(); // Sum of matrix elements + +**Python** + +.. code-block:: python + :linenos: + + C = A.abs() # Matrix with absolute values of all matrix elements + C = A.transpose() # Transpose matrix + C = A.invert() # Invert matrix + v = A.solve(v) # Solve matrix equation x = M x v + s = A.min() # Minimum element of matrix + s = A.max() # Maximum element of matrix + s = A.sum() # Sum of matrix elements + +.. warning:: + The ``invert()`` and ``solve()`` methods are so far only implemented for + sparse matrices. diff --git a/doc/source/users/user_manual/modules/linalg/matrix_storage.rst b/doc/source/users/user_manual/modules/linalg/matrix_storage.rst new file mode 100644 index 0000000..3f29157 --- /dev/null +++ b/doc/source/users/user_manual/modules/linalg/matrix_storage.rst @@ -0,0 +1,248 @@ +Matrix storage classes +====================== + +A matrix is a two-dimensional array of double precision floating point +values, arranged in rows and columns. Three matrix storage classes are +implemented in GammaLib: + +- :doxy:`GMatrix` which explicitly stores all elements of a matrix + +- :doxy:`GMatrixSymmetric` which implements a symmetric matrix and only stores + the lower-left triangle of matrix elements + +- :doxy:`GMatrixSparse` which implements a sparse matrix and only stores the + non-zero elements of a matrix + +All matrix classes derive from the abstract :doxy:`GMatrixBase` class. + +In the most general case, the rows and columns of a matrix are stored in +a continuous array of :math:`{\tt rows} \times {\tt columns}` memory +locations. This storage type is referred to as a *full matrix*, and is +implemented by the class :doxy:`GMatrix`. Operations on full matrixes are in +general relatively fast, but memory requirements may be important to +hold all the elements. In general matrixes are stored column-wise +(or in column-major format). For example, the matrix :: + + 1 2 3 4 5 + 6 7 8 9 10 + 11 12 13 14 15 + +is stored in memory as :: + + | 1 6 11 | 2 7 12 | 3 8 13 | 4 9 14 | 5 10 15 | + +Many physical or mathematical problems treat with a subclass of matrixes +that is symmetric, i.e. for which the element :math:`(row,col)` is identical to +the element :math:`(col,row)`. In this case, the duplicated elements need not to +be stored. The class :doxy:`GMatrixSymmetric` implements such a storage type. +:doxy:`GMatrixSymmetric` stores the lower-left triangle of the matrix in +column-major format. For illustration, the matrix :: + + 1 2 3 4 + 2 5 6 7 + 3 6 8 9 + 4 7 9 10 + +is stored in memory as :: + + | 1 2 3 4 | 5 6 7 | 8 9 | 10 | + +This divides the storage requirements to hold the matrix elements by +almost a factor of two. + +Finally, quite often one has to deal with matrixes that contain a large +number of zeros. Such matrixes are called *sparse matrixes*. If only the +non-zero elements of a sparse matrix are stored the memory requirements +are considerably reduced. This goes however at the expense of matrix +element access, which has become now more complex. In particular, +filling efficiently a sparse matrix is a non-trivial problem (see +:ref:`sec_matrix_filling`). Sparse matrix storage is implemented by +the :doxy:`GMatrixSparse` class. A :doxy:`GMatrixSparse` object contains +three one-dimensional arrays to store the matrix elements: a double type +array that contains in continuous column-major order all non-zero +elements, an int type array that contains for each non-zero element the +row number of its location, and an int type array that contains the +storage location of the first non-zero element for each matrix column. +To illustrate this storage format, the matrix :: + + 1 0 0 7 + 2 5 0 0 + 3 0 6 0 + 4 0 0 8 + +is stored in memory as :: + + | 1 2 3 4 | 5 | 6 | 7 8 | Matrix elements + | 0 1 2 3 | 1 | 2 | 0 3 | Row indices for all elements + | 0 | 4 | 5 | 6 | Storage location of first element of each column + +This example is of course not very economic, since the total number of +Bytes used to store the matrix is +:math:`8 \times 8 + (8 + 4) \times 4 = 112` Bytes, while a full +:math:`4 \times 4` matrix is stored in +:math:`(4 \times 4) \times 8 = 128` Bytes (recall: a double type values +takes 8 Bytes, an int type value takes 4 Bytes). For realistic large +systems, however, the gain in memory space can be dramatical. + +The usage of the :doxy:`GMatrix`, :doxy:`GMatrixSymmetric` and :doxy:`GMatrixSparse` +classes is analoguous in that they implement basically all functions and +methods in an identical way. So from the semantics the user has not to worry +about the storage class. However, matrix element access speeds are not +identical for all storage types, and if performance is an issue (as it +certainly always will be), the user has to consider matrix access more +carefully (see :ref:`sec_matrix_filling`). + +You allocate a matrix using the constructors + +**C++** + +.. code-block:: cpp + :linenos: + + GMatrix A(10,20); // Full 10 x 20 matrix + GMatrixSymmetric B(10,10); // Symmetric 10 x 10 matrix + GMatrixSparse C(1000,10000); // Sparse 1000 x 10000 matrix + + GMatrix D(0,0); // WRONG: empty matrix not allowed + GMatrixSymmetric E(20,22); // WRONG: symmetric matrix requested + +**Python** + +.. code-block:: python + :linenos: + + A = gammalib.GMatrix(10,20) # Full 10 x 20 matrix + B = gammalib.GMatrixSymmetric(10,10) # Symmetric 10 x 10 matrix + C = gammalib.GMatrixSparse(1000,10000) # Sparse 1000 x 10000 matrix + + D = gammalib.GMatrix(0,0) # WRONG: empty matrix not allowed + E = gammalib.GMatrixSymmetric(20,22) # WRONG: symmetric matrix requested + + +In the constructor, the first argument specifies the number of rows, the +second the number of columns: ``A(row,column)``. A symmetric matrix needs of +course an equal number of rows and columns. And an empty matrix is not +allowed. All matrix elements are initialised to 0 by the matrix +allocation. + +Storage class conversion constructors exist for all three classes +to transform one storage class into another: + +**C++** + +.. code-block:: cpp + :linenos: + + GMatrix full(10,10); + GMatrixSymmetric symmetric(plain); + GMatrixSparse sparse(symmetric); + +**Python** + +.. code-block:: python + :linenos: + + full = gammalib.GMatrix(10,10) + symmetric = gammalib.GMatrixSymmetric(full) + sparse = gammalib.GMatrixSparse(symmetric) + +Matrix elements are accessed using + +**C++** + +.. code-block:: cpp + :linenos: + + std::cout << B(5,7) << std::endl; // Matrix access in C++ through () operator + +**Python** + +.. code-block:: python + :linenos: + + print(B[5,7]) # Matrix access in Python through [] operator + +Filling of sparse matrix elements using this operator is possible, although +this can be time consuming due to internal memory management. In general, +each fill of a new non-zero element needs to shift all elements +that are located after that element in memory. + +To reduce the memory management overhead in the filling of a sparse matrix, +methods have been implemented that allow to fill a matrix column wise: + +**C++** + +.. code-block:: cpp + :linenos: + + GMatrixSparse sparse(10,5); + GVector column(10); + column[0] = 1.0; + column[1] = 2.0; + column[5] = 8.0; + sparse.column(0, column); + sparse.add_to_column(0, column); + +**Python** + +.. code-block:: cpp + :linenos: + + sparse = gammalib.GMatrixSparse(10,5) + column = gammalib.GVector(10) + column[0] = 1.0 + column[1] = 2.0 + column[5] = 8.0 + sparse.column(0, column) + sparse.add_to_column(0, column) + +Line 1 allocates a sparse matrix with 10 rows and 5 columns, line 2 +instantiates a vector with 10 elements. In lines 3-5, 3 elements of +the vector are set to specific values, all other elements will default +to 0. In line 6, the elements of the vector are set as the elements +of the first matrix column (column 0). Line 7 differs from line 6 in +that the elements are now not set but added to the existing matrix +elements. + +To further reduce the memory management overhead for the column-wise +fill of a sparse matrix, a "fill-stack" has been implemented. The +"fill-stack" is a buffer that implements a queue for columns that are +to be set or added to the matrix. The columns will be stored in this +"fill-stack" in the order they are provided, and only once the +"fill-stack" is full, or upon request, the "fill-stack" will be flushed +into memory. The "fill-stack" is used as follows: + +**C++** + +.. code-block:: cpp + :linenos: + + sparse.stack_init(size, entries); + ... + sparse.column(0, column); + ... + sparse.stack_flush(); + ... + sparse.stack_destroy(); + +The ``stack_init(size, entries)`` method initialises the "fill-stack", +where ``size`` is the size of the allocated memory buffer and ``entries`` +is the maximum number of columns that will be held by the buffer. +If ``size`` is set to 0 (the default value), a default ``size`` value of +512 is used. If ``entries`` is set to 0 (the default value), the number of +matrix columns is taken as default ``entries`` value. Note that a too large +number of elements will produce some overhead due to "fill-stack" +management, hence ``entries`` should not exceeed a value of the order of +10-100. + +The ``stack_flush()`` method flushes the stack, which is mandatory +before any usage of the matrix. Note that the "fill-stack" **is not +inserted automatically** before any matrix operation, hence manual stack +flushing is needed to make all filled matrix elements available for usage. +The ``stack_destroy()`` method will flush the stack and free all stack +elements. This method should be called once no filling is required anymore. +If ``stack_destroy()`` is called immediately after filling, no call to +``stack_flush()`` is needed as the ``stack_destroy()`` method flushes the +stack before destroying it. The matrix stack is also destroyed by the +sparse matrix destructor, hence manual stack destruction is not +mandatory. diff --git a/doc/source/users/user_manual/modules/linalg/vector.rst b/doc/source/users/user_manual/modules/linalg/vector.rst new file mode 100644 index 0000000..4b8646b --- /dev/null +++ b/doc/source/users/user_manual/modules/linalg/vector.rst @@ -0,0 +1,264 @@ +Vectors +======= + +The :doxy:`GVector` class implements a vector as a one-dimensional array of +successive double precision floating point values. You may allocate +an vectors as follows + +**C++** + +.. code-block:: cpp + :linenos: + + GVector vector1; // Allocates empty vector + GVector vector2(10); // Allocates vector with 10 elements + GVector vector3 = vector1; // Allocates vector by copying + +**Python** + +.. code-block:: python + :linenos: + + vector1 = gammalib.GVector() # Allocates empty vector + vector2 = gammalib.GVector(10) # Allocates vector with 10 elements + vector3 = vector1.copy() # Allocates vector by copying + +On allocation, all elements of a vector are set to 0. + +There are three special constructors which you can use to allocate vectors +with 1, 2 or 3 elements: + +**C++** + +.. code-block:: cpp + :linenos: + + GVector vector1(1.0); // Allocates vector [1.0] + GVector vector2(1.0, 2.0); // Allocates vector [1.0, 2.0] + GVector vector3(1.0, 2.0, 3.0); // Allocates vector [1.0, 2.0, 3.0] + +**Python** + +.. code-block:: python + :linenos: + + vector1 = gammalib.GVector(1.0) # Allocates vector [1.0] + vector2 = gammalib.GVector(1.0, 2.0) # Allocates vector [1.0, 2.0] + vector3 = gammalib.GVector(1.0, 2.0, 3.0) # Allocates vector [1.0, 2.0, 3.0] + +You can access vector elements using the :doxy:`GVector::operator[]` operator. +As usual in C++ and Python, indices start with 0: + +**C++** + +.. code-block:: cpp + :linenos: + + std::cout << vector2[1] << std::endl; // Print 2nd element + 2.0 // Result + +**Python** + +.. code-block:: python + :linenos: + + print(vector2[1]) # Print 2nd element + 2.0 # Result + +Note that in C++ the :doxy:`GVector::operator[]` operator does not check the +validity of the element index. If index checking is needed in C++, use the +:doxy:`GVector::at` method: + +**C++** + +.. code-block:: cpp + :linenos: + + std::cout << vector2.at(1) << std::endl; // Print 2nd element with index checking + 2.0 // Result + +In Python, index checking is always performed. + +All GammaLib objects are printable, and vectors can be printed using + +**C++** + +.. code-block:: cpp + :linenos: + + std::cout << vector3 << std::endl; // Print vector + (1, 2, 3) // Result + +**Python** + +.. code-block:: python + :linenos: + + print(vector3) # Print result + (1, 2, 3) # Result + + +You can handle vectors pretty much the same way you handle floating point +variables. :doxy:`GVector` supports various arithmetic operations: + +**C++** + +.. code-block:: cpp + :linenos: + + GVector a(1.0, 2.0, 3.0); // A vector + GVector b(2.0, 4.0, 6.0); // Another vector + GVector c; // Result vector + double s = 2.0; // A floating point value + c = a + b; // Vector + Vector addition + c += a; // Unary vector addition + c = a + s; // Vector + Scalar addition + c += s; // Unary scalar addition + c = s + b; // Scalar + Vector addition (only C++) + c = a - b; // Vector - Vector subtraction + c -= a; // Unary vector subtraction + c = a - s; // Vector - Scalar subtraction + c -= s; // Unary scalar subtraction + c = s - b; // Scalar - Vector subtraction (only C++) + s = a * b; // Vector * Vector multiplication (dot product) + c = a * s; // Vector * Scalar multiplication + c *= s; // Unary scalar multiplication + c = s * b; // Scalar * Vector multiplication (only C++) + c = a / s; // Vector * Scalar division + c /= s; // Unary scalar division + c = s; // Assigns scalar to all vector elements + c = -a; // Vector negation + +**Python** + +.. code-block:: python + :linenos: + + a = gammalib.GVector(1.0, 2.0, 3.0) # A vector + b = gammalib.GVector(2.0, 4.0, 6.0) # Another vector + c = gammalib.GVector() # Result vector + s = 2.0 # A floating point value + c = a + b # Vector + Vector addition + c += a # Unary vector addition + c = a + s # Vector + Scalar addition + c += s # Unary scalar addition + c = a - b # Vector - Vector subtraction + c -= a # Unary vector subtraction + c = a - s # Vector - Scalar subtraction + c -= s # Unary scalar subtraction + s = a * b # Vector * Vector multiplication (dot product) + c = a * s # Vector * Scalar multiplication + c *= s # Unary scalar multiplication + c = a / s # Vector * Scalar division + c /= s # Unary scalar division + c = -a # Vector negation + +Most of these operations operate element-wise. For example, scalar +additions or subtractions add or subtract a given scalar value from every +vector element. And scalar multiplications and divisions multiply or +divide every vector element by a given value. The dot product implements +the usual formula + +.. math:: + s = \sum_{i=0}^{N-1} a_i b_i + +(where :math:`N` is the number of vector elements). +It is obvious that the dot product, as well as vector addition and +subtraction, require vectors of identical dimensions. If vectors are not +identical, an :doxy:`GException::vector_mismatch` exception will be thrown. + +You can compare vectors using + +Finally, you can use the comparison operators + +**C++** + +.. code-block:: cpp + :linenos: + + int equal = (a == b); // True if all elements equal + int unequal = (a != b); // True if at least one elements unequal + +**Python** + +.. code-block:: python + :linenos: + + equal = (a == b) # True if all elements equal + unequal = (a != b) # True if at least one elements unequal + +In addition to the operators, you can apply the following mathematical functions +to vectors + +**C++** + +.. code-block:: cpp + :linenos: + + c = acos(a); // Vector with acos of all vector elements + c = acosh(a); // Vector with acosh of all vector elements + c = asin(a); // Vector with asin of all vector elements + c = asinh(a); // Vector with asinh of all vector elements + c = atan(a); // Vector with atan of all vector elements + c = atanh(a); // Vector with atanh of all vector elements + c = cos(a); // Vector with cos of all vector elements + c = cosh(a); // Vector with cosh of all vector elements + c = exp(a); // Vector with exp of all vector elements + c = abs(a); // Vector with absolute values of all vector elements + c = log(a); // Vector with natural logarithm of all vector elements + c = log10(a); // Vector with base 10 logarithm of all vector elements + c = sin(a); // Vector with sin of all vector elements + c = sinh(a); // Vector with sinh of all vector elements + c = tan(a); // Vector with tan of all vector elements + c = tanh(a); // Vector with tanh of all vector elements + c = pow(a, 2.7); // Vector with power of 2.7 of all vector elements (only C++) + c = cross(a, b); // Vector with cross product of 3 element vectors + s = norm(a); // Vector norm |a| + s = min(a); // Minimum element of vector + s = max(a); // Maximum element of vector + s = sum(a); // Sum of vector elements + +**Python** + +.. code-block:: python + :linenos: + + c = a.acos() # Vector with acos of all vector elements + c = a.acosh() # Vector with acosh of all vector elements + c = a.asin() # Vector with asin of all vector elements + c = a.asinh() # Vector with asinh of all vector elements + c = a.atan() # Vector with atan of all vector elements + c = a.atanh() # Vector with atanh of all vector elements + c = a.cos() # Vector with cos of all vector elements + c = a.cosh() # Vector with cosh of all vector elements + c = a.exp() # Vector with exp of all vector elements + c = a.abs() # Vector with absolute values of all vector elements + c = a.log() # Vector with natural logarithm of all vector elements + c = a.log10() # Vector with base 10 logarithm of all vector elements + c = a.sin() # Vector with sin of all vector elements + c = a.sinh() # Vector with sinh of all vector elements + c = a.tan() # Vector with tan of all vector elements + c = a.tanh() # Vector with tanh of all vector elements + c = a.cross(b) # Vector with cross product of 3 element vectors + s = a.norm() # Vector norm |a| + s = a.min() # Minimum element of vector + s = a.max() # Maximum element of vector + s = a.sum() # Sum of vector elements + +Finally, the following methods exist: + +**C++** + +.. code-block:: cpp + :linenos: + + int size = a.size(); // Returns number of vector elements + int nonzeros = a.non_zeros(); // Returns number of non-zero vector elements + +**Python** + +.. code-block:: python + :linenos: + + size = a.size(a) # Returns number of vector elements + nonzeros = a.non_zeros() # Returns number of non-zero vector elements -- libgit2 0.21.4