diff --git a/ChangeLog b/ChangeLog
index 10438c3..0200c21 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,8 +1,9 @@
-2016-10-07
+2016-10-08
* Version 1.2.0 released
=========================
+ Add GNdarray class (#1768)
Increase precision of gammalib::str() floating point functions
Add GFits::publish() methods (#1231)
Add GVOTable class (#1231)
diff --git a/NEWS b/NEWS
index efaf416..370994f 100644
--- a/NEWS
+++ b/NEWS
@@ -1,6 +1,6 @@
New Features and Important Changes in GammaLib 1.2.0
-7 October 2016
+8 October 2016
1. Introduction
@@ -16,6 +16,7 @@ The following classes have been added:
- GModelSpectralPlawEnergyFlux
- GModelSpectralMultiplicative
- GVOTable
+- GNdarray
The following classes have been removed:
- ...
@@ -155,7 +156,8 @@ product of spectral model components as a new spectral model (#1861).
9. Numerics module
------------------
-None
+The GNdarray class has been added that implements a n-dimensional double
+precision floating point array (#1768).
10. Observation module
diff --git a/include/GNdarray.hpp b/include/GNdarray.hpp
new file mode 100644
index 0000000..f67e322
--- /dev/null
+++ b/include/GNdarray.hpp
@@ -0,0 +1,532 @@
+/***************************************************************************
+ * GNdarray.hpp - N-dimensional array class *
+ * ----------------------------------------------------------------------- *
+ * copyright (C) 2016 by Juergen Knoedlseder *
+ * ----------------------------------------------------------------------- *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ * *
+ ***************************************************************************/
+/**
+ * @file GNdarray.hpp
+ * @brief N-dimensional array class interface definition
+ * @author Juergen Knoedlseder
+ */
+
+#ifndef GNDARRAY_HPP
+#define GNDARRAY_HPP
+
+/* __ Includes ___________________________________________________________ */
+#include
+#include
+#include "GBase.hpp"
+
+
+/***********************************************************************//**
+ * @class GNdarray
+ *
+ * @brief N-dimensional array class
+ *
+ * This class implement a n-dimensional double precision floating point
+ * array.
+ ***************************************************************************/
+class GNdarray : public GBase {
+
+ // Friend functions
+ friend double min(const GNdarray& array);
+ friend double max(const GNdarray& array);
+ friend double sum(const GNdarray& array);
+ friend GNdarray acos(const GNdarray& array);
+ friend GNdarray acosh(const GNdarray& array);
+ friend GNdarray asin(const GNdarray& array);
+ friend GNdarray asinh(const GNdarray& array);
+ friend GNdarray atan(const GNdarray& array);
+ friend GNdarray atanh(const GNdarray& array);
+ friend GNdarray cos(const GNdarray& array);
+ friend GNdarray cosh(const GNdarray& array);
+ friend GNdarray exp(const GNdarray& array);
+ friend GNdarray abs(const GNdarray& array);
+ friend GNdarray log(const GNdarray& array);
+ friend GNdarray log10(const GNdarray& array);
+ friend GNdarray sin(const GNdarray& array);
+ friend GNdarray sinh(const GNdarray& array);
+ friend GNdarray sqrt(const GNdarray& array);
+ friend GNdarray tan(const GNdarray& array);
+ friend GNdarray tanh(const GNdarray& array);
+ friend GNdarray pow(const GNdarray& array, const double& power);
+
+public:
+ // Constructors and destructors
+ GNdarray(void);
+ explicit GNdarray(const int& nx);
+ GNdarray(const int& nx, const int& ny);
+ GNdarray(const int& nx, const int& ny, const int& nz);
+ GNdarray(const std::vector& n);
+ GNdarray(const GNdarray& array);
+ virtual ~GNdarray(void);
+
+ // Element access operators
+ double& operator()(const int& ix);
+ double& operator()(const int& ix, const int& iy);
+ double& operator()(const int& ix, const int& iy, const int& iz);
+ double& operator()(const std::vector& i);
+ const double& operator()(const int& ix) const;
+ const double& operator()(const int& ix, const int& iy) const;
+ const double& operator()(const int& ix, const int& iy, const int& iz) const;
+ const double& operator()(const std::vector& i) const;
+
+ // Operators
+ GNdarray& operator=(const GNdarray& array);
+ bool operator==(const GNdarray& array) const;
+ bool operator!=(const GNdarray& array) const;
+ GNdarray& operator+=(const GNdarray& array);
+ GNdarray& operator-=(const GNdarray& array);
+ GNdarray& operator+=(const double& value);
+ GNdarray& operator-=(const double& value);
+ GNdarray& operator*=(const double& value);
+ GNdarray& operator/=(const double& value);
+ GNdarray operator-(void) const;
+
+ // Methods
+ void clear(void);
+ GNdarray* clone(void) const;
+ std::string classname(void) const;
+ int dim() const;
+ int size(void) const;
+ const std::vector& shape(void) const;
+ double& at(const int& ix);
+ double& at(const int& ix, const int& iy);
+ double& at(const int& ix, const int& iy, const int& iz);
+ double& at(const std::vector& i);
+ const double& at(const int& ix) const;
+ const double& at(const int& ix, const int& iy) const;
+ const double& at(const int& ix, const int& iy, const int& iz) const;
+ const double& at(const std::vector& i) const;
+ std::string print(const GChatter& chatter = NORMAL) const;
+
+#ifndef SWIG
+protected:
+#endif
+ // Protected methods
+ void init_members(void);
+ void copy_members(const GNdarray& array);
+ void free_members(void);
+ int index(const std::vector& i) const;
+ bool has_same_shape(const GNdarray& array) const;
+ void require_same_shape(const std::string& method, const GNdarray& array) const;
+
+ // Protected members
+ std::vector m_shape; //!< Array dimensions
+ std::vector m_strides; //!< Steps in each dimension when traversing array
+ std::vector m_data; //!< Array data
+};
+
+
+/***********************************************************************//**
+ * @brief Return class name
+ *
+ * @return String containing the class name ("GNdarray").
+ ***************************************************************************/
+inline
+std::string GNdarray::classname(void) const
+{
+ return ("GNdarray");
+}
+
+
+/***********************************************************************//**
+ * @brief 1-dimensional array element access operator
+ *
+ * @param[in] ix Element index [0,...,shape(0)-1].
+ * @return Reference to array element.
+ ***************************************************************************/
+inline
+double& GNdarray::operator()(const int& ix)
+{
+ // Return array element
+ return m_data[ix];
+}
+
+
+/***********************************************************************//**
+ * @brief 2-dimensional array element access operator
+ *
+ * @param[in] ix Index in first dimension [0,...,shape(0)-1].
+ * @param[in] iy Index in second dimension [0,...,shape(1)-1].
+ * @return Reference to array element.
+ ***************************************************************************/
+inline
+double& GNdarray::operator()(const int& ix, const int& iy)
+{
+ // Return array element
+ return m_data[ix+m_strides[1]*iy];
+}
+
+
+/***********************************************************************//**
+ * @brief 3-dimensional array element access operator
+ *
+ * @param[in] ix Index in first dimension [0,...,shape(0)-1].
+ * @param[in] iy Index in second dimension [0,...,shape(1)-1].
+ * @param[in] iz Index in third dimension [0,...,shape(2)-1].
+ * @return Reference to array element.
+ ***************************************************************************/
+inline
+double& GNdarray::operator()(const int& ix, const int& iy, const int& iz)
+{
+ // Return array element
+ return m_data[ix+m_strides[1]*iy+m_strides[2]*iz];
+}
+
+
+/***********************************************************************//**
+ * @brief n-dimensional array element access operator
+ *
+ * @param[in] i Index vector.
+ * @return Reference to array element.
+ ***************************************************************************/
+inline
+double& GNdarray::operator()(const std::vector& i)
+{
+ // Return array element
+ return m_data[index(i)];
+}
+
+
+/***********************************************************************//**
+ * @brief 1-dimensional array element access operator (const variant)
+ *
+ * @param[in] ix Element index [0,...,shape(0)-1].
+ * @return Const reference to array element.
+ ***************************************************************************/
+inline
+const double& GNdarray::operator()(const int& ix) const
+{
+ // Return array element
+ return m_data[ix];
+}
+
+
+/***********************************************************************//**
+ * @brief 2-dimensional array element access operator (const variant)
+ *
+ * @param[in] ix Index in first dimension [0,...,shape(0)-1].
+ * @param[in] iy Index in second dimension [0,...,shape(1)-1].
+ * @return Const reference to array element.
+ ***************************************************************************/
+inline
+const double& GNdarray::operator()(const int& ix, const int& iy) const
+{
+ // Return array element
+ return m_data[ix+m_strides[1]*iy];
+}
+
+
+/***********************************************************************//**
+ * @brief 3-dimensional array element access operator (const variant)
+ *
+ * @param[in] ix Index in first dimension [0,...,shape(0)-1].
+ * @param[in] iy Index in second dimension [0,...,shape(1)-1].
+ * @param[in] iz Index in third dimension [0,...,shape(2)-1].
+ * @return Const reference to array element.
+ ***************************************************************************/
+inline
+const double& GNdarray::operator()(const int& ix, const int& iy, const int& iz) const
+{
+ // Return array element
+ return m_data[ix+m_strides[1]*iy+m_strides[2]*iz];
+}
+
+
+/***********************************************************************//**
+ * @brief n-dimensional array element access operator (const variant)
+ *
+ * @param[in] i Index vector.
+ * @return Const reference to array element.
+ ***************************************************************************/
+inline
+const double& GNdarray::operator()(const std::vector& i) const
+{
+ // Return array element
+ return m_data[index(i)];
+}
+
+
+/***********************************************************************//**
+ * @brief Return dimension of array
+ *
+ * @return Dimension of array
+ *
+ * Returns the dimension of the array.
+ ***************************************************************************/
+inline
+int GNdarray::dim() const
+{
+ return m_shape.size();
+}
+
+
+/***********************************************************************//**
+ * @brief Return number of elements in array
+ *
+ * @return Number of elements in array
+ *
+ * Returns the number of elements in the array.
+ ***************************************************************************/
+inline
+int GNdarray::size() const
+{
+ return m_data.size();
+}
+
+
+/***********************************************************************//**
+ * @brief Return shape of array
+ *
+ * @return Shape of array
+ *
+ * Returns the shape of the array.
+ ***************************************************************************/
+inline
+const std::vector& GNdarray::shape(void) const
+{
+ return m_shape;
+}
+
+
+/***********************************************************************//**
+ * @brief 1-dimensional array element access with range checking
+ *
+ * @param[in] ix Element index [0,...,shape(0)-1].
+ * @return Reference to array element.
+ ***************************************************************************/
+inline
+double& GNdarray::at(const int& ix)
+{
+ return const_cast(static_cast(*this).at(ix));
+}
+
+
+/***********************************************************************//**
+ * @brief 2-dimensional array element access with range checking
+ *
+ * @param[in] ix Index in first dimension [0,...,shape(0)-1].
+ * @param[in] iy Index in second dimension [0,...,shape(1)-1].
+ * @return Reference to array element.
+ ***************************************************************************/
+inline
+double& GNdarray::at(const int& ix, const int& iy)
+{
+ return const_cast(static_cast(*this).at(ix,iy));
+}
+
+
+/***********************************************************************//**
+ * @brief 3-dimensional array element access with range checking
+ *
+ * @param[in] ix Index in first dimension [0,...,shape(0)-1].
+ * @param[in] iy Index in second dimension [0,...,shape(1)-1].
+ * @param[in] iz Index in third dimension [0,...,shape(2)-1].
+ * @return Reference to array element.
+ ***************************************************************************/
+inline
+double& GNdarray::at(const int& ix, const int& iy, const int& iz)
+{
+ return const_cast(static_cast(*this).at(ix,iy,iz));
+}
+
+
+/***********************************************************************//**
+ * @brief n-dimensional array element access with range checking
+ *
+ * @param[in] i Index vector.
+ * @return Reference to array element.
+ ***************************************************************************/
+inline
+double& GNdarray::at(const std::vector& i)
+{
+ return const_cast(static_cast(*this).at(i));
+}
+
+
+/***********************************************************************//**
+ * @brief Return sum of two arrays
+ *
+ * @param[in] a First array.
+ * @param[in] b Second array.
+ * @return Sum of arrays @p a and @p b.
+ *
+ * Returns the sum of arrays @p a and @p b.
+ ***************************************************************************/
+inline
+GNdarray operator+(const GNdarray& a, const GNdarray& b)
+{
+ GNdarray result = a;
+ result += b;
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Add value to array (right addition)
+ *
+ * @param[in] array Array.
+ * @param[in] value Value.
+ * @return Array with @p value added to all elements.
+ *
+ * Returns an array for which the @p value has been added to all elements.
+ ***************************************************************************/
+inline
+GNdarray operator+(const GNdarray& array, const double& value)
+{
+ GNdarray result = array;
+ result += value;
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Add value to array (left addition)
+ *
+ * @param[in] value Value.
+ * @param[in] array Array.
+ * @return Array with @p value added to all elements.
+ *
+ * Returns an array for which the @p value has been added to all elements.
+ ***************************************************************************/
+inline
+GNdarray operator+(const double& value, const GNdarray& array)
+{
+ GNdarray result = array;
+ result += value;
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Return difference of arrays
+ *
+ * @param[in] a First array.
+ * @param[in] b Second array.
+ * @return Difference between array @p a and @p b.
+ *
+ * Returns the difference between array @p a and @p b.
+ ***************************************************************************/
+inline
+GNdarray operator-(const GNdarray& a, const GNdarray& b)
+{
+ GNdarray result = a;
+ result -= b;
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Subtract value from array
+ *
+ * @param[in] array Array.
+ * @param[in] value Value.
+ * @return Array with @p value subtracted from all elements.
+ *
+ * Returns an array for which the @p value has been subtracted from all
+ * elements. For example
+ *
+ * double value = 5.0;
+ * GNdarray result = array - value;
+ ***************************************************************************/
+inline
+GNdarray operator-(const GNdarray& array, const double& value)
+{
+ GNdarray result = array;
+ result -= value;
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Subtract array from value
+ *
+ * @param[in] value Value.
+ * @param[in] array Array.
+ * @return Array with @p value subtracted from all elements.
+ *
+ * Returns an array for which all elements have been subtracted from
+ * @p value. For example
+ *
+ * double value = 5.0;
+ * GNdarray result = value - array;
+ ***************************************************************************/
+inline
+GNdarray operator-(const double& value, const GNdarray& array)
+{
+ GNdarray result = -array;
+ result += value;
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Multiply array by value (right multiplication)
+ *
+ * @param[in] array Array.
+ * @param[in] value Value.
+ * @return Array for which all elements have be multiplied by @p value.
+ *
+ * Returns an array for which all elements have be multiplied by @p value.
+ ***************************************************************************/
+inline
+GNdarray operator*(const GNdarray& array, const double& value)
+{
+ GNdarray result = array;
+ result *= value;
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Multiply array by value (left multiplication)
+ *
+ * @param[in] value Value.
+ * @param[in] array Array.
+ * @return Array for which all elements have be multiplied by @p value.
+ *
+ * Returns an array for which all elements have be multiplied by @p value.
+ ***************************************************************************/
+inline
+GNdarray operator*(const double& value, const GNdarray& array)
+{
+ GNdarray result = array;
+ result *= value;
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Divide array by value
+ *
+ * @param[in] array Array.
+ * @param[in] value Value.
+ * @return Array for which all elements have be divided by @p value.
+ *
+ * Returns an array for which all elements have be divided by @p value.
+ ***************************************************************************/
+inline
+GNdarray operator/(const GNdarray& array, const double& value)
+{
+ GNdarray result = array;
+ result /= value;
+ return result;
+}
+
+#endif /* GNDARRAY_HPP */
diff --git a/include/GammaLib.hpp b/include/GammaLib.hpp
index 159edc0..8bcd064 100644
--- a/include/GammaLib.hpp
+++ b/include/GammaLib.hpp
@@ -63,6 +63,7 @@
#include "GDerivative.hpp"
#include "GFunction.hpp"
#include "GMath.hpp"
+#include "GNdarray.hpp"
/* __ FITS module ________________________________________________________ */
#include "GFits.hpp"
diff --git a/include/Makefile.am b/include/Makefile.am
index 8f29e5b..de1cda3 100644
--- a/include/Makefile.am
+++ b/include/Makefile.am
@@ -44,6 +44,7 @@ pkginclude_HEADERS = GammaLib.hpp \
GDerivative.hpp \
GFunction.hpp \
GMath.hpp \
+ GNdarray.hpp \
GFits.hpp \
GFitsHDU.hpp \
GFitsHeader.hpp \
diff --git a/pyext/GNdarray.i b/pyext/GNdarray.i
new file mode 100755
index 0000000..8d70b7f
--- /dev/null
+++ b/pyext/GNdarray.i
@@ -0,0 +1,223 @@
+/***************************************************************************
+ * GNdarray.i - N-dimensional array class *
+ * ----------------------------------------------------------------------- *
+ * copyright (C) 2016 by Juergen Knoedlseder *
+ * ----------------------------------------------------------------------- *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ * *
+ ***************************************************************************/
+/**
+ * @file GNdarray.i
+ * @brief N-dimensional array class interface definition
+ * @author Juergen Knoedlseder
+ */
+
+%{
+/* Put headers and other declarations here that are needed for compilation */
+#define SWIG
+#include "GNdarray.hpp"
+%}
+
+
+/***********************************************************************//**
+ * @class GNdarray
+ *
+ * @brief N-dimensional array class
+ *
+ * This class implement a n-dimensional double precision floating point
+ * array.
+ ***************************************************************************/
+class GNdarray : public GBase {
+
+public:
+ // Constructors and destructors
+ GNdarray(void);
+ explicit GNdarray(const int& nx);
+ GNdarray(const int& nx, const int& ny);
+ GNdarray(const int& nx, const int& ny, const int& nz);
+ GNdarray(const std::vector& n);
+ GNdarray(const GNdarray& array);
+ virtual ~GNdarray(void);
+
+ // Operators
+ bool operator==(const GNdarray& array) const;
+ bool operator!=(const GNdarray& array) const;
+ GNdarray& operator+=(const GNdarray& array);
+ GNdarray& operator-=(const GNdarray& array);
+ GNdarray& operator+=(const double& value);
+ GNdarray& operator-=(const double& value);
+ GNdarray& operator*=(const double& value);
+ GNdarray operator-(void) const;
+
+ // Methods
+ void clear(void);
+ GNdarray* clone(void) const;
+ std::string classname(void) const;
+ int dim() const;
+ int size(void) const;
+ const std::vector& shape(void) const;
+};
+
+
+/***********************************************************************//**
+ * @brief GNdarray class extension
+ ***************************************************************************/
+%extend GNdarray {
+ double __getitem__(int GTuple[]) {
+ // Handle first the special case of a single index. In Python we want
+ // to support accessing a n-dimensional array with a single index so
+ // that we can iterate over the array
+ if (GTuple[0] == 1) {
+ if (GTuple[1] >= 0 && GTuple[1] < self->size()) {
+ return (self->m_data[GTuple[1]]);
+ }
+ else {
+ throw GException::out_of_range("__getitem__(int)", "Array "
+ "index", GTuple[1], self->size());
+ }
+ }
+ if (GTuple[0] != self->dim()) {
+ throw GException::invalid_value("__getitem__(int)", "Invalid "
+ "access of "+gammalib::str(self->dim())+"-dimensional "
+ "array with "+gammalib::str(GTuple[0])+"-dimensional "
+ "access operator.");
+ }
+ int index = 0;
+ for (int i = 0; i < self->shape().size(); ++i) {
+ if (GTuple[i+1] < 0 || GTuple[i+1] >= self->shape()[i]) {
+ throw GException::out_of_range("__getitem__(int)", "Dimension "+
+ gammalib::str(i)+" index", GTuple[i+1], self->shape()[i]);
+ }
+ index += self->m_strides[i] * GTuple[i+1];
+ }
+ return (self->m_data[index]);
+ }
+ void __setitem__(int GTuple[], const double& value) {
+ if (GTuple[0] != self->dim()) {
+ throw GException::invalid_value("__setitem__(int)", "Invalid "
+ "access of "+gammalib::str(self->dim())+"-dimensional "
+ "array with "+gammalib::str(GTuple[0])+"-dimensional "
+ "access operator.");
+ }
+ int index = 0;
+ for (int i = 0; i < self->shape().size(); ++i) {
+ if (GTuple[i+1] < 0 || GTuple[i+1] >= self->shape()[i]) {
+ throw GException::out_of_range("__setitem__(int)", "Dimension "+
+ gammalib::str(i)+" index", GTuple[i+1], self->shape()[i]);
+ }
+ index += self->m_strides[i] * GTuple[i+1];
+ }
+ self->m_data[index] = value;
+ }
+ GNdarray __add__(const GNdarray &a) {
+ return (*self) + a;
+ }
+ GNdarray __add__(const double &a) {
+ return (*self) + a;
+ }
+ GNdarray __sub__(const GNdarray &a) {
+ return (*self) - a;
+ }
+ GNdarray __sub__(const double &a) {
+ return (*self) - a;
+ }
+ GNdarray __mul__(const double &a) {
+ return (*self) * a;
+ }
+ // Python 2.x
+ GNdarray __div__(const double &a) {
+ return (*self) / a;
+ }
+ // Python 3.x
+ GNdarray __truediv__(const double &a) {
+ return (*self) / a;
+ }
+ // Python 2.x operator/=
+ GNdarray __idiv__(const double& value) {
+ self->operator/=(value);
+ return (*self);
+ }
+ // Python 3.x operator/=
+ GNdarray __itruediv__(const double& value) {
+ self->operator/=(value);
+ return (*self);
+ }
+ int __is__(const GNdarray &a) {
+ return (*self) == a;
+ }
+ GNdarray copy() {
+ return (*self);
+ }
+ double min() {
+ return min(*self);
+ }
+ double max() {
+ return max(*self);
+ }
+ double sum() {
+ return sum(*self);
+ }
+ GNdarray acos() {
+ return acos(*self);
+ }
+ GNdarray acosh() {
+ return acosh(*self);
+ }
+ GNdarray asin() {
+ return asin(*self);
+ }
+ GNdarray asinh() {
+ return asinh(*self);
+ }
+ GNdarray atan() {
+ return atan(*self);
+ }
+ GNdarray atanh() {
+ return atanh(*self);
+ }
+ GNdarray cos() {
+ return cos(*self);
+ }
+ GNdarray cosh() {
+ return cosh(*self);
+ }
+ GNdarray exp() {
+ return exp(*self);
+ }
+ GNdarray abs() {
+ return abs(*self);
+ }
+ GNdarray log() {
+ return log(*self);
+ }
+ GNdarray log10() {
+ return log10(*self);
+ }
+ GNdarray sin() {
+ return sin(*self);
+ }
+ GNdarray sinh() {
+ return sinh(*self);
+ }
+ GNdarray sqrt() {
+ return sqrt(*self);
+ }
+ GNdarray tan() {
+ return tan(*self);
+ }
+ GNdarray tanh() {
+ return tanh(*self);
+ }
+};
diff --git a/pyext/GTypemaps.i b/pyext/GTypemaps.i
index 6bfbbbb..cc1f6f2 100644
--- a/pyext/GTypemaps.i
+++ b/pyext/GTypemaps.i
@@ -78,17 +78,17 @@
/***********************************************************************//**
- * @brief Tuple to index conversion using variable dimensions (1,...,4)
+ * @brief Tuple to index conversion using variable dimensions
*
* The following function provides conversion between a Python tuple and
* an integer array. This allows index access via tuples, such as in
* a[3,5,10] = 10.0 or c = a[2,9].
***************************************************************************/
%{
-static int var_tuple_to_index(PyObject *input, int *ptr) {
+static int var_tuple_to_index(PyObject *input, int *ptr, int dim) {
if (PySequence_Check(input)) {
int size = PyObject_Length(input);
- if (size > 4) {
+ if (size > dim) {
PyErr_SetString(PyExc_ValueError,"Too many arguments in tuple");
return 0;
}
@@ -116,8 +116,8 @@ static int var_tuple_to_index(PyObject *input, int *ptr) {
}
}
%}
-%typemap(in) int GTuple[ANY] (int temp[5]) {
- if (!var_tuple_to_index($input,temp)) {
+%typemap(in) int GTuple[ANY] (int temp[11]) {
+ if (!var_tuple_to_index($input,temp,10)) {
return NULL;
}
$1 = &temp[0];
diff --git a/pyext/gammalib/numerics.i b/pyext/gammalib/numerics.i
index dc2ada8..1d492d3 100644
--- a/pyext/gammalib/numerics.i
+++ b/pyext/gammalib/numerics.i
@@ -38,8 +38,11 @@
/* __ Include standard typemaps for vectors and strings __________________ */
%include stl.i
+%include std_vector.i
+%template(vectori) std::vector;
/* __ Include GammaLib typemaps __________________________________________ */
+%include GTypemaps.i
%include typemap_GFilename.i
/* __ Include interface classes __________________________________________ */
@@ -53,3 +56,4 @@
%include "GFunction.i"
%include "GIntegral.i"
%include "GMath.i"
+%include "GNdarray.i"
diff --git a/src/numerics/GNdarray.cpp b/src/numerics/GNdarray.cpp
new file mode 100644
index 0000000..4310b09
--- /dev/null
+++ b/src/numerics/GNdarray.cpp
@@ -0,0 +1,1267 @@
+/***************************************************************************
+ * GNdarray.cpp - N-dimensional array class *
+ * ----------------------------------------------------------------------- *
+ * copyright (C) 2016 by Juergen Knoedlseder *
+ * ----------------------------------------------------------------------- *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ * *
+ ***************************************************************************/
+/**
+ * @file GNdarray.cpp
+ * @brief N-dimensional array class implementation
+ * @author Juergen Knoedlseder
+ */
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include
+#endif
+#include "GTools.hpp"
+#include "GException.hpp"
+#include "GNdarray.hpp"
+
+
+/* __ Method name definitions ____________________________________________ */
+#define G_OP_ADD "GNdarray::operator+=(GNdarray&)"
+#define G_OP_SUB "GNdarray::operator-=(GNdarray&)"
+#define G_AT1 "GNdarray::at(int&)"
+#define G_AT2 "GNdarray::at(int&, int&)"
+#define G_AT3 "GNdarray::at(int&, int&, int&)"
+#define G_ATN "GNdarray::at(std::vector&)"
+
+
+/*==========================================================================
+ = =
+ = Constructors/destructors =
+ = =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Void constructor
+ ***************************************************************************/
+GNdarray::GNdarray(void)
+{
+ // Initialise class members
+ init_members();
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief 1-dimensional constructor
+ *
+ * @param[in] nx Array dimension.
+ *
+ * Constructs a 1-dimensional array.
+ ***************************************************************************/
+GNdarray::GNdarray(const int& nx)
+{
+ // Initialise class members
+ init_members();
+
+ // Set shape
+ m_shape.push_back(nx);
+
+ // Set strides
+ m_strides.push_back(1);
+
+ // Initialise data array with zeros
+ m_data.assign(nx, 0.0);
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief 2-dimensional constructor
+ *
+ * @param[in] nx Array's first dimension.
+ * @param[in] ny Array's second dimension.
+ *
+ * Constructs a 2-dimensional array.
+ ***************************************************************************/
+GNdarray::GNdarray(const int& nx, const int& ny)
+{
+ // Initialise class members
+ init_members();
+
+ // Compute size of array
+ int size = nx * ny;
+
+ // Set shape
+ m_shape.push_back(nx);
+ m_shape.push_back(ny);
+
+ // Set strides
+ m_strides.push_back(1);
+ m_strides.push_back(nx);
+
+ // Initialise data array with zeros
+ m_data.assign(size, 0.0);
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief 3-dimensional constructor
+ *
+ * @param[in] nx Array's first dimension.
+ * @param[in] ny Array's second dimension.
+ * @param[in] nz Array's third dimension.
+ *
+ * Constructs a 3-dimensional array.
+ ***************************************************************************/
+GNdarray::GNdarray(const int& nx, const int& ny, const int& nz)
+{
+ // Initialise class members
+ init_members();
+
+ // Compute size of array
+ int size = nx * ny * nz;
+
+ // Set shape
+ m_shape.push_back(nx);
+ m_shape.push_back(ny);
+ m_shape.push_back(nz);
+
+ // Set strides
+ m_strides.push_back(1);
+ m_strides.push_back(nx);
+ m_strides.push_back(nx*ny);
+
+ // Initialise data array with zeros
+ m_data.assign(size, 0.0);
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief n-dimensional constructor
+ *
+ * @param[in] n Array dimensions.
+ *
+ * Constructs a n-dimensional array.
+ ***************************************************************************/
+GNdarray::GNdarray(const std::vector& n)
+{
+ // Initialise class members
+ init_members();
+
+ // Compute size of array
+ int size = n[0];
+ for (int i = 1; i < n.size(); ++i) {
+ size *= n[i];
+ }
+
+ // Set shape
+ m_shape = n;
+
+ // Set strides
+ int stride = 1;
+ for (int i = 0; i < n.size(); ++i) {
+ m_strides.push_back(stride);
+ stride *= n[i];
+ }
+
+ // Initialise data array with zeros
+ m_data.assign(size, 0.0);
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy constructor
+ *
+ * @param[in] array Array.
+ ***************************************************************************/
+GNdarray::GNdarray(const GNdarray& array)
+{
+ // Initialise class members
+ init_members();
+
+ // Copy members
+ copy_members(array);
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Destructor
+ ***************************************************************************/
+GNdarray::~GNdarray(void)
+{
+ // Free members
+ free_members();
+
+ // Return
+ return;
+}
+
+
+/*==========================================================================
+ = =
+ = Operators =
+ = =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Assignment operator
+ *
+ * @param[in] array Array.
+ * @return Array.
+ ***************************************************************************/
+GNdarray& GNdarray::operator=(const GNdarray& array)
+{
+ // Execute only if object is not identical
+ if (this != &array) {
+
+ // Free members
+ free_members();
+
+ // Initialise private members
+ init_members();
+
+ // Copy members
+ copy_members(array);
+
+ } // endif: object was not identical
+
+ // Return this object
+ return *this;
+}
+
+
+/***********************************************************************//**
+ * @brief Equality operator
+ *
+ * @param[in] array Array.
+ * @return True if arrays are identical.
+ *
+ * Returns true if the arrays have identical shape and content.
+ ***************************************************************************/
+bool GNdarray::operator==(const GNdarray& array) const
+{
+ // Check of arrays have the same shape
+ bool identity = has_same_shape(array);
+
+ // If the arrays have identical dimensions then check if all elements
+ // are identical. Break as soon as one array element differs.
+ if (identity) {
+ for (int i = 0; i < m_data.size(); ++i) {
+ if (m_data[i] != array.m_data[i]) {
+ identity = false;
+ break;
+ }
+ }
+ }
+
+ // Return result
+ return identity;
+}
+
+
+/***********************************************************************//**
+ * @brief Non-equality operator
+ *
+ * @param[in] array Array.
+ * @return True if arrays are different.
+ *
+ * Returns true if the arrays have different shape or content.
+ ***************************************************************************/
+bool GNdarray::operator!=(const GNdarray& array) const
+{
+ // Get negated result of equality operation
+ bool difference = !(this->operator==(array));
+
+ // Return result
+ return difference;
+}
+
+
+/***********************************************************************//**
+ * @brief Unary addition operator
+ *
+ * @param[in] array Array.
+ * @return Array.
+ *
+ * Adds an array to the current array.
+ ***************************************************************************/
+GNdarray& GNdarray::operator+=(const GNdarray& array)
+{
+ // Throw an exception if the arrays have not the same shape
+ require_same_shape(G_OP_ADD, array);
+
+ // Add elements
+ for (int i = 0; i < m_data.size(); ++i) {
+ m_data[i] += array.m_data[i];
+ }
+
+ // Return array
+ return *this;
+}
+
+
+/***********************************************************************//**
+ * @brief Unary subtraction operator
+ *
+ * @param[in] array Array.
+ * @return Array.
+ *
+ * Subtracts an array from the current array.
+ ***************************************************************************/
+GNdarray& GNdarray::operator-=(const GNdarray& array)
+{
+ // Throw an exception if the arrays have not the same shape
+ require_same_shape(G_OP_SUB, array);
+
+ // Subtract elements
+ for (int i = 0; i < m_data.size(); ++i) {
+ m_data[i] -= array.m_data[i];
+ }
+
+ // Return array
+ return *this;
+}
+
+
+/***********************************************************************//**
+ * @brief Unary value addition operator
+ *
+ * @param[in] value Value.
+ * @return Array.
+ *
+ * Adds a value to all array elements.
+ ***************************************************************************/
+GNdarray& GNdarray::operator+=(const double& value)
+{
+ // Add value to elements
+ for (int i = 0; i < m_data.size(); ++i) {
+ m_data[i] += value;
+ }
+
+ // Return array
+ return *this;
+}
+
+
+/***********************************************************************//**
+ * @brief Unary value subtraction operator
+ *
+ * @param[in] value Value.
+ * @return Array.
+ *
+ * Subtracts a value from all array elements.
+ ***************************************************************************/
+GNdarray& GNdarray::operator-=(const double& value)
+{
+ // Subtract value from elements
+ for (int i = 0; i < m_data.size(); ++i) {
+ m_data[i] -= value;
+ }
+
+ // Return array
+ return *this;
+}
+
+
+/***********************************************************************//**
+ * @brief Unary multiplication operator
+ *
+ * @param[in] value Value.
+ * @return Array.
+ *
+ * Multiply all array elements by a value.
+ ***************************************************************************/
+GNdarray& GNdarray::operator*=(const double& value)
+{
+ // Multiply all elements
+ for (int i = 0; i < m_data.size(); ++i) {
+ m_data[i] *= value;
+ }
+
+ // Return array
+ return *this;
+}
+
+
+/***********************************************************************//**
+ * @brief Unary division operator
+ *
+ * @param[in] value Value.
+ * @return Array.
+ *
+ * Divide all array elements by a value.
+ ***************************************************************************/
+GNdarray& GNdarray::operator/=(const double& value)
+{
+ // Divide all elements
+ for (int i = 0; i < m_data.size(); ++i) {
+ m_data[i] /= value;
+ }
+
+ // Return array
+ return *this;
+}
+
+
+/***********************************************************************//**
+ * @brief Unary minus operator
+ *
+ * @return Array.
+ *
+ * Negate all array elements.
+ ***************************************************************************/
+GNdarray GNdarray::operator-(void) const
+{
+ // Copy array
+ GNdarray result = *this;
+
+ // Negate all elements
+ for (int i = 0; i < m_data.size(); ++i) {
+ result.m_data[i] = -result.m_data[i];
+ }
+
+ // Return array
+ return result;
+}
+
+
+/*==========================================================================
+ = =
+ = Public methods =
+ = =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Clear array
+ ***************************************************************************/
+void GNdarray::clear(void)
+{
+ // Free members
+ free_members();
+
+ // Initialise private members
+ init_members();
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Clone array
+ *
+ * @return Pointer to deep copy of array.
+ ***************************************************************************/
+GNdarray* GNdarray::clone(void) const
+{
+ // Clone array
+ return new GNdarray(*this);
+}
+
+
+/***********************************************************************//**
+ * @brief 1-dimensional array element access with range checking
+ * (const version)
+ *
+ * @param[in] ix Element index [0,...,shape(0)-1].
+ * @return Const reference to array element.
+ *
+ * @exception GException::invalid_value
+ * Array is not 1-dimensional.
+ * @exception GException::out_of_range
+ * Element index is out of range.
+ ***************************************************************************/
+const double& GNdarray::at(const int& ix) const
+{
+ // Throw an exception if the array is not 1-dimensional
+ if (m_shape.size() != 1) {
+ std::string msg = "Invalid access of "+
+ gammalib::str(m_shape.size())+"-dimensional array "
+ "with 1-dimensional access operator.";
+ throw GException::invalid_value(G_AT1, msg);
+ }
+
+ // Throw an exception if the index is outside the valid range
+ if (ix < 0 || ix >= m_shape[0]) {
+ throw GException::out_of_range(G_AT1, "Array index", ix, m_shape[0]);
+ }
+
+ // Return array element
+ return (*this)(ix);
+}
+
+
+/***********************************************************************//**
+ * @brief 2-dimensional array element access with range checking
+ * (const version)
+ *
+ * @param[in] ix Index in first dimension [0,...,shape(0)-1].
+ * @param[in] iy Index in second dimension [0,...,shape(1)-1].
+ * @return Const reference to array element.
+ *
+ * @exception GException::invalid_value
+ * Array is not 2-dimensional.
+ * @exception GException::out_of_range
+ * Element index is out of range.
+ ***************************************************************************/
+const double& GNdarray::at(const int& ix, const int& iy) const
+{
+ // Throw an exception if the array is not 2-dimensional
+ if (m_shape.size() != 2) {
+ std::string msg = "Invalid access of "+
+ gammalib::str(m_shape.size())+"-dimensional array "
+ "with 2-dimensional access operator.";
+ throw GException::invalid_value(G_AT2, msg);
+ }
+
+ // Throw an exception if the indices are outside the valid range
+ if (ix < 0 || ix >= m_shape[0]) {
+ throw GException::out_of_range(G_AT2, "First array index", ix,
+ m_shape[0]);
+ }
+ if (iy < 0 || iy >= m_shape[1]) {
+ throw GException::out_of_range(G_AT2, "Second array index", iy,
+ m_shape[1]);
+ }
+
+ // Return array element
+ return (*this)(ix,iy);
+}
+
+
+/***********************************************************************//**
+ * @brief 3-dimensional array element access with range checking
+ * (const version)
+ *
+ * @param[in] ix Index in first dimension [0,...,shape(0)-1].
+ * @param[in] iy Index in second dimension [0,...,shape(1)-1].
+ * @param[in] iz Index in third dimension [0,...,shape(2)-1].
+ * @return Const reference to array element.
+ *
+ * @exception GException::invalid_value
+ * Array is not 3-dimensional.
+ * @exception GException::out_of_range
+ * Element index is out of range.
+ ***************************************************************************/
+const double& GNdarray::at(const int& ix, const int& iy, const int& iz) const
+{
+ // Throw an exception if the array is not 2-dimensional
+ if (m_shape.size() != 3) {
+ std::string msg = "Invalid access of "+
+ gammalib::str(m_shape.size())+"-dimensional array "
+ "with 3-dimensional access operator.";
+ throw GException::invalid_value(G_AT3, msg);
+ }
+
+ // Throw an exception if the indices are outside the valid range
+ if (ix < 0 || ix >= m_shape[0]) {
+ throw GException::out_of_range(G_AT3, "First array index", ix,
+ m_shape[0]);
+ }
+ if (iy < 0 || iy >= m_shape[1]) {
+ throw GException::out_of_range(G_AT3, "Second array index", iy,
+ m_shape[1]);
+ }
+ if (iz < 0 || iz >= m_shape[2]) {
+ throw GException::out_of_range(G_AT3, "Third array index", iz,
+ m_shape[2]);
+ }
+
+ // Return array element
+ return (*this)(ix,iy,iz);
+}
+
+
+/***********************************************************************//**
+ * @brief n-dimensional array element access with range checking
+ * (const version)
+ *
+ * @param[in] i Index vector.
+ * @return Const reference to array element.
+ *
+ * @exception GException::invalid_value
+ * Array is not of correct dimension.
+ * @exception GException::out_of_range
+ * Element index is out of range.
+ ***************************************************************************/
+const double& GNdarray::at(const std::vector& i) const
+{
+ // Throw an exception if the array is not 2-dimensional
+ if (m_shape.size() != i.size()) {
+ std::string msg = "Invalid access of "+
+ gammalib::str(m_shape.size())+"-dimensional array "
+ "with "+
+ gammalib::str(i.size())+"-dimensional access "
+ "operator.";
+ throw GException::invalid_value(G_ATN, msg);
+ }
+
+ // Throw an exception if the indices are outside the valid range
+ for (int k = 0; k < m_shape.size(); ++k) {
+ if (i[k] < 0 || i[k] >= m_shape[k]) {
+ throw GException::out_of_range(G_ATN, "Dimension "+
+ gammalib::str(k)+" index",
+ i[k], m_shape[k]);
+ }
+ }
+
+ // Return array element
+ return (*this)(i);
+}
+
+
+/***********************************************************************//**
+ * @brief Print array information
+ *
+ * @param[in] chatter Chattiness.
+ * @return String containing array information.
+ ***************************************************************************/
+std::string GNdarray::print(const GChatter& chatter) const
+{
+ // Initialise result string
+ std::string result = "(";
+
+ // Continue only if chatter is not silent
+ if (chatter != SILENT) {
+
+ // Put all elements in string
+ for (int i = 0; i < m_data.size(); ++i) {
+ if (i > 0) {
+ result += ", ";
+ }
+ result += gammalib::str(m_data[i]);
+ }
+
+ // Append closing parentheses
+ result += ")";
+
+ } // endif: chatter was not silent
+
+ // Return result
+ return result;
+}
+
+
+/*==========================================================================
+ = =
+ = Private methods =
+ = =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Initialise class members
+ ***************************************************************************/
+void GNdarray::init_members(void)
+{
+ // Initialise members
+ m_shape.clear();
+ m_strides.clear();
+ m_data.clear();
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy class members
+ *
+ * @param[in] array Array.
+ ***************************************************************************/
+void GNdarray::copy_members(const GNdarray& array)
+{
+ // Copy members
+ m_shape = array.m_shape;
+ m_strides = array.m_strides;
+ m_data = array.m_data;
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Delete class members
+ ***************************************************************************/
+void GNdarray::free_members(void)
+{
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Compute array element index
+ *
+ * @param[in] i Index vector.
+ * @return Array element index.
+ *
+ * Computes the array element index from an index vector. The method assumes
+ * that the index vector has the same dimension as the array. If this
+ * condition is not guaranteed, the condition needs to be checked before
+ * calling the method.
+ ***************************************************************************/
+int GNdarray::index(const std::vector& i) const
+{
+ // Initialise element index
+ int index = 0;
+
+ // Compute element index
+ for (int k = 0; k < m_shape.size(); ++k) {
+ index += m_strides[k] * i[k];
+ }
+
+ // Return index
+ return index;
+}
+
+
+/***********************************************************************//**
+ * @brief Check if array has the same shape
+ *
+ * @param[in] array Array.
+ * @return True if the array has the same shape.
+ ***************************************************************************/
+bool GNdarray::has_same_shape(const GNdarray& array) const
+{
+ // Check if the array has the same dimension
+ bool identity = m_shape.size() == array.m_shape.size();
+
+ // Check if the shape is identical. Break as soon as one shape value
+ // differs.
+ if (identity) {
+ for (int i = 0; i < m_shape.size(); ++i) {
+ if (m_shape[i] != array.m_shape[i]) {
+ identity = false;
+ break;
+ }
+ }
+ }
+
+ // Return result
+ return identity;
+}
+
+
+/***********************************************************************//**
+ * @brief Throw exception if array shapes differ
+ *
+ * @param[in] array Array.
+ ***************************************************************************/
+void GNdarray::require_same_shape(const std::string& method,
+ const GNdarray& array) const
+{
+ // If the shape differs then throw an exception
+ if (!has_same_shape(array)) {
+
+ // Compose exception message
+ std::string msg = "Incompatible array dimensions (";
+ for (int i = 0; i < m_shape.size(); ++i) {
+ if (i > 0) {
+ msg += ", ";
+ }
+ msg += gammalib::str(m_shape[i]);
+ }
+ msg += ") <=> (";
+ for (int i = 0; i < array.m_shape.size(); ++i) {
+ if (i > 0) {
+ msg += ", ";
+ }
+ msg += gammalib::str(array.m_shape[i]);
+ }
+ msg += ").";
+
+ // Throw exception
+ throw GException::invalid_argument(method, msg);
+
+ } // endif: array shapes differed
+
+ // Return
+ return;
+}
+
+
+/*==========================================================================
+ = =
+ = Friends =
+ = =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Computes minimum array element
+ *
+ * @param[in] array Array.
+ * @return Minimum array element.
+ ***************************************************************************/
+double min(const GNdarray& array)
+{
+ // Initialises result
+ double result = 0.0;
+
+ // Continue only if we have elements
+ if (array.m_data.size() > 0) {
+
+ // Search for minimum
+ result = array.m_data[0];
+ for (int i = 1; i < array.m_data.size(); ++i) {
+ if (array.m_data[i] < result) {
+ result = array.m_data[i];
+ }
+ }
+
+ } // endif: there were elements
+
+ // Returns minimum
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes maximum array element
+ *
+ * @param[in] array Array.
+ * @return Maximum array element.
+ ***************************************************************************/
+double max(const GNdarray& array)
+{
+ // Initialises result
+ double result = 0.0;
+
+ // Continue only if we have elements
+ if (array.m_data.size() > 0) {
+
+ // Search for maximum
+ result = array.m_data[0];
+ for (int i = 1; i < array.m_data.size(); ++i) {
+ if (array.m_data[i] > result) {
+ result = array.m_data[i];
+ }
+ }
+
+ } // endif: there were elements
+
+ // Returns maximum
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes array sum
+ *
+ * @param[in] array Array.
+ * @return Sum of array elements.
+ ***************************************************************************/
+double sum(const GNdarray& array)
+{
+ // Compute sum
+ double result = 0.0;
+ for (int i = 0; i < array.m_data.size(); ++i) {
+ result += array.m_data[i];
+ }
+
+ // Returns sum
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes arccos of array elements
+ *
+ * @param[in] array Array.
+ * @return Array containing the arccos of every element.
+ ***************************************************************************/
+GNdarray acos(const GNdarray& array)
+{
+ // Initialise result array
+ GNdarray result = array;
+
+ // Evaluate each array element
+ for (int i = 0; i < array.m_data.size(); ++i) {
+ result.m_data[i] = std::acos(array.m_data[i]);
+ }
+
+ // Return array
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes acosh of array elements
+ *
+ * @param[in] array Array.
+ * @return Array containing the acosh of every element.
+ ***************************************************************************/
+GNdarray acosh(const GNdarray& array)
+{
+ // Initialise result array
+ GNdarray result = array;
+
+ // Evaluate each array element
+ for (int i = 0; i < array.m_data.size(); ++i) {
+ result.m_data[i] = acosh(array.m_data[i]);
+ }
+
+ // Return array
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes arcsin of array elements
+ *
+ * @param[in] array Array.
+ * @return Array containing the arcsin of every element.
+ ***************************************************************************/
+GNdarray asin(const GNdarray& array)
+{
+ // Initialise result array
+ GNdarray result = array;
+
+ // Evaluate each array element
+ for (int i = 0; i < array.m_data.size(); ++i) {
+ result.m_data[i] = std::asin(array.m_data[i]);
+ }
+
+ // Return array
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes asinh of array elements
+ *
+ * @param[in] array Array.
+ * @return Array containing the asinh of every element.
+ ***************************************************************************/
+GNdarray asinh(const GNdarray& array)
+{
+ // Initialise result array
+ GNdarray result = array;
+
+ // Evaluate each array element
+ for (int i = 0; i < array.m_data.size(); ++i) {
+ result.m_data[i] = asinh(array.m_data[i]);
+ }
+
+ // Return array
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes arctan of array elements
+ *
+ * @param[in] array Array.
+ * @return Array containing the arctan of every element.
+ ***************************************************************************/
+GNdarray atan(const GNdarray& array)
+{
+ // Initialise result array
+ GNdarray result = array;
+
+ // Evaluate each array element
+ for (int i = 0; i < array.m_data.size(); ++i) {
+ result.m_data[i] = std::atan(array.m_data[i]);
+ }
+
+ // Return array
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes atanh of array elements
+ *
+ * @param[in] array Array.
+ * @return Array containing the atanh of every element.
+ ***************************************************************************/
+GNdarray atanh(const GNdarray& array)
+{
+ // Initialise result array
+ GNdarray result = array;
+
+ // Evaluate each array element
+ for (int i = 0; i < array.m_data.size(); ++i) {
+ result.m_data[i] = atanh(array.m_data[i]);
+ }
+
+ // Return array
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes cosine of array elements
+ *
+ * @param[in] array Array.
+ * @return Array containing the cosine of every element.
+ ***************************************************************************/
+GNdarray cos(const GNdarray& array)
+{
+ // Initialise result array
+ GNdarray result = array;
+
+ // Evaluate each array element
+ for (int i = 0; i < array.m_data.size(); ++i) {
+ result.m_data[i] = std::cos(array.m_data[i]);
+ }
+
+ // Return array
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes cosh of array elements
+ *
+ * @param[in] array Array.
+ * @return Array containing the cosh of every element.
+ ***************************************************************************/
+GNdarray cosh(const GNdarray& array)
+{
+ // Initialise result array
+ GNdarray result = array;
+
+ // Evaluate each array element
+ for (int i = 0; i < array.m_data.size(); ++i) {
+ result.m_data[i] = std::cosh(array.m_data[i]);
+ }
+
+ // Return array
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes exponential of array elements
+ *
+ * @param[in] array Array.
+ * @return Array containing the exponential of every element.
+ ***************************************************************************/
+GNdarray exp(const GNdarray& array)
+{
+ // Initialise result array
+ GNdarray result = array;
+
+ // Evaluate each array element
+ for (int i = 0; i < array.m_data.size(); ++i) {
+ result.m_data[i] = std::exp(array.m_data[i]);
+ }
+
+ // Return array
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes absolute of array elements
+ *
+ * @param[in] array Array.
+ * @return Array containing the absolute of every element.
+ ***************************************************************************/
+GNdarray abs(const GNdarray& array)
+{
+ // Initialise result array
+ GNdarray result = array;
+
+ // Evaluate each array element
+ for (int i = 0; i < array.m_data.size(); ++i) {
+ result.m_data[i] = std::abs(array.m_data[i]);
+ }
+
+ // Return array
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes natural logarithm of array elements
+ *
+ * @param[in] array Array.
+ * @return Array containing the natural logarithm of every element.
+ ***************************************************************************/
+GNdarray log(const GNdarray& array)
+{
+ // Initialise result array
+ GNdarray result = array;
+
+ // Evaluate each array element
+ for (int i = 0; i < array.m_data.size(); ++i) {
+ result.m_data[i] = std::log(array.m_data[i]);
+ }
+
+ // Return array
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes base10 logarithm of array elements
+ *
+ * @param[in] array Array.
+ * @return Array containing the base10 logarithm of every element.
+ ***************************************************************************/
+GNdarray log10(const GNdarray& array)
+{
+ // Initialise result array
+ GNdarray result = array;
+
+ // Evaluate each array element
+ for (int i = 0; i < array.m_data.size(); ++i) {
+ result.m_data[i] = std::log10(array.m_data[i]);
+ }
+
+ // Return array
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes sine of array elements
+ *
+ * @param[in] array Array.
+ * @return Array containing the sine of every element.
+ ***************************************************************************/
+GNdarray sin(const GNdarray& array)
+{
+ // Initialise result array
+ GNdarray result = array;
+
+ // Evaluate each array element
+ for (int i = 0; i < array.m_data.size(); ++i) {
+ result.m_data[i] = std::sin(array.m_data[i]);
+ }
+
+ // Return array
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes sinh of array elements
+ *
+ * @param[in] array Array.
+ * @return Array containing the sinh of every element.
+ ***************************************************************************/
+GNdarray sinh(const GNdarray& array)
+{
+ // Initialise result array
+ GNdarray result = array;
+
+ // Evaluate each array element
+ for (int i = 0; i < array.m_data.size(); ++i) {
+ result.m_data[i] = std::sinh(array.m_data[i]);
+ }
+
+ // Return array
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes square root of array elements
+ *
+ * @param[in] array Array.
+ * @return Array containing the square root of every element.
+ ***************************************************************************/
+GNdarray sqrt(const GNdarray& array)
+{
+ // Initialise result array
+ GNdarray result = array;
+
+ // Evaluate each array element
+ for (int i = 0; i < array.m_data.size(); ++i) {
+ result.m_data[i] = std::sqrt(array.m_data[i]);
+ }
+
+ // Return array
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes tangens of array elements
+ *
+ * @param[in] array Array.
+ * @return Array containing the tangens of every element.
+ ***************************************************************************/
+GNdarray tan(const GNdarray& array)
+{
+ // Initialise result array
+ GNdarray result = array;
+
+ // Evaluate each array element
+ for (int i = 0; i < array.m_data.size(); ++i) {
+ result.m_data[i] = std::tan(array.m_data[i]);
+ }
+
+ // Return array
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes tanh of array elements
+ *
+ * @param[in] array Array.
+ * @return Array containing the tanh of every element.
+ ***************************************************************************/
+GNdarray tanh(const GNdarray& array)
+{
+ // Initialise result array
+ GNdarray result = array;
+
+ // Evaluate each array element
+ for (int i = 0; i < array.m_data.size(); ++i) {
+ result.m_data[i] = std::tanh(array.m_data[i]);
+ }
+
+ // Return array
+ return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes tanh of array elements
+ *
+ * @param[in] array Array.
+ * @param[in] power Power.
+ * @return Array containing the power of every element.
+ ***************************************************************************/
+GNdarray pow(const GNdarray& array, const double& power)
+{
+ // Initialise result array
+ GNdarray result = array;
+
+ // Evaluate each array element
+ for (int i = 0; i < array.m_data.size(); ++i) {
+ result.m_data[i] = std::pow(array.m_data[i], power);
+ }
+
+ // Return array
+ return result;
+}
diff --git a/src/numerics/Makefile.am b/src/numerics/Makefile.am
index 4ee18f2..63ef3f4 100644
--- a/src/numerics/Makefile.am
+++ b/src/numerics/Makefile.am
@@ -29,6 +29,7 @@ sources = GIntegral.cpp \
GDerivative.cpp \
GFunction.cpp \
GMath.cpp \
+ GNdarray.cpp \
GException_numerics.cpp
# Build libtool library
diff --git a/test/test_GNumerics.cpp b/test/test_GNumerics.cpp
index e955c60..2b1c415 100644
--- a/test/test_GNumerics.cpp
+++ b/test/test_GNumerics.cpp
@@ -1,7 +1,7 @@
/***************************************************************************
* test_GNumerics.cpp - test numerics modules *
* ----------------------------------------------------------------------- *
- * copyright (C) 2010-2014 by Jurgen Knodlseder *
+ * copyright (C) 2010-2016 by Jurgen Knodlseder *
* ----------------------------------------------------------------------- *
* *
* This program is free software: you can redistribute it and/or modify *
@@ -50,10 +50,16 @@ void TestGNumerics::set(void)
m_sigma = 2.5;
// Append tests
- append(static_cast(&TestGNumerics::test_integral),"Test GIntegral");
- append(static_cast(&TestGNumerics::test_romberg_integration),"Test Romberg integration");
- append(static_cast(&TestGNumerics::test_adaptive_simpson_integration),"Test adaptive Simpson integration");
- append(static_cast(&TestGNumerics::test_gauss_kronrod_integration),"Test Gauss-Kronrod integration");
+ append(static_cast(&TestGNumerics::test_ndarray),
+ "Test GNdarray");
+ append(static_cast(&TestGNumerics::test_integral),
+ "Test GIntegral");
+ append(static_cast(&TestGNumerics::test_romberg_integration),
+ "Test Romberg integration");
+ append(static_cast(&TestGNumerics::test_adaptive_simpson_integration),
+ "Test adaptive Simpson integration");
+ append(static_cast(&TestGNumerics::test_gauss_kronrod_integration),
+ "Test Gauss-Kronrod integration");
// Return
return;
@@ -73,6 +79,203 @@ TestGNumerics* TestGNumerics::clone(void) const
/***********************************************************************//**
+ * @brief Test n-dimensional array
+ ***************************************************************************/
+void TestGNumerics::test_ndarray(void)
+{
+ // Test void constructor
+ GNdarray nd0;
+ test_value(nd0.size(), 0);
+ test_value(nd0.dim(), 0);
+
+ // Test 1-dimensional array
+ GNdarray nd1(3);
+ test_value(nd1.size(), 3);
+ test_value(nd1.dim(), 1);
+ test_value(nd1.shape()[0], 3);
+ for (int i = 0; i < 3; ++i) {
+ nd1(i) = double(i);
+ test_value(nd1.at(i), double(i));
+ }
+
+ // Test 2-dimensional array
+ GNdarray nd2(3,2);
+ test_value(nd2.size(), 6);
+ test_value(nd2.dim(), 2);
+ test_value(nd2.shape()[0], 3);
+ test_value(nd2.shape()[1], 2);
+ for (int ix = 0; ix < 3; ++ix) {
+ for (int iy = 0; iy < 2; ++iy) {
+ double ref = double(ix)+double(iy);
+ nd2(ix,iy) = ref;
+ test_value(nd2.at(ix,iy), ref);
+ }
+ }
+ // Test 3-dimensional array
+ GNdarray nd3(2,2,2);
+ test_value(nd3.size(), 8);
+ test_value(nd3.dim(), 3);
+ test_value(nd3.shape()[0], 2);
+ test_value(nd3.shape()[1], 2);
+ test_value(nd3.shape()[2], 2);
+ for (int ix = 0; ix < 2; ++ix) {
+ for (int iy = 0; iy < 2; ++iy) {
+ for (int iz = 0; iz < 2; ++iz) {
+ double ref = double(ix)+double(iy)+double(iz);
+ nd3(ix,iy,iz) = ref;
+ test_value(nd3.at(ix,iy,iz), ref);
+ }
+ }
+ }
+
+ // Test n-dimensional array
+ std::vector n(5,2);
+ GNdarray ndn(n);
+ test_value(ndn.size(), 32);
+ test_value(ndn.dim(), 5);
+ for (int i = 0; i < 5; ++i) {
+ test_value(ndn.shape()[i], 2);
+ }
+ std::vector itest(5,1);
+ ndn(itest) = 3.1415;
+ test_value(ndn.at(itest), 3.1415);
+
+ // Setup arrays for operator tests
+ GNdarray nda(1,2);
+ nda(0,0) = 3.0;
+ nda(0,1) = 5.5;
+ GNdarray ndb = nda;
+ GNdarray ndc = nda + ndb;
+
+ // Test operator=
+ test_value(ndb(0,0), 3.0);
+ test_value(ndb(0,1), 5.5);
+
+ // Test operator== and operator!=
+ test_assert(nda == ndb, "Test operator==");
+ test_assert(nda != ndc, "Test operator!=");
+
+ // Test array operator+=
+ nda += ndb;
+ test_value(nda(0,0), 6.0);
+ test_value(nda(0,1), 11.0);
+
+ // Test array operator-=
+ nda -= ndb;
+ test_value(nda(0,0), 3.0);
+ test_value(nda(0,1), 5.5);
+
+ // Test value operator+=
+ nda += 1.5;
+ test_value(nda(0,0), 4.5);
+ test_value(nda(0,1), 7.0);
+
+ // Test value operator-=
+ nda -= 1.5;
+ test_value(nda(0,0), 3.0);
+ test_value(nda(0,1), 5.5);
+
+ // Test value operator*=
+ nda *= 3.0;
+ test_value(nda(0,0), 9.0);
+ test_value(nda(0,1), 16.5);
+
+ // Test value operator/=
+ nda /= 3.0;
+ test_value(nda(0,0), 3.0);
+ test_value(nda(0,1), 5.5);
+
+ // Test operator-
+ GNdarray ndd = -nda;
+ test_value(ndd(0,0), -3.0);
+ test_value(ndd(0,1), -5.5);
+
+ // Test operator+(GNdarray,GNdarray)
+ ndd = nda + ndb;
+ test_value(ndd(0,0), 6.0);
+ test_value(ndd(0,1), 11.0);
+
+ // Test operator+(GNdarray,double) and operator+(double,GNdarray)
+ ndd = nda + 1.5;
+ test_value(ndd(0,0), 4.5);
+ test_value(ndd(0,1), 7.0);
+ ndd = 1.5 + nda;
+ test_value(ndd(0,0), 4.5);
+ test_value(ndd(0,1), 7.0);
+
+ // Test operator-(GNdarray,GNdarray)
+ ndd = nda - ndb;
+ test_value(ndd(0,0), 0.0);
+ test_value(ndd(0,1), 0.0);
+
+ // Test operator-(GNdarray,double) and operator-(double,GNdarray)
+ ndd = nda - 1.5;
+ test_value(ndd(0,0), 1.5);
+ test_value(ndd(0,1), 4.0);
+ ndd = 1.5 - nda;
+ test_value(ndd(0,0), -1.5);
+ test_value(ndd(0,1), -4.0);
+
+ // Test operator*(GNdarray,double) and operator*(double,GNdarray)
+ ndd = nda * 3.0;
+ test_value(ndd(0,0), 9.0);
+ test_value(ndd(0,1), 16.5);
+ ndd = 3.0 * nda;
+ test_value(ndd(0,0), 9.0);
+ test_value(ndd(0,1), 16.5);
+
+ // Test operator/(GNdarray,double)
+ ndd = nda / 3.0;
+ test_value(ndd(0,0), 1.0);
+ test_value(ndd(0,1), 5.5/3.0);
+
+ // Test friend functions
+ test_value(min(nda), 3.0);
+ test_value(max(nda), 5.5);
+ test_value(sum(nda), 8.5);
+ test_value(acos(nda)(0,0), std::acos(3.0));
+ test_value(acos(nda)(0,1), std::acos(5.5));
+ test_value(acosh(nda)(0,0), acosh(3.0));
+ test_value(acosh(nda)(0,1), acosh(5.5));
+ test_value(asin(nda)(0,0), std::asin(3.0));
+ test_value(asin(nda)(0,1), std::asin(5.5));
+ test_value(asinh(nda)(0,0), asinh(3.0));
+ test_value(asinh(nda)(0,1), asinh(5.5));
+ test_value(atan(nda)(0,0), std::atan(3.0));
+ test_value(atan(nda)(0,1), std::atan(5.5));
+ test_value(atanh(nda)(0,0), atanh(3.0));
+ test_value(atanh(nda)(0,1), atanh(5.5));
+ test_value(cos(nda)(0,0), std::cos(3.0));
+ test_value(cos(nda)(0,1), std::cos(5.5));
+ test_value(cosh(nda)(0,0), std::cosh(3.0));
+ test_value(cosh(nda)(0,1), std::cosh(5.5));
+ test_value(exp(nda)(0,0), std::exp(3.0));
+ test_value(exp(nda)(0,1), std::exp(5.5));
+ test_value(abs(nda)(0,0), std::abs(3.0));
+ test_value(abs(nda)(0,1), std::abs(5.5));
+ test_value(log(nda)(0,0), std::log(3.0));
+ test_value(log(nda)(0,1), std::log(5.5));
+ test_value(log10(nda)(0,0), std::log10(3.0));
+ test_value(log10(nda)(0,1), std::log10(5.5));
+ test_value(sin(nda)(0,0), std::sin(3.0));
+ test_value(sin(nda)(0,1), std::sin(5.5));
+ test_value(sinh(nda)(0,0), std::sinh(3.0));
+ test_value(sinh(nda)(0,1), std::sinh(5.5));
+ test_value(sqrt(nda)(0,0), std::sqrt(3.0));
+ test_value(sqrt(nda)(0,1), std::sqrt(5.5));
+ test_value(tan(nda)(0,0), std::tan(3.0));
+ test_value(tan(nda)(0,1), std::tan(5.5));
+ test_value(tanh(nda)(0,0), std::tanh(3.0));
+ test_value(tanh(nda)(0,1), std::tanh(5.5));
+ test_value(pow(nda,3.0)(0,0), std::pow(3.0,3.0));
+ test_value(pow(nda,3.0)(0,1), std::pow(5.5,3.0));
+
+ // Exit test
+ return;
+}
+
+
+/***********************************************************************//**
* @brief Test model parameter handling
***************************************************************************/
void TestGNumerics::test_integral(void)
diff --git a/test/test_GNumerics.hpp b/test/test_GNumerics.hpp
index e348782..4e82bba 100644
--- a/test/test_GNumerics.hpp
+++ b/test/test_GNumerics.hpp
@@ -1,7 +1,7 @@
/***************************************************************************
* test_GNumerics.hpp - test numerics modules *
* ----------------------------------------------------------------------- *
- * copyright (C) 2010-2014 by Juergen Knoedlseder *
+ * copyright (C) 2010-2016 by Juergen Knoedlseder *
* ----------------------------------------------------------------------- *
* *
* This program is free software: you can redistribute it and/or modify *
@@ -68,6 +68,7 @@ public:
virtual void set(void);
virtual TestGNumerics* clone(void) const;
virtual std::string classname(void) const { return "TestGNumerics"; }
+ void test_ndarray(void);
void test_integral(void);
void test_romberg_integration(void);
void test_adaptive_simpson_integration(void);