Olga Sergijenko / gammalib

Commit ca7c6322b56e2a9ab4dbe43f84267289dedbd234

Authored by Jürgen Knödlseder
1 parent 6fe87dfe

Use GVector instead of GNdarray in GFunctions and GIntegrals classes

Inline Side-by-side
Showing 5 changed files with 66 additions and 66 deletions
include/GFunctions.hpp
  View file @ca7c632
... ... @@ -30,7 +30,7 @@
30 30 /* __ Includes ___________________________________________________________ */
31 31  
32 32 /* __ Forward declarations _______________________________________________ */
33   -class GNdarray;
  33 +class GVector;
34 34  
35 35  
36 36 /***********************************************************************//**
... ... @@ -60,8 +60,8 @@ public:
60 60 GFunctions& operator=(const GFunctions& functions);
61 61  
62 62 // Methods
63   - virtual const GNdarray& array(void) const = 0;
64   - virtual GNdarray eval(const double& x) = 0;
  63 + virtual int size(void) const = 0;
  64 + virtual GVector eval(const double& x) = 0;
65 65  
66 66 protected:
67 67 // Protected methods
... ...
include/GIntegrals.hpp
  View file @ca7c632
... ... @@ -31,7 +31,7 @@
31 31 #include <string>
32 32 #include <vector>
33 33 #include "GBase.hpp"
34   -#include "GNdarray.hpp"
  34 +#include "GVector.hpp"
35 35  
36 36 /* __ Forward declarations _______________________________________________ */
37 37 class GFunctions;
... ... @@ -75,15 +75,15 @@ public:
75 75 const std::string& message(void) const;
76 76 void kernels(GFunctions* kernels);
77 77 const GFunctions* kernels(void) const;
78   - GNdarray romberg(std::vector<double> bounds,
  78 + GVector romberg(std::vector<double> bounds,
79 79 const int& order = 5);
80   - GNdarray romberg(const double& a,
  80 + GVector romberg(const double& a,
81 81 const double& b,
82 82 const int& order = 5);
83   - GNdarray trapzd(const double& a,
  83 + GVector trapzd(const double& a,
84 84 const double& b,
85 85 const int& n,
86   - GNdarray result);
  86 + GVector result);
87 87 std::string print(const GChatter& chatter = NORMAL) const;
88 88  
89 89 protected:
... ... @@ -91,12 +91,12 @@ protected:
91 91 void init_members(void);
92 92 void copy_members(const GIntegrals& integral);
93 93 void free_members(void);
94   - double polint(const double* xa,
95   - const GNdarray* ya,
96   - const int& n,
97   - const int& index,
98   - const double& x,
99   - double* dy);
  94 + double polint(const double* xa,
  95 + const GVector* ya,
  96 + const int& n,
  97 + const int& index,
  98 + const double& x,
  99 + double* dy);
100 100  
101 101 // Protected data area
102 102 GFunctions* m_kernels; //!< Pointer to function kernels
... ... @@ -111,8 +111,8 @@ protected:
111 111 mutable bool m_isvalid; //!< Integration result valid (true=yes)
112 112 mutable bool m_has_abserr; //!< Has absolute integration errors
113 113 mutable bool m_has_relerr; //!< Has relative integration errors
114   - mutable GNdarray m_abserr; //!< Absolute integration errors
115   - mutable GNdarray m_relerr; //!< Absolute integration errors
  114 + mutable GVector m_abserr; //!< Absolute integration errors
  115 + mutable GVector m_relerr; //!< Absolute integration errors
116 116 mutable std::string m_message; //!< Status message (if result is invalid)
117 117 };
118 118  
... ...
pyext/GFunctions.i
  View file @ca7c632
... ... @@ -26,7 +26,7 @@
26 26 %{
27 27 /* Put headers and other declarations here that are needed for compilation */
28 28 #include "GFunctions.hpp"
29   -#include "GNdarray.hpp"
  29 +#include "GVector.hpp"
30 30 %}
31 31  
32 32  
... ... @@ -43,6 +43,6 @@ public:
43 43 virtual ~GFunction(void);
44 44  
45 45 // Methods
46   - virtual const GNdarray& array(void) const = 0;
47   - virtual GNdarray eval(const double& x) = 0;
  46 + virtual int size(void) const = 0;
  47 + virtual GVector eval(const double& x) = 0;
48 48 };
... ...
pyext/GIntegrals.i
  View file @ca7c632
