Michael Mayer / gammalib

Commit 32eaa5e9aed263014ebe02a551172f44e6c8f086

Authored by Jürgen Knödlseder
1 parent e654f290

Add GNdarray class (#1768) 

An initial implementation of the GNdarray class has been added that implements a n-dimensional double precision floating point array.
Inline Side-by-side
Showing 12 changed files with 2250 additions and 14 deletions
ChangeLog
  View file @32eaa5e
1   -2016-10-07
  1 +2016-10-08
2 2  
3 3 * Version 1.2.0 released
4 4 =========================
5 5  
  6 + Add GNdarray class (#1768)
6 7 Increase precision of gammalib::str() floating point functions
7 8 Add GFits::publish() methods (#1231)
8 9 Add GVOTable class (#1231)
... ...
1 1 New Features and Important Changes in GammaLib 1.2.0
2 2  
3   -7 October 2016
  3 +8 October 2016
4 4  
5 5  
6 6 1. Introduction
... ... @@ -16,6 +16,7 @@ The following classes have been added:
16 16 - GModelSpectralPlawEnergyFlux
17 17 - GModelSpectralMultiplicative
18 18 - GVOTable
  19 +- GNdarray
19 20  
20 21 The following classes have been removed:
21 22 - ...
... ... @@ -155,7 +156,8 @@ product of spectral model components as a new spectral model (#1861).
155 156  
156 157 9. Numerics module
157 158 ------------------
158   -None
  159 +The GNdarray class has been added that implements a n-dimensional double
  160 +precision floating point array (#1768).
159 161  
160 162  
161 163 10. Observation module
... ...
include/GNdarray.hpp 0 → 100644
  View file @32eaa5e
  1 +/***************************************************************************
  2 + * GNdarray.hpp - N-dimensional array class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GNdarray.hpp
  23 + * @brief N-dimensional array class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GNDARRAY_HPP
  28 +#define GNDARRAY_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include <vector>
  33 +#include "GBase.hpp"
  34 +
  35 +
  36 +/***********************************************************************//**
  37 + * @class GNdarray
  38 + *
  39 + * @brief N-dimensional array class
  40 + *
  41 + * This class implement a n-dimensional double precision floating point
  42 + * array.
  43 + ***************************************************************************/
  44 +class GNdarray : public GBase {
  45 +
  46 + // Friend functions
  47 + friend double min(const GNdarray& array);
  48 + friend double max(const GNdarray& array);
  49 + friend double sum(const GNdarray& array);
  50 + friend GNdarray acos(const GNdarray& array);
  51 + friend GNdarray acosh(const GNdarray& array);
  52 + friend GNdarray asin(const GNdarray& array);
  53 + friend GNdarray asinh(const GNdarray& array);
  54 + friend GNdarray atan(const GNdarray& array);
  55 + friend GNdarray atanh(const GNdarray& array);
  56 + friend GNdarray cos(const GNdarray& array);
  57 + friend GNdarray cosh(const GNdarray& array);
  58 + friend GNdarray exp(const GNdarray& array);
  59 + friend GNdarray abs(const GNdarray& array);
  60 + friend GNdarray log(const GNdarray& array);
  61 + friend GNdarray log10(const GNdarray& array);
  62 + friend GNdarray sin(const GNdarray& array);
  63 + friend GNdarray sinh(const GNdarray& array);
  64 + friend GNdarray sqrt(const GNdarray& array);
  65 + friend GNdarray tan(const GNdarray& array);
  66 + friend GNdarray tanh(const GNdarray& array);
  67 + friend GNdarray pow(const GNdarray& array, const double& power);
  68 +
  69 +public:
  70 + // Constructors and destructors
  71 + GNdarray(void);
  72 + explicit GNdarray(const int& nx);
  73 + GNdarray(const int& nx, const int& ny);
  74 + GNdarray(const int& nx, const int& ny, const int& nz);
  75 + GNdarray(const std::vector<int>& n);
  76 + GNdarray(const GNdarray& array);
  77 + virtual ~GNdarray(void);
  78 +
  79 + // Element access operators
  80 + double& operator()(const int& ix);
  81 + double& operator()(const int& ix, const int& iy);
  82 + double& operator()(const int& ix, const int& iy, const int& iz);
  83 + double& operator()(const std::vector<int>& i);
  84 + const double& operator()(const int& ix) const;
  85 + const double& operator()(const int& ix, const int& iy) const;
  86 + const double& operator()(const int& ix, const int& iy, const int& iz) const;
  87 + const double& operator()(const std::vector<int>& i) const;
  88 +
  89 + // Operators
  90 + GNdarray& operator=(const GNdarray& array);
  91 + bool operator==(const GNdarray& array) const;
  92 + bool operator!=(const GNdarray& array) const;
  93 + GNdarray& operator+=(const GNdarray& array);
  94 + GNdarray& operator-=(const GNdarray& array);
  95 + GNdarray& operator+=(const double& value);
  96 + GNdarray& operator-=(const double& value);
  97 + GNdarray& operator*=(const double& value);
  98 + GNdarray& operator/=(const double& value);
  99 + GNdarray operator-(void) const;
  100 +
  101 + // Methods
  102 + void clear(void);
  103 + GNdarray* clone(void) const;
  104 + std::string classname(void) const;
  105 + int dim() const;
  106 + int size(void) const;
  107 + const std::vector<int>& shape(void) const;
  108 + double& at(const int& ix);
  109 + double& at(const int& ix, const int& iy);
  110 + double& at(const int& ix, const int& iy, const int& iz);
  111 + double& at(const std::vector<int>& i);
  112 + const double& at(const int& ix) const;
  113 + const double& at(const int& ix, const int& iy) const;
  114 + const double& at(const int& ix, const int& iy, const int& iz) const;
  115 + const double& at(const std::vector<int>& i) const;
  116 + std::string print(const GChatter& chatter = NORMAL) const;
  117 +
  118 +#ifndef SWIG
  119 +protected:
  120 +#endif
  121 + // Protected methods
  122 + void init_members(void);
  123 + void copy_members(const GNdarray& array);
  124 + void free_members(void);
  125 + int index(const std::vector<int>& i) const;
  126 + bool has_same_shape(const GNdarray& array) const;
  127 + void require_same_shape(const std::string& method, const GNdarray& array) const;
  128 +
  129 + // Protected members
  130 + std::vector<int> m_shape; //!< Array dimensions
  131 + std::vector<int> m_strides; //!< Steps in each dimension when traversing array
  132 + std::vector<double> m_data; //!< Array data
  133 +};
  134 +
  135 +
  136 +/***********************************************************************//**
  137 + * @brief Return class name
  138 + *
  139 + * @return String containing the class name ("GNdarray").
  140 + ***************************************************************************/
  141 +inline
  142 +std::string GNdarray::classname(void) const
  143 +{
  144 + return ("GNdarray");
  145 +}
  146 +
  147 +
  148 +/***********************************************************************//**
  149 + * @brief 1-dimensional array element access operator
  150 + *
  151 + * @param[in] ix Element index [0,...,shape(0)-1].
  152 + * @return Reference to array element.
  153 + ***************************************************************************/
  154 +inline
  155 +double& GNdarray::operator()(const int& ix)
  156 +{
  157 + // Return array element
  158 + return m_data[ix];
  159 +}
  160 +
  161 +
  162 +/***********************************************************************//**
  163 + * @brief 2-dimensional array element access operator
  164 + *
  165 + * @param[in] ix Index in first dimension [0,...,shape(0)-1].
  166 + * @param[in] iy Index in second dimension [0,...,shape(1)-1].
  167 + * @return Reference to array element.
  168 + ***************************************************************************/
  169 +inline
  170 +double& GNdarray::operator()(const int& ix, const int& iy)
  171 +{
  172 + // Return array element
  173 + return m_data[ix+m_strides[1]*iy];
  174 +}
  175 +
  176 +
  177 +/***********************************************************************//**
  178 + * @brief 3-dimensional array element access operator
  179 + *
  180 + * @param[in] ix Index in first dimension [0,...,shape(0)-1].
  181 + * @param[in] iy Index in second dimension [0,...,shape(1)-1].
  182 + * @param[in] iz Index in third dimension [0,...,shape(2)-1].
  183 + * @return Reference to array element.
  184 + ***************************************************************************/
  185 +inline
  186 +double& GNdarray::operator()(const int& ix, const int& iy, const int& iz)
  187 +{
  188 + // Return array element
  189 + return m_data[ix+m_strides[1]*iy+m_strides[2]*iz];
  190 +}
  191 +
  192 +
  193 +/***********************************************************************//**
  194 + * @brief n-dimensional array element access operator
  195 + *
  196 + * @param[in] i Index vector.
  197 + * @return Reference to array element.
  198 + ***************************************************************************/
  199 +inline
  200 +double& GNdarray::operator()(const std::vector<int>& i)
  201 +{
  202 + // Return array element
  203 + return m_data[index(i)];
  204 +}
  205 +
  206 +
  207 +/***********************************************************************//**
  208 + * @brief 1-dimensional array element access operator (const variant)
  209 + *
  210 + * @param[in] ix Element index [0,...,shape(0)-1].
  211 + * @return Const reference to array element.
  212 + ***************************************************************************/
  213 +inline
  214 +const double& GNdarray::operator()(const int& ix) const
  215 +{
  216 + // Return array element
  217 + return m_data[ix];
  218 +}
  219 +
  220 +
  221 +/***********************************************************************//**
  222 + * @brief 2-dimensional array element access operator (const variant)
  223 + *
  224 + * @param[in] ix Index in first dimension [0,...,shape(0)-1].
  225 + * @param[in] iy Index in second dimension [0,...,shape(1)-1].
  226 + * @return Const reference to array element.
  227 + ***************************************************************************/
  228 +inline
  229 +const double& GNdarray::operator()(const int& ix, const int& iy) const
  230 +{
  231 + // Return array element
  232 + return m_data[ix+m_strides[1]*iy];
  233 +}
  234 +
  235 +
  236 +/***********************************************************************//**
  237 + * @brief 3-dimensional array element access operator (const variant)
  238 + *
  239 + * @param[in] ix Index in first dimension [0,...,shape(0)-1].
  240 + * @param[in] iy Index in second dimension [0,...,shape(1)-1].
  241 + * @param[in] iz Index in third dimension [0,...,shape(2)-1].
  242 + * @return Const reference to array element.
  243 + ***************************************************************************/
  244 +inline
  245 +const double& GNdarray::operator()(const int& ix, const int& iy, const int& iz) const
  246 +{
  247 + // Return array element
  248 + return m_data[ix+m_strides[1]*iy+m_strides[2]*iz];
  249 +}
  250 +
  251 +
  252 +/***********************************************************************//**
  253 + * @brief n-dimensional array element access operator (const variant)
  254 + *
  255 + * @param[in] i Index vector.
  256 + * @return Const reference to array element.
  257 + ***************************************************************************/
  258 +inline
  259 +const double& GNdarray::operator()(const std::vector<int>& i) const
  260 +{
  261 + // Return array element
  262 + return m_data[index(i)];
  263 +}
  264 +
  265 +
  266 +/***********************************************************************//**
  267 + * @brief Return dimension of array
  268 + *
  269 + * @return Dimension of array
  270 + *
  271 + * Returns the dimension of the array.
  272 + ***************************************************************************/
  273 +inline
  274 +int GNdarray::dim() const
  275 +{
  276 + return m_shape.size();
  277 +}
  278 +
  279 +
  280 +/***********************************************************************//**
  281 + * @brief Return number of elements in array
  282 + *
  283 + * @return Number of elements in array
  284 + *
  285 + * Returns the number of elements in the array.
  286 + ***************************************************************************/
  287 +inline
  288 +int GNdarray::size() const
  289 +{
  290 + return m_data.size();
  291 +}
  292 +
  293 +
  294 +/***********************************************************************//**
  295 + * @brief Return shape of array
  296 + *
  297 + * @return Shape of array
  298 + *
  299 + * Returns the shape of the array.
  300 + ***************************************************************************/
  301 +inline
  302 +const std::vector<int>& GNdarray::shape(void) const
  303 +{
  304 + return m_shape;
  305 +}
  306 +
  307 +
  308 +/***********************************************************************//**
  309 + * @brief 1-dimensional array element access with range checking
  310 + *
  311 + * @param[in] ix Element index [0,...,shape(0)-1].
  312 + * @return Reference to array element.
  313 + ***************************************************************************/
  314 +inline
  315 +double& GNdarray::at(const int& ix)
  316 +{
  317 + return const_cast<double &>(static_cast<const GNdarray &>(*this).at(ix));
  318 +}
  319 +
  320 +
  321 +/***********************************************************************//**
  322 + * @brief 2-dimensional array element access with range checking
  323 + *
  324 + * @param[in] ix Index in first dimension [0,...,shape(0)-1].
  325 + * @param[in] iy Index in second dimension [0,...,shape(1)-1].
  326 + * @return Reference to array element.
  327 + ***************************************************************************/
  328 +inline
  329 +double& GNdarray::at(const int& ix, const int& iy)
  330 +{
  331 + return const_cast<double &>(static_cast<const GNdarray &>(*this).at(ix,iy));
  332 +}
  333 +
  334 +
  335 +/***********************************************************************//**
  336 + * @brief 3-dimensional array element access with range checking
  337 + *
  338 + * @param[in] ix Index in first dimension [0,...,shape(0)-1].
  339 + * @param[in] iy Index in second dimension [0,...,shape(1)-1].
  340 + * @param[in] iz Index in third dimension [0,...,shape(2)-1].
  341 + * @return Reference to array element.
  342 + ***************************************************************************/
  343 +inline
  344 +double& GNdarray::at(const int& ix, const int& iy, const int& iz)
  345 +{
  346 + return const_cast<double &>(static_cast<const GNdarray &>(*this).at(ix,iy,iz));
  347 +}
  348 +
  349 +
  350 +/***********************************************************************//**
  351 + * @brief n-dimensional array element access with range checking
  352 + *
  353 + * @param[in] i Index vector.
  354 + * @return Reference to array element.
  355 + ***************************************************************************/
  356 +inline
  357 +double& GNdarray::at(const std::vector<int>& i)
  358 +{
  359 + return const_cast<double &>(static_cast<const GNdarray &>(*this).at(i));
  360 +}
  361 +
  362 +
  363 +/***********************************************************************//**
  364 + * @brief Return sum of two arrays
  365 + *
  366 + * @param[in] a First array.
  367 + * @param[in] b Second array.
  368 + * @return Sum of arrays @p a and @p b.
  369 + *
  370 + * Returns the sum of arrays @p a and @p b.
  371 + ***************************************************************************/
  372 +inline
  373 +GNdarray operator+(const GNdarray& a, const GNdarray& b)
  374 +{
  375 + GNdarray result = a;
  376 + result += b;
  377 + return result;
  378 +}
  379 +
  380 +
  381 +/***********************************************************************//**
  382 + * @brief Add value to array (right addition)
  383 + *
  384 + * @param[in] array Array.
  385 + * @param[in] value Value.
  386 + * @return Array with @p value added to all elements.
  387 + *
  388 + * Returns an array for which the @p value has been added to all elements.
  389 + ***************************************************************************/
  390 +inline
  391 +GNdarray operator+(const GNdarray& array, const double& value)
  392 +{
  393 + GNdarray result = array;
  394 + result += value;
  395 + return result;
  396 +}
  397 +
  398 +
  399 +/***********************************************************************//**
  400 + * @brief Add value to array (left addition)
  401 + *
  402 + * @param[in] value Value.
  403 + * @param[in] array Array.
  404 + * @return Array with @p value added to all elements.
  405 + *
  406 + * Returns an array for which the @p value has been added to all elements.
  407 + ***************************************************************************/
  408 +inline
  409 +GNdarray operator+(const double& value, const GNdarray& array)
  410 +{
  411 + GNdarray result = array;
  412 + result += value;
  413 + return result;
  414 +}
  415 +
  416 +
  417 +/***********************************************************************//**
  418 + * @brief Return difference of arrays
  419 + *
  420 + * @param[in] a First array.
  421 + * @param[in] b Second array.
  422 + * @return Difference between array @p a and @p b.
  423 + *
  424 + * Returns the difference between array @p a and @p b.
  425 + ***************************************************************************/
  426 +inline
  427 +GNdarray operator-(const GNdarray& a, const GNdarray& b)
  428 +{
  429 + GNdarray result = a;
  430 + result -= b;
  431 + return result;
  432 +}
  433 +
  434 +
  435 +/***********************************************************************//**
  436 + * @brief Subtract value from array
  437 + *
  438 + * @param[in] array Array.
  439 + * @param[in] value Value.
  440 + * @return Array with @p value subtracted from all elements.
  441 + *
  442 + * Returns an array for which the @p value has been subtracted from all
  443 + * elements. For example
  444 + *
  445 + * double value = 5.0;
  446 + * GNdarray result = array - value;
  447 + ***************************************************************************/
  448 +inline
  449 +GNdarray operator-(const GNdarray& array, const double& value)
  450 +{
  451 + GNdarray result = array;
  452 + result -= value;
  453 + return result;
  454 +}
  455 +
  456 +
  457 +/***********************************************************************//**
  458 + * @brief Subtract array from value
  459 + *
  460 + * @param[in] value Value.
  461 + * @param[in] array Array.
  462 + * @return Array with @p value subtracted from all elements.
  463 + *
  464 + * Returns an array for which all elements have been subtracted from
  465 + * @p value. For example
  466 + *
  467 + * double value = 5.0;
  468 + * GNdarray result = value - array;
  469 + ***************************************************************************/
  470 +inline
  471 +GNdarray operator-(const double& value, const GNdarray& array)
  472 +{
  473 + GNdarray result = -array;
  474 + result += value;
  475 + return result;
  476 +}
  477 +
  478 +
  479 +/***********************************************************************//**
  480 + * @brief Multiply array by value (right multiplication)
  481 + *
  482 + * @param[in] array Array.
  483 + * @param[in] value Value.
  484 + * @return Array for which all elements have be multiplied by @p value.
  485 + *
  486 + * Returns an array for which all elements have be multiplied by @p value.
  487 + ***************************************************************************/
  488 +inline
  489 +GNdarray operator*(const GNdarray& array, const double& value)
  490 +{
  491 + GNdarray result = array;
  492 + result *= value;
  493 + return result;
  494 +}
  495 +
  496 +
  497 +/***********************************************************************//**
  498 + * @brief Multiply array by value (left multiplication)
  499 + *
  500 + * @param[in] value Value.
  501 + * @param[in] array Array.
  502 + * @return Array for which all elements have be multiplied by @p value.
  503 + *
  504 + * Returns an array for which all elements have be multiplied by @p value.
  505 + ***************************************************************************/
  506 +inline
  507 +GNdarray operator*(const double& value, const GNdarray& array)
  508 +{
  509 + GNdarray result = array;
  510 + result *= value;
  511 + return result;
  512 +}
  513 +
  514 +
  515 +/***********************************************************************//**
  516 + * @brief Divide array by value
  517 + *
  518 + * @param[in] array Array.
  519 + * @param[in] value Value.
  520 + * @return Array for which all elements have be divided by @p value.
  521 + *
  522 + * Returns an array for which all elements have be divided by @p value.
  523 + ***************************************************************************/
  524 +inline
  525 +GNdarray operator/(const GNdarray& array, const double& value)
  526 +{
  527 + GNdarray result = array;
  528 + result /= value;
  529 + return result;
  530 +}
  531 +
  532 +#endif /* GNDARRAY_HPP */
... ...
include/GammaLib.hpp
  View file @32eaa5e
... ... @@ -63,6 +63,7 @@
63 63 #include "GDerivative.hpp"
64 64 #include "GFunction.hpp"
65 65 #include "GMath.hpp"
  66 +#include "GNdarray.hpp"
66 67  
67 68 /* __ FITS module ________________________________________________________ */
68 69 #include "GFits.hpp"
... ...
include/Makefile.am
  View file @32eaa5e
... ... @@ -44,6 +44,7 @@ pkginclude_HEADERS = GammaLib.hpp \
44 44 GDerivative.hpp \
45 45 GFunction.hpp \
46 46 GMath.hpp \
  47 + GNdarray.hpp \
47 48 GFits.hpp \
48 49 GFitsHDU.hpp \
49 50 GFitsHeader.hpp \
... ...
pyext/GNdarray.i 0 → 100755
  View file @32eaa5e
  1 +/***************************************************************************
  2 + * GNdarray.i - N-dimensional array class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GNdarray.i
  23 + * @brief N-dimensional array class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +%{
  28 +/* Put headers and other declarations here that are needed for compilation */
  29 +#define SWIG
  30 +#include "GNdarray.hpp"
  31 +%}
  32 +
  33 +
  34 +/***********************************************************************//**
  35 + * @class GNdarray
  36 + *
  37 + * @brief N-dimensional array class
  38 + *
  39 + * This class implement a n-dimensional double precision floating point
  40 + * array.
  41 + ***************************************************************************/
  42 +class GNdarray : public GBase {
  43 +
  44 +public:
  45 + // Constructors and destructors
  46 + GNdarray(void);
  47 + explicit GNdarray(const int& nx);
  48 + GNdarray(const int& nx, const int& ny);
  49 + GNdarray(const int& nx, const int& ny, const int& nz);
  50 + GNdarray(const std::vector<int>& n);
  51 + GNdarray(const GNdarray& array);
  52 + virtual ~GNdarray(void);
  53 +
  54 + // Operators
  55 + bool operator==(const GNdarray& array) const;
  56 + bool operator!=(const GNdarray& array) const;
  57 + GNdarray& operator+=(const GNdarray& array);
  58 + GNdarray& operator-=(const GNdarray& array);
  59 + GNdarray& operator+=(const double& value);
  60 + GNdarray& operator-=(const double& value);
  61 + GNdarray& operator*=(const double& value);
  62 + GNdarray operator-(void) const;
  63 +
  64 + // Methods
  65 + void clear(void);
  66 + GNdarray* clone(void) const;
  67 + std::string classname(void) const;
  68 + int dim() const;
  69 + int size(void) const;
  70 + const std::vector<int>& shape(void) const;
  71 +};
  72 +
  73 +
  74 +/***********************************************************************//**
  75 + * @brief GNdarray class extension
  76 + ***************************************************************************/
  77 +%extend GNdarray {
  78 + double __getitem__(int GTuple[]) {
  79 + // Handle first the special case of a single index. In Python we want
  80 + // to support accessing a n-dimensional array with a single index so
  81 + // that we can iterate over the array
  82 + if (GTuple[0] == 1) {
  83 + if (GTuple[1] >= 0 && GTuple[1] < self->size()) {
  84 + return (self->m_data[GTuple[1]]);
  85 + }
  86 + else {
  87 + throw GException::out_of_range("__getitem__(int)", "Array "
  88 + "index", GTuple[1], self->size());
  89 + }
  90 + }
  91 + if (GTuple[0] != self->dim()) {
  92 + throw GException::invalid_value("__getitem__(int)", "Invalid "
  93 + "access of "+gammalib::str(self->dim())+"-dimensional "
  94 + "array with "+gammalib::str(GTuple[0])+"-dimensional "
  95 + "access operator.");
  96 + }
  97 + int index = 0;
  98 + for (int i = 0; i < self->shape().size(); ++i) {
  99 + if (GTuple[i+1] < 0 || GTuple[i+1] >= self->shape()[i]) {
  100 + throw GException::out_of_range("__getitem__(int)", "Dimension "+
  101 + gammalib::str(i)+" index", GTuple[i+1], self->shape()[i]);
  102 + }
  103 + index += self->m_strides[i] * GTuple[i+1];
  104 + }
  105 + return (self->m_data[index]);
  106 + }
  107 + void __setitem__(int GTuple[], const double& value) {
  108 + if (GTuple[0] != self->dim()) {
  109 + throw GException::invalid_value("__setitem__(int)", "Invalid "
  110 + "access of "+gammalib::str(self->dim())+"-dimensional "
  111 + "array with "+gammalib::str(GTuple[0])+"-dimensional "
  112 + "access operator.");
  113 + }
  114 + int index = 0;
  115 + for (int i = 0; i < self->shape().size(); ++i) {
  116 + if (GTuple[i+1] < 0 || GTuple[i+1] >= self->shape()[i]) {
  117 + throw GException::out_of_range("__setitem__(int)", "Dimension "+
  118 + gammalib::str(i)+" index", GTuple[i+1], self->shape()[i]);
  119 + }
  120 + index += self->m_strides[i] * GTuple[i+1];
  121 + }
  122 + self->m_data[index] = value;
  123 + }
  124 + GNdarray __add__(const GNdarray &a) {
  125 + return (*self) + a;
  126 + }
  127 + GNdarray __add__(const double &a) {
  128 + return (*self) + a;
  129 + }
  130 + GNdarray __sub__(const GNdarray &a) {
  131 + return (*self) - a;
  132 + }
  133 + GNdarray __sub__(const double &a) {
  134 + return (*self) - a;
  135 + }
  136 + GNdarray __mul__(const double &a) {
  137 + return (*self) * a;
  138 + }
  139 + // Python 2.x
  140 + GNdarray __div__(const double &a) {
  141 + return (*self) / a;
  142 + }
  143 + // Python 3.x
  144 + GNdarray __truediv__(const double &a) {
  145 + return (*self) / a;
  146 + }
  147 + // Python 2.x operator/=
  148 + GNdarray __idiv__(const double& value) {
  149 + self->operator/=(value);
  150 + return (*self);
  151 + }
  152 + // Python 3.x operator/=
  153 + GNdarray __itruediv__(const double& value) {
  154 + self->operator/=(value);
  155 + return (*self);
  156 + }
  157 + int __is__(const GNdarray &a) {
  158 + return (*self) == a;
  159 + }
  160 + GNdarray copy() {
  161 + return (*self);
  162 + }
  163 + double min() {
  164 + return min(*self);
  165 + }
  166 + double max() {
  167 + return max(*self);
  168 + }
  169 + double sum() {
  170 + return sum(*self);
  171 + }
  172 + GNdarray acos() {
  173 + return acos(*self);
  174 + }
  175 + GNdarray acosh() {
  176 + return acosh(*self);
  177 + }
  178 + GNdarray asin() {
  179 + return asin(*self);
  180 + }
  181 + GNdarray asinh() {
  182 + return asinh(*self);
  183 + }
  184 + GNdarray atan() {
  185 + return atan(*self);
  186 + }
  187 + GNdarray atanh() {
  188 + return atanh(*self);
  189 + }
  190 + GNdarray cos() {
  191 + return cos(*self);
  192 + }
  193 + GNdarray cosh() {
  194 + return cosh(*self);
  195 + }
  196 + GNdarray exp() {
  197 + return exp(*self);
  198 + }
  199 + GNdarray abs() {
  200 + return abs(*self);
  201 + }
  202 + GNdarray log() {
  203 + return log(*self);
  204 + }
  205 + GNdarray log10() {
  206 + return log10(*self);
  207 + }
  208 + GNdarray sin() {
  209 + return sin(*self);
  210 + }
  211 + GNdarray sinh() {
  212 + return sinh(*self);
  213 + }
  214 + GNdarray sqrt() {
  215 + return sqrt(*self);
  216 + }
  217 + GNdarray tan() {
  218 + return tan(*self);
  219 + }
  220 + GNdarray tanh() {
  221 + return tanh(*self);
  222 + }
  223 +};
... ...
pyext/GTypemaps.i
  View file @32eaa5e
... ... @@ -78,17 +78,17 @@
78 78  
79 79  
80 80 /***********************************************************************//**
81   - * @brief Tuple to index conversion using variable dimensions (1,...,4)
  81 + * @brief Tuple to index conversion using variable dimensions
82 82 *
83 83 * The following function provides conversion between a Python tuple and
84 84 * an integer array. This allows index access via tuples, such as in
85 85 * a[3,5,10] = 10.0 or c = a[2,9].
86 86 ***************************************************************************/
87 87 %{
88   -static int var_tuple_to_index(PyObject *input, int *ptr) {
  88 +static int var_tuple_to_index(PyObject *input, int *ptr, int dim) {
89 89 if (PySequence_Check(input)) {
90 90 int size = PyObject_Length(input);
91   - if (size > 4) {
  91 + if (size > dim) {
92 92 PyErr_SetString(PyExc_ValueError,"Too many arguments in tuple");
93 93 return 0;
94 94 }
... ... @@ -116,8 +116,8 @@ static int var_tuple_to_index(PyObject *input, int *ptr) {
116 116 }
117 117 }
118 118 %}
119   -%typemap(in) int GTuple[ANY] (int temp[5]) {
120   - if (!var_tuple_to_index($input,temp)) {
  119 +%typemap(in) int GTuple[ANY] (int temp[11]) {
  120 + if (!var_tuple_to_index($input,temp,10)) {
121 121 return NULL;
122 122 }
123 123 $1 = &temp[0];
... ...
pyext/gammalib/numerics.i
  View file @32eaa5e
... ... @@ -38,8 +38,11 @@
38 38  
39 39 /* __ Include standard typemaps for vectors and strings __________________ */
40 40 %include stl.i
  41 +%include std_vector.i
  42 +%template(vectori) std::vector<int>;
41 43  
42 44 /* __ Include GammaLib typemaps __________________________________________ */
  45 +%include GTypemaps.i
43 46 %include typemap_GFilename.i
44 47  
45 48 /* __ Include interface classes __________________________________________ */
... ... @@ -53,3 +56,4 @@
53 56 %include "GFunction.i"
54 57 %include "GIntegral.i"
55 58 %include "GMath.i"
  59 +%include "GNdarray.i"
... ...
src/numerics/GNdarray.cpp 0 → 100644
  View file @32eaa5e
  1 +/***************************************************************************
  2 + * GNdarray.cpp - N-dimensional array class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GNdarray.cpp
  23 + * @brief N-dimensional array class implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GTools.hpp"
  32 +#include "GException.hpp"
  33 +#include "GNdarray.hpp"
  34 +
  35 +
  36 +/* __ Method name definitions ____________________________________________ */
  37 +#define G_OP_ADD "GNdarray::operator+=(GNdarray&)"
  38 +#define G_OP_SUB "GNdarray::operator-=(GNdarray&)"
  39 +#define G_AT1 "GNdarray::at(int&)"
  40 +#define G_AT2 "GNdarray::at(int&, int&)"
  41 +#define G_AT3 "GNdarray::at(int&, int&, int&)"
  42 +#define G_ATN "GNdarray::at(std::vector<int>&)"
  43 +
  44 +
  45 +/*==========================================================================
  46 + = =
  47 + = Constructors/destructors =
  48 + = =
  49 + ==========================================================================*/
  50 +
  51 +/***********************************************************************//**
  52 + * @brief Void constructor
  53 + ***************************************************************************/
  54 +GNdarray::GNdarray(void)
  55 +{
  56 + // Initialise class members
  57 + init_members();
  58 +
  59 + // Return
  60 + return;
  61 +}
  62 +
  63 +
  64 +/***********************************************************************//**
  65 + * @brief 1-dimensional constructor
  66 + *
  67 + * @param[in] nx Array dimension.
  68 + *
  69 + * Constructs a 1-dimensional array.
  70 + ***************************************************************************/
  71 +GNdarray::GNdarray(const int& nx)
  72 +{
  73 + // Initialise class members
  74 + init_members();
  75 +
  76 + // Set shape
  77 + m_shape.push_back(nx);
  78 +
  79 + // Set strides
  80 + m_strides.push_back(1);
  81 +
  82 + // Initialise data array with zeros
  83 + m_data.assign(nx, 0.0);
  84 +
  85 + // Return
  86 + return;
  87 +}
  88 +
  89 +
  90 +/***********************************************************************//**
  91 + * @brief 2-dimensional constructor
  92 + *
  93 + * @param[in] nx Array's first dimension.
  94 + * @param[in] ny Array's second dimension.
  95 + *
  96 + * Constructs a 2-dimensional array.
  97 + ***************************************************************************/
  98 +GNdarray::GNdarray(const int& nx, const int& ny)
  99 +{
  100 + // Initialise class members
  101 + init_members();
  102 +
  103 + // Compute size of array
  104 + int size = nx * ny;
  105 +
  106 + // Set shape
  107 + m_shape.push_back(nx);
  108 + m_shape.push_back(ny);
  109 +
  110 + // Set strides
  111 + m_strides.push_back(1);
  112 + m_strides.push_back(nx);
  113 +
  114 + // Initialise data array with zeros
  115 + m_data.assign(size, 0.0);
  116 +
  117 + // Return
  118 + return;
  119 +}
  120 +
  121 +
  122 +/***********************************************************************//**
  123 + * @brief 3-dimensional constructor
  124 + *
  125 + * @param[in] nx Array's first dimension.
  126 + * @param[in] ny Array's second dimension.
  127 + * @param[in] nz Array's third dimension.
  128 + *
  129 + * Constructs a 3-dimensional array.
  130 + ***************************************************************************/
  131 +GNdarray::GNdarray(const int& nx, const int& ny, const int& nz)
  132 +{
  133 + // Initialise class members
  134 + init_members();
  135 +
  136 + // Compute size of array
  137 + int size = nx * ny * nz;
  138 +
  139 + // Set shape
  140 + m_shape.push_back(nx);
  141 + m_shape.push_back(ny);
  142 + m_shape.push_back(nz);
  143 +
  144 + // Set strides
  145 + m_strides.push_back(1);
  146 + m_strides.push_back(nx);
  147 + m_strides.push_back(nx*ny);
  148 +
  149 + // Initialise data array with zeros
  150 + m_data.assign(size, 0.0);
  151 +
  152 + // Return
  153 + return;
  154 +}
  155 +
  156 +
  157 +/***********************************************************************//**
  158 + * @brief n-dimensional constructor
  159 + *
  160 + * @param[in] n Array dimensions.
  161 + *
  162 + * Constructs a n-dimensional array.
  163 + ***************************************************************************/
  164 +GNdarray::GNdarray(const std::vector<int>& n)
  165 +{
  166 + // Initialise class members
  167 + init_members();
  168 +
  169 + // Compute size of array
  170 + int size = n[0];
  171 + for (int i = 1; i < n.size(); ++i) {
  172 + size *= n[i];
  173 + }
  174 +
  175 + // Set shape
  176 + m_shape = n;
  177 +
  178 + // Set strides
  179 + int stride = 1;
  180 + for (int i = 0; i < n.size(); ++i) {
  181 + m_strides.push_back(stride);
  182 + stride *= n[i];
  183 + }
  184 +
  185 + // Initialise data array with zeros
  186 + m_data.assign(size, 0.0);
  187 +
  188 + // Return
  189 + return;
  190 +}
  191 +
  192 +
  193 +/***********************************************************************//**
  194 + * @brief Copy constructor
  195 + *
  196 + * @param[in] array Array.
  197 + ***************************************************************************/
  198 +GNdarray::GNdarray(const GNdarray& array)
  199 +{
  200 + // Initialise class members
  201 + init_members();
  202 +
  203 + // Copy members
  204 + copy_members(array);
  205 +
  206 + // Return
  207 + return;
  208 +}
  209 +
  210 +
  211 +/***********************************************************************//**
  212 + * @brief Destructor
  213 + ***************************************************************************/
  214 +GNdarray::~GNdarray(void)
  215 +{
  216 + // Free members
  217 + free_members();
  218 +
  219 + // Return
  220 + return;
  221 +}
  222 +
  223 +
  224 +/*==========================================================================
  225 + = =
  226 + = Operators =
  227 + = =
  228 + ==========================================================================*/
  229 +
  230 +/***********************************************************************//**
  231 + * @brief Assignment operator
  232 + *
  233 + * @param[in] array Array.
  234 + * @return Array.
  235 + ***************************************************************************/
  236 +GNdarray& GNdarray::operator=(const GNdarray& array)
  237 +{
  238 + // Execute only if object is not identical
  239 + if (this != &array) {
  240 +
  241 + // Free members
  242 + free_members();
  243 +
  244 + // Initialise private members
  245 + init_members();
  246 +
  247 + // Copy members
  248 + copy_members(array);
  249 +
  250 + } // endif: object was not identical
  251 +
  252 + // Return this object
  253 + return *this;
  254 +}
  255 +
  256 +
  257 +/***********************************************************************//**
  258 + * @brief Equality operator
  259 + *
  260 + * @param[in] array Array.
  261 + * @return True if arrays are identical.
  262 + *
  263 + * Returns true if the arrays have identical shape and content.
  264 + ***************************************************************************/
  265 +bool GNdarray::operator==(const GNdarray& array) const
  266 +{
  267 + // Check of arrays have the same shape
  268 + bool identity = has_same_shape(array);
  269 +
  270 + // If the arrays have identical dimensions then check if all elements
  271 + // are identical. Break as soon as one array element differs.
  272 + if (identity) {
  273 + for (int i = 0; i < m_data.size(); ++i) {
  274 + if (m_data[i] != array.m_data[i]) {
  275 + identity = false;
  276 + break;
  277 + }
  278 + }
  279 + }
  280 +
  281 + // Return result
  282 + return identity;
  283 +}
  284 +
  285 +
  286 +/***********************************************************************//**
  287 + * @brief Non-equality operator
  288 + *
  289 + * @param[in] array Array.
  290 + * @return True if arrays are different.
  291 + *
  292 + * Returns true if the arrays have different shape or content.
  293 + ***************************************************************************/
  294 +bool GNdarray::operator!=(const GNdarray& array) const
  295 +{
  296 + // Get negated result of equality operation
  297 + bool difference = !(this->operator==(array));
  298 +
  299 + // Return result
  300 + return difference;
  301 +}
  302 +
  303 +
  304 +/***********************************************************************//**
  305 + * @brief Unary addition operator
  306 + *
  307 + * @param[in] array Array.
  308 + * @return Array.
  309 + *
  310 + * Adds an array to the current array.
  311 + ***************************************************************************/
  312 +GNdarray& GNdarray::operator+=(const GNdarray& array)
  313 +{
  314 + // Throw an exception if the arrays have not the same shape
  315 + require_same_shape(G_OP_ADD, array);
  316 +
  317 + // Add elements
  318 + for (int i = 0; i < m_data.size(); ++i) {
  319 + m_data[i] += array.m_data[i];
  320 + }
  321 +
  322 + // Return array
  323 + return *this;
  324 +}
  325 +
  326 +
  327 +/***********************************************************************//**
  328 + * @brief Unary subtraction operator
  329 + *
  330 + * @param[in] array Array.
  331 + * @return Array.
  332 + *
  333 + * Subtracts an array from the current array.
  334 + ***************************************************************************/
  335 +GNdarray& GNdarray::operator-=(const GNdarray& array)
  336 +{
  337 + // Throw an exception if the arrays have not the same shape
  338 + require_same_shape(G_OP_SUB, array);
  339 +
  340 + // Subtract elements
  341 + for (int i = 0; i < m_data.size(); ++i) {
  342 + m_data[i] -= array.m_data[i];
  343 + }
  344 +
  345 + // Return array
  346 + return *this;
  347 +}
  348 +
  349 +
  350 +/***********************************************************************//**
  351 + * @brief Unary value addition operator
  352 + *
  353 + * @param[in] value Value.
  354 + * @return Array.
  355 + *
  356 + * Adds a value to all array elements.
  357 + ***************************************************************************/
  358 +GNdarray& GNdarray::operator+=(const double& value)
  359 +{
  360 + // Add value to elements
  361 + for (int i = 0; i < m_data.size(); ++i) {
  362 + m_data[i] += value;
  363 + }
  364 +
  365 + // Return array
  366 + return *this;
  367 +}
  368 +
  369 +
  370 +/***********************************************************************//**
  371 + * @brief Unary value subtraction operator
  372 + *
  373 + * @param[in] value Value.
  374 + * @return Array.
  375 + *
  376 + * Subtracts a value from all array elements.
  377 + ***************************************************************************/
  378 +GNdarray& GNdarray::operator-=(const double& value)
  379 +{
  380 + // Subtract value from elements
  381 + for (int i = 0; i < m_data.size(); ++i) {
  382 + m_data[i] -= value;
  383 + }
  384 +
  385 + // Return array
  386 + return *this;
  387 +}
  388 +
  389 +
  390 +/***********************************************************************//**
  391 + * @brief Unary multiplication operator
  392 + *
  393 + * @param[in] value Value.
  394 + * @return Array.
  395 + *
  396 + * Multiply all array elements by a value.
  397 + ***************************************************************************/
  398 +GNdarray& GNdarray::operator*=(const double& value)
  399 +{
  400 + // Multiply all elements
  401 + for (int i = 0; i < m_data.size(); ++i) {
  402 + m_data[i] *= value;
  403 + }
  404 +
  405 + // Return array
  406 + return *this;
  407 +}
  408 +
  409 +
  410 +/***********************************************************************//**
  411 + * @brief Unary division operator
  412 + *
  413 + * @param[in] value Value.
  414 + * @return Array.
  415 + *
  416 + * Divide all array elements by a value.
  417 + ***************************************************************************/
  418 +GNdarray& GNdarray::operator/=(const double& value)
  419 +{
  420 + // Divide all elements
  421 + for (int i = 0; i < m_data.size(); ++i) {
  422 + m_data[i] /= value;
  423 + }
  424 +
  425 + // Return array
  426 + return *this;
  427 +}
  428 +
  429 +
  430 +/***********************************************************************//**
  431 + * @brief Unary minus operator
  432 + *
  433 + * @return Array.
  434 + *
  435 + * Negate all array elements.
  436 + ***************************************************************************/
  437 +GNdarray GNdarray::operator-(void) const
  438 +{
  439 + // Copy array
  440 + GNdarray result = *this;
  441 +
  442 + // Negate all elements
  443 + for (int i = 0; i < m_data.size(); ++i) {
  444 + result.m_data[i] = -result.m_data[i];
  445 + }
  446 +
  447 + // Return array
  448 + return result;
  449 +}
  450 +
  451 +
  452 +/*==========================================================================
  453 + = =
  454 + = Public methods =
  455 + = =
  456 + ==========================================================================*/
  457 +
  458 +/***********************************************************************//**
  459 + * @brief Clear array
  460 + ***************************************************************************/
  461 +void GNdarray::clear(void)
  462 +{
  463 + // Free members
  464 + free_members();
  465 +
  466 + // Initialise private members
  467 + init_members();
  468 +
  469 + // Return
  470 + return;
  471 +}
  472 +
  473 +
  474 +/***********************************************************************//**
  475 + * @brief Clone array
  476 + *
  477 + * @return Pointer to deep copy of array.
  478 + ***************************************************************************/
  479 +GNdarray* GNdarray::clone(void) const
  480 +{
  481 + // Clone array
  482 + return new GNdarray(*this);
  483 +}
  484 +
  485 +
  486 +/***********************************************************************//**
  487 + * @brief 1-dimensional array element access with range checking
  488 + * (const version)
  489 + *
  490 + * @param[in] ix Element index [0,...,shape(0)-1].
  491 + * @return Const reference to array element.
  492 + *
  493 + * @exception GException::invalid_value
  494 + * Array is not 1-dimensional.
  495 + * @exception GException::out_of_range
  496 + * Element index is out of range.
  497 + ***************************************************************************/
  498 +const double& GNdarray::at(const int& ix) const
  499 +{
  500 + // Throw an exception if the array is not 1-dimensional
  501 + if (m_shape.size() != 1) {
  502 + std::string msg = "Invalid access of "+
  503 + gammalib::str(m_shape.size())+"-dimensional array "
  504 + "with 1-dimensional access operator.";
  505 + throw GException::invalid_value(G_AT1, msg);
  506 + }
  507 +
  508 + // Throw an exception if the index is outside the valid range
  509 + if (ix < 0 || ix >= m_shape[0]) {
  510 + throw GException::out_of_range(G_AT1, "Array index", ix, m_shape[0]);
  511 + }
  512 +
  513 + // Return array element
  514 + return (*this)(ix);
  515 +}
  516 +
  517 +
  518 +/***********************************************************************//**
  519 + * @brief 2-dimensional array element access with range checking
  520 + * (const version)
  521 + *
  522 + * @param[in] ix Index in first dimension [0,...,shape(0)-1].
  523 + * @param[in] iy Index in second dimension [0,...,shape(1)-1].
  524 + * @return Const reference to array element.
  525 + *
  526 + * @exception GException::invalid_value
  527 + * Array is not 2-dimensional.
  528 + * @exception GException::out_of_range
  529 + * Element index is out of range.
  530 + ***************************************************************************/
  531 +const double& GNdarray::at(const int& ix, const int& iy) const
  532 +{
  533 + // Throw an exception if the array is not 2-dimensional
  534 + if (m_shape.size() != 2) {
  535 + std::string msg = "Invalid access of "+
  536 + gammalib::str(m_shape.size())+"-dimensional array "
  537 + "with 2-dimensional access operator.";
  538 + throw GException::invalid_value(G_AT2, msg);
  539 + }
  540 +
  541 + // Throw an exception if the indices are outside the valid range
  542 + if (ix < 0 || ix >= m_shape[0]) {
  543 + throw GException::out_of_range(G_AT2, "First array index", ix,
  544 + m_shape[0]);
  545 + }
  546 + if (iy < 0 || iy >= m_shape[1]) {
  547 + throw GException::out_of_range(G_AT2, "Second array index", iy,
  548 + m_shape[1]);
  549 + }
  550 +
  551 + // Return array element
  552 + return (*this)(ix,iy);
  553 +}
  554 +
  555 +
  556 +/***********************************************************************//**
  557 + * @brief 3-dimensional array element access with range checking
  558 + * (const version)
  559 + *
  560 + * @param[in] ix Index in first dimension [0,...,shape(0)-1].
  561 + * @param[in] iy Index in second dimension [0,...,shape(1)-1].
  562 + * @param[in] iz Index in third dimension [0,...,shape(2)-1].
  563 + * @return Const reference to array element.
  564 + *
  565 + * @exception GException::invalid_value
  566 + * Array is not 3-dimensional.
  567 + * @exception GException::out_of_range
  568 + * Element index is out of range.
  569 + ***************************************************************************/
  570 +const double& GNdarray::at(const int& ix, const int& iy, const int& iz) const
  571 +{
  572 + // Throw an exception if the array is not 2-dimensional
  573 + if (m_shape.size() != 3) {
  574 + std::string msg = "Invalid access of "+
  575 + gammalib::str(m_shape.size())+"-dimensional array "
  576 + "with 3-dimensional access operator.";
  577 + throw GException::invalid_value(G_AT3, msg);
  578 + }
  579 +
  580 + // Throw an exception if the indices are outside the valid range
  581 + if (ix < 0 || ix >= m_shape[0]) {
  582 + throw GException::out_of_range(G_AT3, "First array index", ix,
  583 + m_shape[0]);
  584 + }
  585 + if (iy < 0 || iy >= m_shape[1]) {
  586 + throw GException::out_of_range(G_AT3, "Second array index", iy,
  587 + m_shape[1]);
  588 + }
  589 + if (iz < 0 || iz >= m_shape[2]) {
  590 + throw GException::out_of_range(G_AT3, "Third array index", iz,
  591 + m_shape[2]);
  592 + }
  593 +
  594 + // Return array element
  595 + return (*this)(ix,iy,iz);
  596 +}
  597 +
  598 +
  599 +/***********************************************************************//**
  600 + * @brief n-dimensional array element access with range checking
  601 + * (const version)
  602 + *
  603 + * @param[in] i Index vector.
  604 + * @return Const reference to array element.
  605 + *
  606 + * @exception GException::invalid_value
  607 + * Array is not of correct dimension.
  608 + * @exception GException::out_of_range
  609 + * Element index is out of range.
  610 + ***************************************************************************/
  611 +const double& GNdarray::at(const std::vector<int>& i) const
  612 +{
  613 + // Throw an exception if the array is not 2-dimensional
  614 + if (m_shape.size() != i.size()) {
  615 + std::string msg = "Invalid access of "+
  616 + gammalib::str(m_shape.size())+"-dimensional array "
  617 + "with "+
  618 + gammalib::str(i.size())+"-dimensional access "
  619 + "operator.";
  620 + throw GException::invalid_value(G_ATN, msg);
  621 + }
  622 +
  623 + // Throw an exception if the indices are outside the valid range
  624 + for (int k = 0; k < m_shape.size(); ++k) {
  625 + if (i[k] < 0 || i[k] >= m_shape[k]) {
  626 + throw GException::out_of_range(G_ATN, "Dimension "+
  627 + gammalib::str(k)+" index",
  628 + i[k], m_shape[k]);
  629 + }
  630 + }
  631 +
  632 + // Return array element
  633 + return (*this)(i);
  634 +}
  635 +
  636 +
  637 +/***********************************************************************//**
  638 + * @brief Print array information
  639 + *
  640 + * @param[in] chatter Chattiness.
  641 + * @return String containing array information.
  642 + ***************************************************************************/
  643 +std::string GNdarray::print(const GChatter& chatter) const
  644 +{
  645 + // Initialise result string
  646 + std::string result = "(";
  647 +
  648 + // Continue only if chatter is not silent
  649 + if (chatter != SILENT) {
  650 +
  651 + // Put all elements in string
  652 + for (int i = 0; i < m_data.size(); ++i) {
  653 + if (i > 0) {
  654 + result += ", ";
  655 + }
  656 + result += gammalib::str(m_data[i]);
  657 + }
  658 +
  659 + // Append closing parentheses
  660 + result += ")";
  661 +
  662 + } // endif: chatter was not silent
  663 +
  664 + // Return result
  665 + return result;
  666 +}
  667 +
  668 +
  669 +/*==========================================================================
  670 + = =
  671 + = Private methods =
  672 + = =
  673 + ==========================================================================*/
  674 +
  675 +/***********************************************************************//**
  676 + * @brief Initialise class members
  677 + ***************************************************************************/
  678 +void GNdarray::init_members(void)
  679 +{
  680 + // Initialise members
  681 + m_shape.clear();
  682 + m_strides.clear();
  683 + m_data.clear();
  684 +
  685 + // Return
  686 + return;
  687 +}
  688 +
  689 +
  690 +/***********************************************************************//**
  691 + * @brief Copy class members
  692 + *
  693 + * @param[in] array Array.
  694 + ***************************************************************************/
  695 +void GNdarray::copy_members(const GNdarray& array)
  696 +{
  697 + // Copy members
  698 + m_shape = array.m_shape;
  699 + m_strides = array.m_strides;
  700 + m_data = array.m_data;
  701 +
  702 + // Return
  703 + return;
  704 +}
  705 +
  706 +
  707 +/***********************************************************************//**
  708 + * @brief Delete class members
  709 + ***************************************************************************/
  710 +void GNdarray::free_members(void)
  711 +{
  712 + // Return
  713 + return;
  714 +}
  715 +
  716 +
  717 +/***********************************************************************//**
  718 + * @brief Compute array element index
  719 + *
  720 + * @param[in] i Index vector.
  721 + * @return Array element index.
  722 + *
  723 + * Computes the array element index from an index vector. The method assumes
  724 + * that the index vector has the same dimension as the array. If this
  725 + * condition is not guaranteed, the condition needs to be checked before
  726 + * calling the method.
  727 + ***************************************************************************/
  728 +int GNdarray::index(const std::vector<int>& i) const
  729 +{
  730 + // Initialise element index
  731 + int index = 0;
  732 +
  733 + // Compute element index
  734 + for (int k = 0; k < m_shape.size(); ++k) {
  735 + index += m_strides[k] * i[k];
  736 + }
  737 +
  738 + // Return index
  739 + return index;
  740 +}
  741 +
  742 +
  743 +/***********************************************************************//**
  744 + * @brief Check if array has the same shape
  745 + *
  746 + * @param[in] array Array.
  747 + * @return True if the array has the same shape.
  748 + ***************************************************************************/
  749 +bool GNdarray::has_same_shape(const GNdarray& array) const
  750 +{
  751 + // Check if the array has the same dimension
  752 + bool identity = m_shape.size() == array.m_shape.size();
  753 +
  754 + // Check if the shape is identical. Break as soon as one shape value
  755 + // differs.
  756 + if (identity) {
  757 + for (int i = 0; i < m_shape.size(); ++i) {
  758 + if (m_shape[i] != array.m_shape[i]) {
  759 + identity = false;
  760 + break;
  761 + }
  762 + }
  763 + }
  764 +
  765 + // Return result
  766 + return identity;
  767 +}
  768 +
  769 +
  770 +/***********************************************************************//**
  771 + * @brief Throw exception if array shapes differ
  772 + *
  773 + * @param[in] array Array.
  774 + ***************************************************************************/
  775 +void GNdarray::require_same_shape(const std::string& method,
  776 + const GNdarray& array) const
  777 +{
  778 + // If the shape differs then throw an exception
  779 + if (!has_same_shape(array)) {
  780 +
  781 + // Compose exception message
  782 + std::string msg = "Incompatible array dimensions (";
  783 + for (int i = 0; i < m_shape.size(); ++i) {
  784 + if (i > 0) {
  785 + msg += ", ";
  786 + }
  787 + msg += gammalib::str(m_shape[i]);
  788 + }
  789 + msg += ") <=> (";
  790 + for (int i = 0; i < array.m_shape.size(); ++i) {
  791 + if (i > 0) {
  792 + msg += ", ";
  793 + }
  794 + msg += gammalib::str(array.m_shape[i]);
  795 + }
  796 + msg += ").";
  797 +
  798 + // Throw exception
  799 + throw GException::invalid_argument(method, msg);
  800 +
  801 + } // endif: array shapes differed
  802 +
  803 + // Return
  804 + return;
  805 +}
  806 +
  807 +
  808 +/*==========================================================================
  809 + = =
  810 + = Friends =
  811 + = =
  812 + ==========================================================================*/
  813 +
  814 +/***********************************************************************//**
  815 + * @brief Computes minimum array element
  816 + *
  817 + * @param[in] array Array.
  818 + * @return Minimum array element.
  819 + ***************************************************************************/
  820 +double min(const GNdarray& array)
  821 +{
  822 + // Initialises result
  823 + double result = 0.0;
  824 +
  825 + // Continue only if we have elements
  826 + if (array.m_data.size() > 0) {
  827 +
  828 + // Search for minimum
  829 + result = array.m_data[0];
  830 + for (int i = 1; i < array.m_data.size(); ++i) {
  831 + if (array.m_data[i] < result) {
  832 + result = array.m_data[i];
  833 + }
  834 + }
  835 +
  836 + } // endif: there were elements
  837 +
  838 + // Returns minimum
  839 + return result;
  840 +}
  841 +
  842 +
  843 +/***********************************************************************//**
  844 + * @brief Computes maximum array element
  845 + *
  846 + * @param[in] array Array.
  847 + * @return Maximum array element.
  848 + ***************************************************************************/
  849 +double max(const GNdarray& array)
  850 +{
  851 + // Initialises result
  852 + double result = 0.0;
  853 +
  854 + // Continue only if we have elements
  855 + if (array.m_data.size() > 0) {
  856 +
  857 + // Search for maximum
  858 + result = array.m_data[0];
  859 + for (int i = 1; i < array.m_data.size(); ++i) {
  860 + if (array.m_data[i] > result) {
  861 + result = array.m_data[i];
  862 + }
  863 + }
  864 +
  865 + } // endif: there were elements
  866 +
  867 + // Returns maximum
  868 + return result;
  869 +}
  870 +
  871 +
  872 +/***********************************************************************//**
  873 + * @brief Computes array sum
  874 + *
  875 + * @param[in] array Array.
  876 + * @return Sum of array elements.
  877 + ***************************************************************************/
  878 +double sum(const GNdarray& array)
  879 +{
  880 + // Compute sum
  881 + double result = 0.0;
  882 + for (int i = 0; i < array.m_data.size(); ++i) {
  883 + result += array.m_data[i];
  884 + }
  885 +
  886 + // Returns sum
  887 + return result;
  888 +}
  889 +
  890 +
  891 +/***********************************************************************//**
  892 + * @brief Computes arccos of array elements
  893 + *
  894 + * @param[in] array Array.
  895 + * @return Array containing the arccos of every element.
  896 + ***************************************************************************/
  897 +GNdarray acos(const GNdarray& array)
  898 +{
  899 + // Initialise result array
  900 + GNdarray result = array;
  901 +
  902 + // Evaluate each array element
  903 + for (int i = 0; i < array.m_data.size(); ++i) {
  904 + result.m_data[i] = std::acos(array.m_data[i]);
  905 + }
  906 +
  907 + // Return array
  908 + return result;
  909 +}
  910 +
  911 +
  912 +/***********************************************************************//**
  913 + * @brief Computes acosh of array elements
  914 + *
  915 + * @param[in] array Array.
  916 + * @return Array containing the acosh of every element.
  917 + ***************************************************************************/
  918 +GNdarray acosh(const GNdarray& array)
  919 +{
  920 + // Initialise result array
  921 + GNdarray result = array;
  922 +
  923 + // Evaluate each array element
  924 + for (int i = 0; i < array.m_data.size(); ++i) {
  925 + result.m_data[i] = acosh(array.m_data[i]);
  926 + }
  927 +
  928 + // Return array
  929 + return result;
  930 +}
  931 +
  932 +
  933 +/***********************************************************************//**
  934 + * @brief Computes arcsin of array elements
  935 + *
  936 + * @param[in] array Array.
  937 + * @return Array containing the arcsin of every element.
  938 + ***************************************************************************/
  939 +GNdarray asin(const GNdarray& array)
  940 +{
  941 + // Initialise result array
  942 + GNdarray result = array;
  943 +
  944 + // Evaluate each array element
  945 + for (int i = 0; i < array.m_data.size(); ++i) {
  946 + result.m_data[i] = std::asin(array.m_data[i]);
  947 + }
  948 +
  949 + // Return array
  950 + return result;
  951 +}
  952 +
  953 +
  954 +/***********************************************************************//**
  955 + * @brief Computes asinh of array elements
  956 + *
  957 + * @param[in] array Array.
  958 + * @return Array containing the asinh of every element.
  959 + ***************************************************************************/
  960 +GNdarray asinh(const GNdarray& array)
  961 +{
  962 + // Initialise result array
  963 + GNdarray result = array;
  964 +
  965 + // Evaluate each array element
  966 + for (int i = 0; i < array.m_data.size(); ++i) {
  967 + result.m_data[i] = asinh(array.m_data[i]);
  968 + }
  969 +
  970 + // Return array
  971 + return result;
  972 +}
  973 +
  974 +
  975 +/***********************************************************************//**
  976 + * @brief Computes arctan of array elements
  977 + *
  978 + * @param[in] array Array.
  979 + * @return Array containing the arctan of every element.
  980 + ***************************************************************************/
  981 +GNdarray atan(const GNdarray& array)
  982 +{
  983 + // Initialise result array
  984 + GNdarray result = array;
  985 +
  986 + // Evaluate each array element
  987 + for (int i = 0; i < array.m_data.size(); ++i) {
  988 + result.m_data[i] = std::atan(array.m_data[i]);
  989 + }
  990 +
  991 + // Return array
  992 + return result;
  993 +}
  994 +
  995 +
  996 +/***********************************************************************//**
  997 + * @brief Computes atanh of array elements
  998 + *
  999 + * @param[in] array Array.
  1000 + * @return Array containing the atanh of every element.
  1001 + ***************************************************************************/
  1002 +GNdarray atanh(const GNdarray& array)
  1003 +{
  1004 + // Initialise result array
  1005 + GNdarray result = array;
  1006 +
  1007 + // Evaluate each array element
  1008 + for (int i = 0; i < array.m_data.size(); ++i) {
  1009 + result.m_data[i] = atanh(array.m_data[i]);
  1010 + }
  1011 +
  1012 + // Return array
  1013 + return result;
  1014 +}
  1015 +
  1016 +
  1017 +/***********************************************************************//**
  1018 + * @brief Computes cosine of array elements
  1019 + *
  1020 + * @param[in] array Array.
  1021 + * @return Array containing the cosine of every element.
  1022 + ***************************************************************************/
  1023 +GNdarray cos(const GNdarray& array)
  1024 +{
  1025 + // Initialise result array
  1026 + GNdarray result = array;
  1027 +
  1028 + // Evaluate each array element
  1029 + for (int i = 0; i < array.m_data.size(); ++i) {
  1030 + result.m_data[i] = std::cos(array.m_data[i]);
  1031 + }
  1032 +
  1033 + // Return array
  1034 + return result;
  1035 +}
  1036 +
  1037 +
  1038 +/***********************************************************************//**
  1039 + * @brief Computes cosh of array elements
  1040 + *
  1041 + * @param[in] array Array.
  1042 + * @return Array containing the cosh of every element.
  1043 + ***************************************************************************/
  1044 +GNdarray cosh(const GNdarray& array)
  1045 +{
  1046 + // Initialise result array
  1047 + GNdarray result = array;
  1048 +
  1049 + // Evaluate each array element
  1050 + for (int i = 0; i < array.m_data.size(); ++i) {
  1051 + result.m_data[i] = std::cosh(array.m_data[i]);
  1052 + }
  1053 +
  1054 + // Return array
  1055 + return result;
  1056 +}
  1057 +
  1058 +
  1059 +/***********************************************************************//**
  1060 + * @brief Computes exponential of array elements
  1061 + *
  1062 + * @param[in] array Array.
  1063 + * @return Array containing the exponential of every element.
  1064 + ***************************************************************************/
  1065 +GNdarray exp(const GNdarray& array)
  1066 +{
  1067 + // Initialise result array
  1068 + GNdarray result = array;
  1069 +
  1070 + // Evaluate each array element
  1071 + for (int i = 0; i < array.m_data.size(); ++i) {
  1072 + result.m_data[i] = std::exp(array.m_data[i]);
  1073 + }
  1074 +
  1075 + // Return array
  1076 + return result;
  1077 +}
  1078 +
  1079 +
  1080 +/***********************************************************************//**
  1081 + * @brief Computes absolute of array elements
  1082 + *
  1083 + * @param[in] array Array.
  1084 + * @return Array containing the absolute of every element.
  1085 + ***************************************************************************/
  1086 +GNdarray abs(const GNdarray& array)
  1087 +{
  1088 + // Initialise result array
  1089 + GNdarray result = array;
  1090 +
  1091 + // Evaluate each array element
  1092 + for (int i = 0; i < array.m_data.size(); ++i) {
  1093 + result.m_data[i] = std::abs(array.m_data[i]);
  1094 + }
  1095 +
  1096 + // Return array
  1097 + return result;
  1098 +}
  1099 +
  1100 +
  1101 +/***********************************************************************//**
  1102 + * @brief Computes natural logarithm of array elements
  1103 + *
  1104 + * @param[in] array Array.
  1105 + * @return Array containing the natural logarithm of every element.
  1106 + ***************************************************************************/
  1107 +GNdarray log(const GNdarray& array)
  1108 +{
  1109 + // Initialise result array
  1110 + GNdarray result = array;
  1111 +
  1112 + // Evaluate each array element
  1113 + for (int i = 0; i < array.m_data.size(); ++i) {
  1114 + result.m_data[i] = std::log(array.m_data[i]);
  1115 + }
  1116 +
  1117 + // Return array
  1118 + return result;
  1119 +}
  1120 +
  1121 +
  1122 +/***********************************************************************//**
  1123 + * @brief Computes base10 logarithm of array elements
  1124 + *
  1125 + * @param[in] array Array.
  1126 + * @return Array containing the base10 logarithm of every element.
  1127 + ***************************************************************************/
  1128 +GNdarray log10(const GNdarray& array)
  1129 +{
  1130 + // Initialise result array
  1131 + GNdarray result = array;
  1132 +
  1133 + // Evaluate each array element
  1134 + for (int i = 0; i < array.m_data.size(); ++i) {
  1135 + result.m_data[i] = std::log10(array.m_data[i]);
  1136 + }
  1137 +
  1138 + // Return array
  1139 + return result;
  1140 +}
  1141 +
  1142 +
  1143 +/***********************************************************************//**
  1144 + * @brief Computes sine of array elements
  1145 + *
  1146 + * @param[in] array Array.
  1147 + * @return Array containing the sine of every element.
  1148 + ***************************************************************************/
  1149 +GNdarray sin(const GNdarray& array)
  1150 +{
  1151 + // Initialise result array
  1152 + GNdarray result = array;
  1153 +
  1154 + // Evaluate each array element
  1155 + for (int i = 0; i < array.m_data.size(); ++i) {
  1156 + result.m_data[i] = std::sin(array.m_data[i]);
  1157 + }
  1158 +
  1159 + // Return array
  1160 + return result;
  1161 +}
  1162 +
  1163 +
  1164 +/***********************************************************************//**
  1165 + * @brief Computes sinh of array elements
  1166 + *
  1167 + * @param[in] array Array.
  1168 + * @return Array containing the sinh of every element.
  1169 + ***************************************************************************/
  1170 +GNdarray sinh(const GNdarray& array)
  1171 +{
  1172 + // Initialise result array
  1173 + GNdarray result = array;
  1174 +
  1175 + // Evaluate each array element
  1176 + for (int i = 0; i < array.m_data.size(); ++i) {
  1177 + result.m_data[i] = std::sinh(array.m_data[i]);
  1178 + }
  1179 +
  1180 + // Return array
  1181 + return result;
  1182 +}
  1183 +
  1184 +
  1185 +/***********************************************************************//**
  1186 + * @brief Computes square root of array elements
  1187 + *
  1188 + * @param[in] array Array.
  1189 + * @return Array containing the square root of every element.
  1190 + ***************************************************************************/
  1191 +GNdarray sqrt(const GNdarray& array)
  1192 +{
  1193 + // Initialise result array
  1194 + GNdarray result = array;
  1195 +
  1196 + // Evaluate each array element
  1197 + for (int i = 0; i < array.m_data.size(); ++i) {
  1198 + result.m_data[i] = std::sqrt(array.m_data[i]);
  1199 + }
  1200 +
  1201 + // Return array
  1202 + return result;
  1203 +}
  1204 +
  1205 +
  1206 +/***********************************************************************//**
  1207 + * @brief Computes tangens of array elements
  1208 + *
  1209 + * @param[in] array Array.
  1210 + * @return Array containing the tangens of every element.
  1211 + ***************************************************************************/
  1212 +GNdarray tan(const GNdarray& array)
  1213 +{
  1214 + // Initialise result array
  1215 + GNdarray result = array;
  1216 +
  1217 + // Evaluate each array element
  1218 + for (int i = 0; i < array.m_data.size(); ++i) {
  1219 + result.m_data[i] = std::tan(array.m_data[i]);
  1220 + }
  1221 +
  1222 + // Return array
  1223 + return result;
  1224 +}
  1225 +
  1226 +
  1227 +/***********************************************************************//**
  1228 + * @brief Computes tanh of array elements
  1229 + *
  1230 + * @param[in] array Array.
  1231 + * @return Array containing the tanh of every element.
  1232 + ***************************************************************************/
  1233 +GNdarray tanh(const GNdarray& array)
  1234 +{
  1235 + // Initialise result array
  1236 + GNdarray result = array;
  1237 +
  1238 + // Evaluate each array element
  1239 + for (int i = 0; i < array.m_data.size(); ++i) {
  1240 + result.m_data[i] = std::tanh(array.m_data[i]);
  1241 + }
  1242 +
  1243 + // Return array
  1244 + return result;
  1245 +}
  1246 +
  1247 +
  1248 +/***********************************************************************//**
  1249 + * @brief Computes tanh of array elements
  1250 + *
  1251 + * @param[in] array Array.
  1252 + * @param[in] power Power.
  1253 + * @return Array containing the power of every element.
  1254 + ***************************************************************************/
  1255 +GNdarray pow(const GNdarray& array, const double& power)
  1256 +{
  1257 + // Initialise result array
  1258 + GNdarray result = array;
  1259 +
  1260 + // Evaluate each array element
  1261 + for (int i = 0; i < array.m_data.size(); ++i) {
  1262 + result.m_data[i] = std::pow(array.m_data[i], power);
  1263 + }
  1264 +
  1265 + // Return array
  1266 + return result;
  1267 +}
... ...
src/numerics/Makefile.am
  View file @32eaa5e
... ... @@ -29,6 +29,7 @@ sources = GIntegral.cpp \
29 29 GDerivative.cpp \
30 30 GFunction.cpp \
31 31 GMath.cpp \
  32 + GNdarray.cpp \
32 33 GException_numerics.cpp
33 34  
34 35 # Build libtool library
... ...
test/test_GNumerics.cpp
  View file @32eaa5e
1 1 /***************************************************************************
2 2 * test_GNumerics.cpp - test numerics modules *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2014 by Jurgen Knodlseder *
  4 + * copyright (C) 2010-2016 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -50,10 +50,16 @@ void TestGNumerics::set(void)
50 50 m_sigma = 2.5;
51 51  
52 52 // Append tests
53   - append(static_cast<pfunction>(&TestGNumerics::test_integral),"Test GIntegral");
54   - append(static_cast<pfunction>(&TestGNumerics::test_romberg_integration),"Test Romberg integration");
55   - append(static_cast<pfunction>(&TestGNumerics::test_adaptive_simpson_integration),"Test adaptive Simpson integration");
56   - append(static_cast<pfunction>(&TestGNumerics::test_gauss_kronrod_integration),"Test Gauss-Kronrod integration");
  53 + append(static_cast<pfunction>(&TestGNumerics::test_ndarray),
  54 + "Test GNdarray");
  55 + append(static_cast<pfunction>(&TestGNumerics::test_integral),
  56 + "Test GIntegral");
  57 + append(static_cast<pfunction>(&TestGNumerics::test_romberg_integration),
  58 + "Test Romberg integration");
  59 + append(static_cast<pfunction>(&TestGNumerics::test_adaptive_simpson_integration),
  60 + "Test adaptive Simpson integration");
  61 + append(static_cast<pfunction>(&TestGNumerics::test_gauss_kronrod_integration),
  62 + "Test Gauss-Kronrod integration");
57 63  
58 64 // Return
59 65 return;
... ... @@ -73,6 +79,203 @@ TestGNumerics* TestGNumerics::clone(void) const
73 79  
74 80  
75 81 /***********************************************************************//**
  82 + * @brief Test n-dimensional array
  83 + ***************************************************************************/
  84 +void TestGNumerics::test_ndarray(void)
  85 +{
  86 + // Test void constructor
  87 + GNdarray nd0;
  88 + test_value(nd0.size(), 0);
  89 + test_value(nd0.dim(), 0);
  90 +
  91 + // Test 1-dimensional array
  92 + GNdarray nd1(3);
  93 + test_value(nd1.size(), 3);
  94 + test_value(nd1.dim(), 1);
  95 + test_value(nd1.shape()[0], 3);
  96 + for (int i = 0; i < 3; ++i) {
  97 + nd1(i) = double(i);
  98 + test_value(nd1.at(i), double(i));
  99 + }
  100 +
  101 + // Test 2-dimensional array
  102 + GNdarray nd2(3,2);
  103 + test_value(nd2.size(), 6);
  104 + test_value(nd2.dim(), 2);
  105 + test_value(nd2.shape()[0], 3);
  106 + test_value(nd2.shape()[1], 2);
  107 + for (int ix = 0; ix < 3; ++ix) {
  108 + for (int iy = 0; iy < 2; ++iy) {
  109 + double ref = double(ix)+double(iy);
  110 + nd2(ix,iy) = ref;
  111 + test_value(nd2.at(ix,iy), ref);
  112 + }
  113 + }
  114 + // Test 3-dimensional array
  115 + GNdarray nd3(2,2,2);
  116 + test_value(nd3.size(), 8);
  117 + test_value(nd3.dim(), 3);
  118 + test_value(nd3.shape()[0], 2);
  119 + test_value(nd3.shape()[1], 2);
  120 + test_value(nd3.shape()[2], 2);
  121 + for (int ix = 0; ix < 2; ++ix) {
  122 + for (int iy = 0; iy < 2; ++iy) {
  123 + for (int iz = 0; iz < 2; ++iz) {
  124 + double ref = double(ix)+double(iy)+double(iz);
  125 + nd3(ix,iy,iz) = ref;
  126 + test_value(nd3.at(ix,iy,iz), ref);
  127 + }
  128 + }
  129 + }
  130 +
  131 + // Test n-dimensional array
  132 + std::vector<int> n(5,2);
  133 + GNdarray ndn(n);
  134 + test_value(ndn.size(), 32);
  135 + test_value(ndn.dim(), 5);
  136 + for (int i = 0; i < 5; ++i) {
  137 + test_value(ndn.shape()[i], 2);
  138 + }
  139 + std::vector<int> itest(5,1);
  140 + ndn(itest) = 3.1415;
  141 + test_value(ndn.at(itest), 3.1415);
  142 +
  143 + // Setup arrays for operator tests
  144 + GNdarray nda(1,2);
  145 + nda(0,0) = 3.0;
  146 + nda(0,1) = 5.5;
  147 + GNdarray ndb = nda;
  148 + GNdarray ndc = nda + ndb;
  149 +
  150 + // Test operator=
  151 + test_value(ndb(0,0), 3.0);
  152 + test_value(ndb(0,1), 5.5);
  153 +
  154 + // Test operator== and operator!=
  155 + test_assert(nda == ndb, "Test operator==");
  156 + test_assert(nda != ndc, "Test operator!=");
  157 +
  158 + // Test array operator+=
  159 + nda += ndb;
  160 + test_value(nda(0,0), 6.0);
  161 + test_value(nda(0,1), 11.0);
  162 +
  163 + // Test array operator-=
  164 + nda -= ndb;
  165 + test_value(nda(0,0), 3.0);
  166 + test_value(nda(0,1), 5.5);
  167 +
  168 + // Test value operator+=
  169 + nda += 1.5;
  170 + test_value(nda(0,0), 4.5);
  171 + test_value(nda(0,1), 7.0);
  172 +
  173 + // Test value operator-=
  174 + nda -= 1.5;
  175 + test_value(nda(0,0), 3.0);
  176 + test_value(nda(0,1), 5.5);
  177 +
  178 + // Test value operator*=
  179 + nda *= 3.0;
  180 + test_value(nda(0,0), 9.0);
  181 + test_value(nda(0,1), 16.5);
  182 +
  183 + // Test value operator/=
  184 + nda /= 3.0;
  185 + test_value(nda(0,0), 3.0);
  186 + test_value(nda(0,1), 5.5);
  187 +
  188 + // Test operator-
  189 + GNdarray ndd = -nda;
  190 + test_value(ndd(0,0), -3.0);
  191 + test_value(ndd(0,1), -5.5);
  192 +
  193 + // Test operator+(GNdarray,GNdarray)
  194 + ndd = nda + ndb;
  195 + test_value(ndd(0,0), 6.0);
  196 + test_value(ndd(0,1), 11.0);
  197 +
  198 + // Test operator+(GNdarray,double) and operator+(double,GNdarray)
  199 + ndd = nda + 1.5;
  200 + test_value(ndd(0,0), 4.5);
  201 + test_value(ndd(0,1), 7.0);
  202 + ndd = 1.5 + nda;
  203 + test_value(ndd(0,0), 4.5);
  204 + test_value(ndd(0,1), 7.0);
  205 +
  206 + // Test operator-(GNdarray,GNdarray)
  207 + ndd = nda - ndb;
  208 + test_value(ndd(0,0), 0.0);
  209 + test_value(ndd(0,1), 0.0);
  210 +
  211 + // Test operator-(GNdarray,double) and operator-(double,GNdarray)
  212 + ndd = nda - 1.5;
  213 + test_value(ndd(0,0), 1.5);
  214 + test_value(ndd(0,1), 4.0);
  215 + ndd = 1.5 - nda;
  216 + test_value(ndd(0,0), -1.5);
  217 + test_value(ndd(0,1), -4.0);
  218 +
  219 + // Test operator*(GNdarray,double) and operator*(double,GNdarray)
  220 + ndd = nda * 3.0;
  221 + test_value(ndd(0,0), 9.0);
  222 + test_value(ndd(0,1), 16.5);
  223 + ndd = 3.0 * nda;
  224 + test_value(ndd(0,0), 9.0);
  225 + test_value(ndd(0,1), 16.5);
  226 +
  227 + // Test operator/(GNdarray,double)
  228 + ndd = nda / 3.0;
  229 + test_value(ndd(0,0), 1.0);
  230 + test_value(ndd(0,1), 5.5/3.0);
  231 +
  232 + // Test friend functions
  233 + test_value(min(nda), 3.0);
  234 + test_value(max(nda), 5.5);
  235 + test_value(sum(nda), 8.5);
  236 + test_value(acos(nda)(0,0), std::acos(3.0));
  237 + test_value(acos(nda)(0,1), std::acos(5.5));
  238 + test_value(acosh(nda)(0,0), acosh(3.0));
  239 + test_value(acosh(nda)(0,1), acosh(5.5));
  240 + test_value(asin(nda)(0,0), std::asin(3.0));
  241 + test_value(asin(nda)(0,1), std::asin(5.5));
  242 + test_value(asinh(nda)(0,0), asinh(3.0));
  243 + test_value(asinh(nda)(0,1), asinh(5.5));
  244 + test_value(atan(nda)(0,0), std::atan(3.0));
  245 + test_value(atan(nda)(0,1), std::atan(5.5));
  246 + test_value(atanh(nda)(0,0), atanh(3.0));
  247 + test_value(atanh(nda)(0,1), atanh(5.5));
  248 + test_value(cos(nda)(0,0), std::cos(3.0));
  249 + test_value(cos(nda)(0,1), std::cos(5.5));
  250 + test_value(cosh(nda)(0,0), std::cosh(3.0));
  251 + test_value(cosh(nda)(0,1), std::cosh(5.5));
  252 + test_value(exp(nda)(0,0), std::exp(3.0));
  253 + test_value(exp(nda)(0,1), std::exp(5.5));
  254 + test_value(abs(nda)(0,0), std::abs(3.0));
  255 + test_value(abs(nda)(0,1), std::abs(5.5));
  256 + test_value(log(nda)(0,0), std::log(3.0));
  257 + test_value(log(nda)(0,1), std::log(5.5));
  258 + test_value(log10(nda)(0,0), std::log10(3.0));
  259 + test_value(log10(nda)(0,1), std::log10(5.5));
  260 + test_value(sin(nda)(0,0), std::sin(3.0));
  261 + test_value(sin(nda)(0,1), std::sin(5.5));
  262 + test_value(sinh(nda)(0,0), std::sinh(3.0));
  263 + test_value(sinh(nda)(0,1), std::sinh(5.5));
  264 + test_value(sqrt(nda)(0,0), std::sqrt(3.0));
  265 + test_value(sqrt(nda)(0,1), std::sqrt(5.5));
  266 + test_value(tan(nda)(0,0), std::tan(3.0));
  267 + test_value(tan(nda)(0,1), std::tan(5.5));
  268 + test_value(tanh(nda)(0,0), std::tanh(3.0));
  269 + test_value(tanh(nda)(0,1), std::tanh(5.5));
  270 + test_value(pow(nda,3.0)(0,0), std::pow(3.0,3.0));
  271 + test_value(pow(nda,3.0)(0,1), std::pow(5.5,3.0));
  272 +
  273 + // Exit test
  274 + return;
  275 +}
  276 +
  277 +
  278 +/***********************************************************************//**
76 279 * @brief Test model parameter handling
77 280 ***************************************************************************/
78 281 void TestGNumerics::test_integral(void)
... ...
test/test_GNumerics.hpp
  View file @32eaa5e
1 1 /***************************************************************************
2 2 * test_GNumerics.hpp - test numerics modules *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -68,6 +68,7 @@ public:
68 68 virtual void set(void);
69 69 virtual TestGNumerics* clone(void) const;
70 70 virtual std::string classname(void) const { return "TestGNumerics"; }
  71 + void test_ndarray(void);
71 72 void test_integral(void);
72 73 void test_romberg_integration(void);
73 74 void test_adaptive_simpson_integration(void);
... ...