... ... @@ -64,15 +64,15 @@ public:
64 64 const std::string& message(void) const;
65 65 void kernels(GFunctions* kernels);
66 66 const GFunctions* kernels(void) const;
67   - GNdarray romberg(std::vector<double> bounds,
  67 + GVector romberg(std::vector<double> bounds,
68 68 const int& order = 5);
69   - GNdarray romberg(const double& a,
  69 + GVector romberg(const double& a,
70 70 const double& b,
71 71 const int& order = 5);
72   - GNdarray trapzd(const double& a,
  72 + GVector trapzd(const double& a,
73 73 const double& b,
74 74 const int& n,
75   - GNdarray result);
  75 + GVector result);
76 76 };
77 77  
78 78  
... ...
src/numerics/GIntegrals.cpp
  View file @ca7c632
... ... @@ -37,9 +37,9 @@
37 37 #define G_ROMBERG1 "GIntegrals::romberg(std::vector<double>, int&)"
38 38 #define G_ROMBERG2 "GIntegrals::romberg(double&, double&, int&)"
39 39 #define G_GAUSS_KRONROD "GIntegrals::gauss_kronrod(double&, double&)"
40   -#define G_TRAPZD "GIntegrals::trapzd(double&, double&, int&, GNdarray)"
41   -#define G_POLINT "GIntegrals::polint(double*, GNdarray*, int, int, double, "\
42   - "GNdarray*)"
  40 +#define G_TRAPZD "GIntegrals::trapzd(double&, double&, int&, GVector)"
  41 +#define G_POLINT "GIntegrals::polint(double*, GVector*, int, int, double, "\
  42 + "GVector*)"
43 43  
44 44 /* __ Macros _____________________________________________________________ */
45 45  
... ... @@ -192,7 +192,7 @@ GIntegrals* GIntegrals::clone(void) const
192 192 *
193 193 * @param[in] bounds Integration boundaries.
194 194 * @param[in] order Integration order (default: 5)
195   - * @return Array of integration results.
  195 + * @return Vector of integration results.
196 196 *
197 197 * Returns the integral of the integrand, computed over a number of
198 198 * intervals [a0,a1], [a1,a2], ... that are given as an unordered vector
... ... @@ -207,7 +207,7 @@ GIntegrals* GIntegrals::clone(void) const
207 207 * The number of iterations is limited by m_max_iter. m_eps specifies the
208 208 * requested fractional accuracy. By default it is set to 1e-6.
209 209 ***************************************************************************/
210   -GNdarray GIntegrals::romberg(std::vector<double> bounds, const int& order)
  210 +GVector GIntegrals::romberg(std::vector<double> bounds, const int& order)
211 211 {
212 212 // Throw an exception if the instance has no kernels
213 213 if (m_kernels == NULL) {
... ... @@ -220,7 +220,7 @@ GNdarray GIntegrals::romberg(std::vector&lt;double&gt; bounds, const int&amp; order)
220 220 std::sort(bounds.begin(), bounds.end());
221 221  
222 222 // Initialise result
223   - GNdarray result(m_kernels->array().shape());
  223 + GVector result(m_kernels->size());
224 224  
225 225 // Initialise integration status information
226 226 int calls = 0;
... ... @@ -245,7 +245,7 @@ GNdarray GIntegrals::romberg(std::vector&lt;double&gt; bounds, const int&amp; order)
245 245 * @param[in] a Left integration boundary.
246 246 * @param[in] b Right integration boundary.
247 247 * @param[in] order Integration order (default: 5)
248   - * @return Array of integration results.
  248 + * @return Vector of integration results.
249 249 *
250 250 * @exception GException::invalid_value
251 251 * Function kernels not set.
... ... @@ -262,7 +262,7 @@ GNdarray GIntegrals::romberg(std::vector&lt;double&gt; bounds, const int&amp; order)
262 262 * The number of iterations is limited by m_max_iter. m_eps specifies the
263 263 * requested fractional accuracy. By default it is set to 1e-6.
264 264 ***************************************************************************/
265   -GNdarray GIntegrals::romberg(const double& a, const double& b, const int& order)
  265 +GVector GIntegrals::romberg(const double& a, const double& b, const int& order)
266 266 {
267 267 // Throw an exception if the instance has no kernels
268 268 if (m_kernels == NULL) {
... ... @@ -272,7 +272,7 @@ GNdarray GIntegrals::romberg(const double&amp; a, const double&amp; b, const int&amp; order)
272 272 }
273 273  
274 274 // Initialise result
275   - GNdarray result(m_kernels->array().shape());
  275 + GVector result(m_kernels->size());
276 276  
277 277 // Initialise integration status information
278 278 m_isvalid = true;
... ... @@ -284,8 +284,8 @@ GNdarray GIntegrals::romberg(const double&amp; a, const double&amp; b, const int&amp; order)
284 284 if (b > a) {
285 285  
286 286 // Initialise variables
287   - bool converged = false;
288   - GNdarray dss(m_kernels->array().shape());
  287 + bool converged = false;
  288 + GVector dss(m_kernels->size());
289 289  
290 290 // Determine (maximum) number of iterations
291 291 int max_iter = (m_fix_iter > 0) ? m_fix_iter : m_max_iter;
... ... @@ -303,12 +303,12 @@ GNdarray GIntegrals::romberg(const double&amp; a, const double&amp; b, const int&amp; order)
303 303 }
304 304  
305 305 // Allocate temporal storage
306   - double* h = new double[max_iter+2];
307   - GNdarray* s = new GNdarray[max_iter+2];
  306 + double* h = new double[max_iter+2];
  307 + GVector* s = new GVector[max_iter+2];
308 308  
309 309 // Initialise step size and initial result
310 310 h[1] = 1.0;
311   - s[0] = GNdarray(m_kernels->array().shape());
  311 + s[0] = GVector(m_kernels->size());
312 312  
313 313 // Iterative loop
314 314 for (m_iter = 1; m_iter <= max_iter; ++m_iter) {
... ... @@ -335,9 +335,9 @@ GNdarray GIntegrals::romberg(const double&amp; a, const double&amp; b, const int&amp; order)
335 335 // Compute result using polynom interpolation
336 336 for (int i = 0; i < result.size(); ++i) {
337 337 double dy = 0;
338   - result(i) = polint(&h[m_iter-order], &s[m_iter-order],
  338 + result[i] = polint(&h[m_iter-order], &s[m_iter-order],
339 339 order, i, 0.0, &dy);
340   - dss(i) = dy;
  340 + dss[i] = dy;
341 341 }
342 342  
343 343 // If a fixed number of iterations has been requested and if
... ... @@ -354,7 +354,7 @@ GNdarray GIntegrals::romberg(const double&amp; a, const double&amp; b, const int&amp; order)
354 354 else {
355 355 converged = true;
356 356 for (int i = 0; i < result.size(); ++i) {
357   - if (std::abs(dss(i)) > m_eps * std::abs(result(i))) {
  357 + if (std::abs(dss[i]) > m_eps * std::abs(result[i])) {
358 358 converged = false;
359 359 break;
360 360 }
... ... @@ -371,13 +371,13 @@ GNdarray GIntegrals::romberg(const double&amp; a, const double&amp; b, const int&amp; order)
371 371 m_abserr = abs(dss);
372 372 m_relerr = m_abserr;
373 373 for (int i = 0; i < result.size(); ++i) {
374   - double abs_result = std::abs(result(i));
  374 + double abs_result = std::abs(result[i]);
375 375 if (abs_result > 0.0) {
376   - m_relerr(i) /= abs_result;
  376 + m_relerr[i] /= abs_result;
377 377 m_has_relerr = true;
378 378 }
379 379 else {
380   - m_relerr(i) = 0.0;
  380 + m_relerr[i] = 0.0;
381 381 }
382 382 }
383 383 break;
... ... @@ -424,8 +424,8 @@ GNdarray GIntegrals::romberg(const double&amp; a, const double&amp; b, const int&amp; order)
424 424 * @param[in] a Left integration boundary.
425 425 * @param[in] b Right integration boundary.
426 426 * @param[in] n Number of steps.
427   - * @param[in] result Result array from a previous trapezoidal integration step.
428   - * @return Array of integration results.
  427 + * @param[in] result Result vector from a previous trapezoidal integration step.
  428 + * @return Vector of integration results.
429 429 *
430 430 * @exception GException::invalid_value
431 431 * Function kernels not set.
... ... @@ -451,10 +451,10 @@ GNdarray GIntegrals::romberg(const double&amp; a, const double&amp; b, const int&amp; order)
451 451 * where \f${\tt result}_j\f$ is the integration result for function \f$j\f$
452 452 * from a previous call to the method with @p n = @p n - 1.
453 453 ***************************************************************************/
454   -GNdarray GIntegrals::trapzd(const double& a,
455   - const double& b,
456   - const int& n,
457   - GNdarray result)
  454 +GVector GIntegrals::trapzd(const double& a,
  455 + const double& b,
  456 + const int& n,
  457 + GVector result)
458 458 {
459 459 // Throw an exception if the instance has no kernel
460 460 if (m_kernels == NULL) {
... ... @@ -465,7 +465,7 @@ GNdarray GIntegrals::trapzd(const double&amp; a,
465 465  
466 466 // Handle case of identical boundaries
467 467 if (a == b) {
468   - result = GNdarray(result.shape());
  468 + result = GVector(result.size());
469 469 }
470 470  
471 471 // ... otherwise use trapeziodal rule
... ... @@ -475,9 +475,9 @@ GNdarray GIntegrals::trapzd(const double&amp; a,
475 475 if (n == 1) {
476 476  
477 477 // Evaluate integrand at boundaries
478   - GNdarray y_a = m_kernels->eval(a);
479   - GNdarray y_b = m_kernels->eval(b);
480   - m_calls += 2;
  478 + GVector y_a = m_kernels->eval(a);
  479 + GVector y_b = m_kernels->eval(b);
  480 + m_calls += 2;
481 481  
482 482 // Compute result
483 483 result = 0.5 * (b-a) * (y_a + y_b);
... ... @@ -528,8 +528,8 @@ GNdarray GIntegrals::trapzd(const double&amp; a,
528 528 }
529 529  
530 530 // Sum up values
531   - double x = a + 0.5 * del;
532   - GNdarray sum(result.shape());
  531 + double x = a + 0.5 * del;
  532 + GVector sum(result.size());
533 533 for (int j = 0; j < it; ++j, x += del) {
534 534  
535 535 // Evaluate and add integrand
... ... @@ -692,23 +692,23 @@ void GIntegrals::free_members(void)
692 692 * @brief Perform Polynomial interpolation
693 693 *
694 694 * @param[in] xa Pointer to array of X values.
695   - * @param[in] ya Pointer to Ndarray of Y values.
  695 + * @param[in] ya Pointer to GVector of Y values.
696 696 * @param[in] n Number of elements in arrays.
697   - * @param[in] index Ndarray index for with interpolations should be performed.
  697 + * @param[in] index GVector index for with interpolations should be performed.
698 698 * @param[in] x X value for which interpolations should be performed.
699   - * @param[out] dy Pointer to Ndarray of error estimates for interpolated values.
  699 + * @param[out] dy Pointer to array of error estimates for interpolated values.
700 700 * @return Interpolated value.
701 701 *
702 702 * Given arrays @p xa[1,..,n] and @p ya[1,..,n], and given a value @p x, this
703 703 * method returns a value y, and an error estimate @p dy. If P(x) is the
704 704 * polynomial of degree n-1, then the returned value y=P(x).
705 705 ***************************************************************************/
706   -double GIntegrals::polint(const double* xa,
707   - const GNdarray* ya,
708   - const int& n,
709   - const int& index,
710   - const double& x,
711   - double* dy)
  706 +double GIntegrals::polint(const double* xa,
  707 + const GVector* ya,
  708 + const int& n,
  709 + const int& index,
  710 + const double& x,
  711 + double* dy)
712 712 {
713 713 // Initialise result
714 714 double y = 0.0;
... ... @@ -728,13 +728,13 @@ double GIntegrals::polint(const double* xa,
728 728 ns = i;
729 729 dif = dift;
730 730 }
731   - double value = ya[i+1](index);
  731 + double value = ya[i+1][index];
732 732 c[i] = value;
733 733 d[i] = value;
734 734 }
735 735  
736 736 // Get initial approximation to y
737   - y = ya[ns+1](index);
  737 + y = ya[ns+1][index];
738 738 ns--;
739 739  
740 740 // Loop over each column of the tableau
... ...