GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GNodeArray.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GNodeArray.cpp - Array of nodes class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2008-2021 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GNodeArray.cpp
23  * @brief Node array class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cmath>
32 #include "GTools.hpp"
33 #include "GException.hpp"
34 #include "GFilename.hpp"
35 #include "GNodeArray.hpp"
36 #include "GVector.hpp"
37 #include "GFits.hpp"
38 #include "GFitsTable.hpp"
39 #include "GFitsBinTable.hpp"
40 #include "GFitsTableCol.hpp"
41 #include "GFitsTableDoubleCol.hpp"
42 
43 /* __ Method name definitions ____________________________________________ */
44 #define G_AT "GNodeArray::at(int&)"
45 #define G_INSERT "GNodeArray::insert(int&, double&)"
46 #define G_REMOVE "GNodeArray::remove(int&)"
47 #define G_INTERPOLATE "GNodeArray::interpolate(double&,"\
48  " std::vector<double>&)"
49 #define G_SET_VALUE "GNodeArray::set_value(double&)"
50 
51 /* __ Macros _____________________________________________________________ */
52 
53 /* __ Coding definitions _________________________________________________ */
54 
55 /* __ Debug definitions __________________________________________________ */
56 
57 /* __ Constants __________________________________________________________ */
58 
59 
60 /*==========================================================================
61  = =
62  = Constructors/destructors =
63  = =
64  ==========================================================================*/
65 
66 /***********************************************************************//**
67  * @brief Void constructor
68  ***************************************************************************/
70 {
71  // Initialise class members for clean destruction
72  init_members();
73 
74  // Return
75  return;
76 }
77 
78 
79 /***********************************************************************//**
80  * @brief FITS file constructor
81  *
82  * @param[in] filename FITS file name.
83  *
84  * Constructs node array from a FITS file.
85  ***************************************************************************/
87 {
88  // Initialise members
89  init_members();
90 
91  // Load FITS file
92  load(filename);
93 
94  // Return
95  return;
96 }
97 
98 
99 /***********************************************************************//**
100  * @brief Element constructor
101  *
102  * @param[in] num Number of elements.
103  * @param[in] array Array of elements.
104  *
105  * Constructs a node array from an @p array of double precision values of
106  * length @p num.
107  ***************************************************************************/
108 GNodeArray::GNodeArray(const int& num, const double* array)
109 {
110  // Initialise class members for clean destruction
111  init_members();
112 
113  // Use elements to initialise nodes
114  nodes(num, array);
115 
116  // Return
117  return;
118 }
119 
120 
121 /***********************************************************************//**
122  * @brief GVector constructor
123  *
124  * @param[in] vector Vector.
125  *
126  * Constructs a node array from the elements in a @p vector, represented by
127  * a GVector object.
128  ***************************************************************************/
130 {
131  // Initialise class members for clean destruction
132  init_members();
133 
134  // Use vector to initialise nodes
135  nodes(vector);
136 
137  // Return
138  return;
139 }
140 
141 
142 /***********************************************************************//**
143  * @brief Vector constructor
144  *
145  * @param[in] vector Vector.
146  *
147  * Constructs a node array from the elements in a double precision @p vector.
148  ***************************************************************************/
149 GNodeArray::GNodeArray(const std::vector<double>& vector)
150 {
151  // Initialise class members for clean destruction
152  init_members();
153 
154  // Use vector to initialise nodes
155  nodes(vector);
156 
157  // Return
158  return;
159 }
160 
161 
162 /***********************************************************************//**
163  * @brief Copy constructor
164  *
165  * @param[in] array Node array.
166  ***************************************************************************/
168 {
169  // Initialise class members for clean destruction
170  init_members();
171 
172  // Copy members
173  copy_members(array);
174 
175  // Return
176  return;
177 }
178 
179 
180 /***********************************************************************//**
181  * @brief Destructor
182  ***************************************************************************/
184 {
185  // Free members
186  free_members();
187 
188  // Return
189  return;
190 }
191 
192 
193 /*==========================================================================
194  = =
195  = Operators =
196  = =
197  ==========================================================================*/
198 
199 /***********************************************************************//**
200  * @brief Assignment operator
201  *
202  * @param[in] array Node array.
203  * @return Node array.
204  ***************************************************************************/
206 {
207  // Execute only if object is not identical
208  if (this != &array) {
209 
210  // Free members
211  free_members();
212 
213  // Initialise private members
214  init_members();
215 
216  // Copy members
217  copy_members(array);
218 
219  } // endif: object was not identical
220 
221  // Return this object
222  return *this;
223 }
224 
225 
226 /*==========================================================================
227  = =
228  = Public methods =
229  = =
230  ==========================================================================*/
231 
232 /***********************************************************************//**
233  * @brief Clear node array
234  ***************************************************************************/
236 {
237  // Free class members
238  free_members();
239 
240  // Initialise members
241  init_members();
242 
243  // Return
244  return;
245 }
246 
247 
248 /***********************************************************************//**
249  * @brief Clone node array
250  *
251  * @return Pointer to deep copy of node array.
252  ***************************************************************************/
254 {
255  return new GNodeArray(*this);
256 }
257 
258 
259 /***********************************************************************//**
260  * @brief Node access operator
261  *
262  * @param[in] index Node index [0,...,size()-1].
263  * @return Node value.
264  *
265  * @exception GException::out_of_range
266  * Node index is out of range.
267  *
268  * Returns a reference to the node with the specified @p index. The @p index
269  * is checked on its validity.
270  ***************************************************************************/
271 double& GNodeArray::at(const int& index)
272 {
273  // Compile option: raise an exception if index is out of range
274  #if defined(G_RANGE_CHECK)
275  if (index < 0 || index >= size()) {
276  throw GException::out_of_range(G_AT, "Node index", index, size());
277  }
278  #endif
279 
280  // Signal that setup needs to be called
281  m_need_setup = false;
282 
283  // Return node
284  return m_node[index];
285 }
286 
287 
288 /***********************************************************************//**
289  * @brief Node access operator (const version)
290  *
291  * @param[in] index Node index [0,...,size()-1].
292  * @return Node value.
293  *
294  * @exception GException::out_of_range
295  * Node index is out of range.
296  *
297  * Returns a reference to the node with the specified @p index. The @p index
298  * is checked on its validity.
299  ***************************************************************************/
300 const double& GNodeArray::at(const int& index) const
301 {
302  // Compile option: raise an exception if index is out of range
303  #if defined(G_RANGE_CHECK)
304  if (index < 0 || index >= size()) {
305  throw GException::out_of_range(G_AT, "Node index", index, size());
306  }
307  #endif
308 
309  // Signal that setup needs to be called
310  m_need_setup = false;
311 
312  // Return node
313  return m_node[index];
314 }
315 
316 
317 /***********************************************************************//**
318  * @brief Set node array
319  *
320  * @param[in] num Number of nodes
321  * @param[in] array Node values \f$x_i\f$.
322  *
323  * Setup node array from an array of double precision values.
324  ***************************************************************************/
325 void GNodeArray::nodes(const int& num, const double* array)
326 {
327  // Free any existing memory
328  free_members();
329 
330  // Initialise members
331  init_members();
332 
333  // Set node values
334  for (int i = 0; i < num; ++i) {
335  m_node.push_back(array[i]);
336  }
337 
338  // Setup node distances and linear array handling
339  setup();
340 
341  // Return
342  return;
343 }
344 
345 
346 /***********************************************************************//**
347  * @brief Append one node to array.
348  *
349  * @param[in] node Node.
350  ***************************************************************************/
351 void GNodeArray::append(const double& node)
352 {
353  // Add node
354  m_node.push_back(node);
355 
356  // Setup node distances and linear array handling
357  setup();
358 
359  // Return
360  return;
361 }
362 
363 
364 /***********************************************************************//**
365  * @brief Insert one node into array.
366  *
367  * @param[in] index Node index [0,...,size()-1].
368  * @param[in] node Node.
369  *
370  * @exception GException::out_of_range
371  * Node index is out of range.
372  *
373  * Inserts a @p node into the node array before the node with the specified
374  * @p index.
375  ***************************************************************************/
376 void GNodeArray::insert(const int& index, const double& node)
377 {
378  // Compile option: raise exception if index is out of range
379  #if defined(G_RANGE_CHECK)
380  if (is_empty()) {
381  if (index > 0) {
382  throw GException::out_of_range(G_INSERT, "Node index", index, size());
383  }
384  }
385  else {
386  if (index < 0 || index >= size()) {
387  throw GException::out_of_range(G_INSERT, "Node index", index, size());
388  }
389  }
390  #endif
391 
392  // Inserts node
393  m_node.insert(m_node.begin()+index, node);
394 
395  // Setup node distances and linear array handling
396  setup();
397 
398  // Return
399  return;
400 }
401 
402 
403 /***********************************************************************//**
404  * @brief Remove one node into array.
405  *
406  * @param[in] index Node index [0,...,size()-1].
407  *
408  * @exception GException::out_of_range
409  * Node index is out of range.
410  *
411  * Remove node of specified @p index from node array.
412  ***************************************************************************/
413 void GNodeArray::remove(const int& index)
414 {
415  // Compile option: raise exception if index is out of range
416  #if defined(G_RANGE_CHECK)
417  if (index < 0 || index >= size()) {
418  throw GException::out_of_range(G_REMOVE, "Node index", index, size());
419  }
420  #endif
421 
422  // Remove node
423  m_node.erase(m_node.begin() + index);
424 
425  // Setup node distances and linear array handling
426  setup();
427 
428  // Return
429  return;
430 }
431 
432 
433 /***********************************************************************//**
434  * @brief Append node array
435  *
436  * @param[in] nodes Node array.
437  *
438  * Append a node array to the node array.
439  ***************************************************************************/
440 void GNodeArray::extend(const GNodeArray& nodes)
441 {
442  // Do nothing if node array is empty
443  if (!nodes.is_empty()) {
444 
445  // Get size. Note that we extract the size first to avoid an
446  // endless loop that arises when a container is appended to
447  // itself.
448  int num = nodes.size();
449 
450  // Reserve enough space
451  reserve(size() + num);
452 
453  // Append all nodes
454  for (int i = 0; i < num; ++i) {
455  m_node.push_back(nodes[i]);
456  }
457 
458  // Setup node distances and linear array handling
459  setup();
460 
461  } // endif: node array was not empty
462 
463  // Return
464  return;
465 }
466 
467 
468 /***********************************************************************//**
469  * @brief Set node array from vector
470  *
471  * @param[in] vector Vector from which node array will be built.
472  *
473  * Setup node array from a vector of values.
474  ***************************************************************************/
475 void GNodeArray::nodes(const GVector& vector)
476 {
477  // Free any existing memory
478  free_members();
479 
480  // Initialise members
481  init_members();
482 
483  // Set node values
484  for (int i = 0; i < vector.size(); ++i) {
485  m_node.push_back(vector[i]);
486  }
487 
488  // Setup node distances and linear array handling
489  setup();
490 
491  // Return
492  return;
493 }
494 
495 
496 /***********************************************************************//**
497  * @brief Set node array from vector
498  *
499  * @param[in] vector Vector from which node array will be built.
500  *
501  * Setup node array from a vector of double precision values.
502  ***************************************************************************/
503 void GNodeArray::nodes(const std::vector<double>& vector)
504 {
505  // Free any existing memory
506  free_members();
507 
508  // Initialise members
509  init_members();
510 
511  // Set node values
512  m_node = vector;
513 
514  // Setup node distances and linear array handling
515  setup();
516 
517  // Return
518  return;
519 }
520 
521 
522 /***********************************************************************//**
523  * @brief Interpolate value.
524  *
525  * @param[in] value Value \f$x\f$ at which interpolation should be done.
526  * @param[in] vector Vector \f$y_i\f$ that should be interpolated.
527  *
528  * @exception GException::invalid_value
529  * Not enough nodes for interpolation in node array.
530  * @exception GException::invalid_argument
531  * Size of node vector does not match the size of vector argument.
532  *
533  * This method performs a linear interpolation of values \f$y_i\f$. The
534  * corresponding values \f$x_i\f$ are stored in the node array.
535  ***************************************************************************/
536 double GNodeArray::interpolate(const double& value,
537  const std::vector<double>& vector) const
538 {
539  // Throw exception if there are not enough nodes
540  if (m_node.size() < 2) {
541  std::string msg = "There are "+gammalib::str(m_node.size())+" nodes "
542  "in node array, yet at least 2 nodes are required "
543  "for interpolation. Please specify more nodes "
544  "before calling the method.";
546  }
547 
548  // Throw exception if vectors have not the same size
549  if (m_node.size() != vector.size()) {
550  std::string msg = "Vector length "+gammalib::str(vector.size())+
551  " differs from size "+gammalib::str(m_node.size())+
552  " of node array. Please specify a vector of "
553  "length "+gammalib::str(m_node.size())+".";
555  }
556 
557  // Set interpolation value
558  set_value(value);
559 
560  // Interpolate
561  double y = vector[inx_left()] * wgt_left() +
562  vector[inx_right()] * wgt_right();
563 
564  // Return
565  return y;
566 }
567 
568 
569 /***********************************************************************//**
570  * @brief Set indices and weighting factors for interpolation
571  *
572  * @param[in] value Value for which the interpolation should be done.
573  *
574  * @exception GException::invalid_value
575  * No nodes are available for interpolation.
576  *
577  * Set the indices that bound the specified value and the corresponding
578  * weighting factors for linear interpolation. If the array has a linear
579  * form (i.e. the nodes are equidistant), an analytic formula is used to
580  * determine the boundary indices. If the nodes are not equidistant the
581  * boundary indices are searched by bisection. If there is only a single
582  * node, no interpolation is done and the index of this node is returned.
583  *
584  * Note that this method needs to be called after changing the node array
585  * to setup the corrected interpolation indices and weights.
586  ***************************************************************************/
587 void GNodeArray::set_value(const double& value) const
588 {
589  // Get number of nodes
590  int nodes = m_node.size();
591 
592  // Throw an exception if there are no nodes
593  if (nodes < 1) {
594  std::string msg = "Attempting to set interpolating value without "
595  "having any nodes. Interpolation can only be "
596  "done if nodes are available.";
598  }
599 
600  // Initialize computation flag
601  bool compute = true;
602 
603  // Update cache if required. If cache was updated, computation is
604  // enforced. Otherwise we check if the value has changed. The value
605  // check is only done when a last value has been recorded.
606  if (m_need_setup) {
607  setup();
608  }
609  else {
610  if (m_has_last_value && value == m_last_value) {
611  compute = false;
612  }
613  }
614 
615  // Continue only if computation is required
616  if (compute) {
617 
618  // Handle special case of a single node
619  if (nodes == 1) {
620  m_inx_left = 0;
621  m_inx_right = 0;
622  m_wgt_left = 1.0;
623  m_wgt_right = 0.0;
624  m_wgt_grad_left = 0.0;
625  m_wgt_grad_right = 0.0;
626  }
627 
628  // Handle all other cases
629  else {
630 
631  // If array is linear then get left index from analytic formula
632  if (m_is_linear) {
633 
634  // Set left index
635  m_inx_left = int(m_linear_slope * value + m_linear_offset);
636 
637  // Keep index in valid range
638  if (m_inx_left < 0) {
639  m_inx_left = 0;
640  }
641  else if (m_inx_left >= nodes-1) {
642  m_inx_left = nodes - 2;
643  }
644 
645  } // endif: array is linear
646 
647  // ... otherwise search the relevant indices by bisection
648  else {
649 
650  // Set left index if value is before first node
651  if (value < m_node[0]) {
652  m_inx_left = 0;
653  }
654 
655  // Set left index if value is after last node
656  else if (value > m_node[nodes-1]) {
657  m_inx_left = nodes - 2;
658  }
659 
660  // Set left index by bisection
661  else {
662  int low = 0;
663  int high = nodes - 1;
664  while ((high - low) > 1) {
665  int mid = (low+high) / 2;
666  if (m_node[mid] > value) {
667  high = mid;
668  }
669  else {
670  low = mid;
671  }
672  }
673  m_inx_left = low;
674  } // endelse: did bisection
675  }
676 
677  // Set right index
678  m_inx_right = m_inx_left + 1;
679 
680  // Set weighting factors and gradients
684  m_wgt_left = 1.0 - m_wgt_right;
685 
686  } // endelse: more than one node was present
687 
688  } // endif: computation was required
689 
690  // Store last value and signal availability
691  m_last_value = value;
692  m_has_last_value = true;
693 
694  // Return
695  return;
696 }
697 
698 
699 /***********************************************************************//**
700  * @brief Load nodes from FITS file
701  *
702  * @param[in] filename FITS filename.
703  *
704  * Loads the node array from a FITS file.
705  *
706  * If no extension name is provided, the node array is loaded from the
707  * "NODES" extension.
708  ***************************************************************************/
709 void GNodeArray::load(const GFilename& filename)
710 {
711  // Open FITS file
712  GFits fits(filename);
713 
714  // Get nodes table
715  const GFitsTable& table = *fits.table(filename.extname("NODES"));
716 
717  // Read nodes from table
718  read(table);
719 
720  // Close FITS file
721  fits.close();
722 
723  // Return
724  return;
725 }
726 
727 
728 /***********************************************************************//**
729  * @brief Save node array into FITS file
730  *
731  * @param[in] filename FITS filename.
732  * @param[in] clobber Overwrite an existing node array extension?
733  *
734  * Saves node array into a FITS file. If a file with the given @p filename
735  * does not yet exist it will be created, otherwise the method opens the
736  * existing file. Node arrays can only be appended to an existing file if the
737  * @p clobber flag is set to `true` (otherwise an exception is thrown).
738  *
739  * The method will append a binary FITS table containing the node array to
740  * the FITS file. The extension name can be specified as part of the
741  * @p filename. For example the @p filename
742  *
743  * myfile.fits[NODE ARRAY]
744  *
745  * will save the node array in the `NODE ARRAY` extension of the
746  * `myfile.fits` file. If the extension exists already in the file it will be
747  * replaced, otherwise a new extension will be created. If no extension name
748  * is provided, the method will use `NODES` as the default extension name
749  * for the node array.
750  ***************************************************************************/
751 void GNodeArray::save(const GFilename& filename, const bool& clobber) const
752 {
753  // Open or create FITS file (without extension name since the requested
754  // extension may not yet exist in the file)
755  GFits fits(filename.url(), true);
756 
757  // Write node array to FITS object
758  write(fits, filename.extname("NODES"));
759 
760  // Save to file
761  fits.save(clobber);
762 
763  // Return
764  return;
765 }
766 
767 
768 /***********************************************************************//**
769  * @brief Read nodes from FITS table
770  *
771  * @param[in] table FITS table.
772  *
773  * Reads the nodes from a FITS @p table.
774  ***************************************************************************/
775 void GNodeArray::read(const GFitsTable& table)
776 {
777  // Clear nodes
778  clear();
779 
780  // Extract number of nodes in FITS file
781  int num = table.nrows();
782 
783  // Continue only if there are nodes bins
784  if (num > 0) {
785 
786  // Get the column with the name "Value"
787  const GFitsTableCol* col_nodes = table["Value"];
788 
789  // Set nodes
790  for (int i = 0; i < num; ++i) {
791  append(col_nodes->real(i));
792  }
793 
794  // Setup node distances and linear array handling
795  setup();
796 
797  } // endif: there were nodes
798 
799  // Return
800  return;
801 }
802 
803 
804 /***********************************************************************//**
805  * @brief Write nodes into FITS object
806  *
807  * @param[in] fits FITS file.
808  * @param[in] extname Nodes extension name (defaults to "NODES")
809  *
810  * Writes nodes into FITS object.
811  ***************************************************************************/
812 void GNodeArray::write(GFits& fits, const std::string& extname) const
813 {
814  // Set number of nodes
815  int num = m_node.size();
816 
817  // Create node column
818  GFitsTableDoubleCol col_nodes("Value", num);
819 
820  // Fill node column
821  for (int i = 0; i < num; ++i) {
822  col_nodes(i) = m_node[i];
823  }
824 
825  // Create nodes table
826  GFitsBinTable table(num);
827  table.append(col_nodes);
828  table.extname(extname);
829 
830  // If the FITS object contains already an extension with the same
831  // name then remove now this extension
832  if (fits.contains(extname)) {
833  fits.remove(extname);
834  }
835 
836  // Append NODES table to FITS file
837  fits.append(table);
838 
839  // Return
840  return;
841 }
842 
843 
844 /***********************************************************************//**
845  * @brief Print nodes
846  *
847  * @param[in] chatter Chattiness (defaults to NORMAL).
848  * @return String containing nodes information.
849  ***************************************************************************/
850 std::string GNodeArray::print(const GChatter& chatter) const
851 {
852  // Initialise result string
853  std::string result;
854 
855  // Continue only if chatter is not silent
856  if (chatter != SILENT) {
857 
858  // Append header
859  result.append("=== GNodeArray ===");
860 
861  // Append array type
862  result.append("\n"+gammalib::parformat("Number of nodes in array"));
863  result.append(gammalib::str(size()));
864  if (m_is_linear) {
865  result.append("\n"+gammalib::parformat("Array type")+"linear");
866  result.append("\n"+gammalib::parformat("Linear slope"));
867  result.append(gammalib::str(m_linear_slope));
868  result.append("\n"+gammalib::parformat("Linear offset"));
869  result.append(gammalib::str(m_linear_offset));
870  }
871  else {
872  result.append("\n"+gammalib::parformat("Array type")+"nonlinear");
873  }
874 
875  // EXPLICIT: Append indices and weights
876  if (chatter >= EXPLICIT) {
877  result.append("\n"+gammalib::parformat("Indices and weights"));
878  result.append("("+gammalib::str(m_inx_left)+",");
879  result.append(gammalib::str(m_inx_right)+")=");
880  result.append("("+gammalib::str(m_wgt_left)+",");
881  result.append(gammalib::str(m_wgt_right)+")");
882  result.append(" grad=");
883  result.append("("+gammalib::str(m_wgt_grad_left)+",");
884  result.append(gammalib::str(m_wgt_grad_right)+")");
885  }
886 
887  // VERBOSE: Append nodes
888  if (chatter >= VERBOSE) {
889  for (int i = 0; i < size(); ++i) {
890  result.append("\n"+gammalib::parformat("Node "+gammalib::str(i)));
891  result.append(gammalib::str(m_node[i]));
892  if (i < m_step.size()) {
893  result.append(" (delta="+gammalib::str(m_step[i])+")");
894  }
895  }
896  }
897 
898  } // endif: chatter was not silent
899 
900  // Return result
901  return result;
902 }
903 
904 
905 /*==========================================================================
906  = =
907  = Private methods =
908  = =
909  ==========================================================================*/
910 
911 /***********************************************************************//**
912  * @brief Initialise class members
913  ***************************************************************************/
915 {
916  // Initialise members
917  m_node.clear();
918  m_step.clear();
919  m_is_linear = false;
920  m_has_last_value = false;
921  m_last_value = 0.0;
922  m_linear_slope = 0.0;
923  m_linear_offset = 0.0;
924  m_inx_left = 0;
925  m_inx_right = 0;
926  m_wgt_left = 0.0;
927  m_wgt_right = 0.0;
928  m_wgt_grad_left = 0.0;
929  m_wgt_grad_right = 0.0;
930  m_need_setup = false;
931 
932  // Return
933  return;
934 }
935 
936 
937 /***********************************************************************//**
938  * @brief Copy class members
939  *
940  * @param[in] array Node array.
941  ***************************************************************************/
943 {
944  // Copy number of bins
945  m_node = array.m_node;
946  m_step = array.m_step;
947  m_is_linear = array.m_is_linear;
949  m_last_value = array.m_last_value;
952  m_inx_left = array.m_inx_left;
953  m_inx_right = array.m_inx_right;
954  m_wgt_left = array.m_wgt_left;
955  m_wgt_right = array.m_wgt_right;
958  m_need_setup = array.m_need_setup;
959 
960  // Return
961  return;
962 }
963 
964 
965 /***********************************************************************//**
966  * @brief Delete class members
967  ***************************************************************************/
969 {
970  // Return
971  return;
972 }
973 
974 
975 /***********************************************************************//**
976  * @brief Compute distance array and linear slope/offset
977  *
978  * Precomputes values for fast interpolation. The precomputation requires
979  * at least 2 nodes to be present in the node array. If less than two
980  * nodes are present, the distance vector m_step will be empty and no
981  * computation is done.
982  ***************************************************************************/
983 void GNodeArray::setup(void) const
984 {
985  // Reset distance vector
986  m_step.clear();
987 
988  // Get number of nodes
989  int nodes = m_node.size();
990 
991  // Continue only if there are at least 2 nodes
992  if (nodes > 1) {
993 
994  // Setup distance array between subsequent nodes
995  for (int i = 0; i < nodes-1; ++i) {
996  m_step.push_back(m_node[i+1] - m_node[i]);
997  }
998 
999  // Evaluate linear slope and offset
1000  m_linear_slope = double(nodes-1) / (m_node[nodes-1] - m_node[0]);
1002 
1003  // Check if nodes form a linear array
1004  m_is_linear = true;
1005  for (int i = 0; i < nodes-1; ++i) {
1006  double eps = m_linear_slope * m_node[i] + m_linear_offset - double(i);
1007  if (std::abs(eps) > 1.0e-6) {
1008  m_is_linear = false;
1009  break;
1010  }
1011  }
1012 
1013  } // endif: there were at least two nodes
1014 
1015  // Signal that setup has been called
1016  m_need_setup = false;
1017 
1018  // Return
1019  return;
1020 }
GNodeArray & operator=(const GNodeArray &array)
Assignment operator.
Definition: GNodeArray.cpp:205
int size(void) const
Return number of nodes in node array.
Definition: GNodeArray.hpp:192
FITS table double column class interface definition.
virtual ~GNodeArray(void)
Destructor.
Definition: GNodeArray.cpp:183
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:482
bool contains(const int &extno) const
Check if HDU exists in FITS file.
Definition: GFits.hpp:282
Node array class.
Definition: GNodeArray.hpp:60
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1253
void free_members(void)
Delete class members.
Definition: GNodeArray.cpp:968
void load(const GFilename &filename)
Load nodes from FITS file.
Definition: GNodeArray.cpp:709
bool m_has_last_value
Last value is valid.
Definition: GNodeArray.hpp:123
double m_linear_offset
Offset for linear array.
Definition: GNodeArray.hpp:127
const double & wgt_left(void) const
Returns left node weight.
Definition: GNodeArray.hpp:264
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
void reserve(const int &num)
Reserves space for nodes in node array.
Definition: GNodeArray.hpp:220
std::string extname(const std::string &defaultname="") const
Return extension name.
Definition: GFilename.cpp:385
GFitsTableCol * append(const GFitsTableCol &column)
Append column to the table.
Definition: GFitsTable.hpp:147
Gammalib tools definition.
FITS file class.
Definition: GFits.hpp:63
FITS file class interface definition.
#define G_SET_VALUE
Definition: GNodeArray.cpp:49
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
Definition: GNodeArray.cpp:587
double m_wgt_grad_left
Weight gradient for left node.
Definition: GNodeArray.hpp:132
double m_linear_slope
Slope for linear array.
Definition: GNodeArray.hpp:126
FITS table column abstract base class definition.
std::vector< double > m_node
Array of nodes.
Definition: GNodeArray.hpp:118
bool is_empty(void) const
Signals if there are no nodes in node array.
Definition: GNodeArray.hpp:206
void save(const GFilename &filename, const bool &clobber=false) const
Save node array into FITS file.
Definition: GNodeArray.cpp:751
void write(GFits &fits, const std::string &extname="NODES") const
Write nodes into FITS object.
Definition: GNodeArray.cpp:812
void init_members(void)
Initialise class members.
Definition: GNodeArray.cpp:914
double m_wgt_right
Weight for right node for linear interpolation.
Definition: GNodeArray.hpp:131
double m_last_value
Last requested value.
Definition: GNodeArray.hpp:125
Node array class interface definition.
void remove(const int &index)
Remove one node into array.
Definition: GNodeArray.cpp:413
double m_wgt_left
Weight for left node for linear interpolation.
Definition: GNodeArray.hpp:130
void setup(void) const
Compute distance array and linear slope/offset.
Definition: GNodeArray.cpp:983
const double & wgt_right(void) const
Returns right node weight.
Definition: GNodeArray.hpp:279
void remove(const int &extno)
Remove HDU from FITS file.
Definition: GFits.cpp:862
Filename class.
Definition: GFilename.hpp:62
Abstract interface for FITS table column.
double & at(const int &index)
Node access operator.
Definition: GNodeArray.cpp:271
std::vector< double > m_step
Distance to next node.
Definition: GNodeArray.hpp:124
double interpolate(const double &value, const std::vector< double > &vector) const
Interpolate value.
Definition: GNodeArray.cpp:536
void extend(const GNodeArray &nodes)
Append node array.
Definition: GNodeArray.cpp:440
#define G_INSERT
Definition: GNodeArray.cpp:45
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
GChatter
Definition: GTypemaps.hpp:33
const int & inx_left(void) const
Returns left node index.
Definition: GNodeArray.hpp:235
bool m_need_setup
Call of setup is required.
Definition: GNodeArray.hpp:121
Vector class interface definition.
GNodeArray(void)
Void constructor.
Definition: GNodeArray.cpp:69
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
int m_inx_left
Index of left node for linear interpolation.
Definition: GNodeArray.hpp:128
int m_inx_right
Index of right node for linear interpolation.
Definition: GNodeArray.hpp:129
void nodes(const int &num, const double *array)
Set node array.
Definition: GNodeArray.cpp:325
const int & inx_right(void) const
Returns right node index.
Definition: GNodeArray.hpp:249
const int & nrows(void) const
Return number of rows in table.
Definition: GFitsTable.hpp:119
double m_wgt_grad_right
Weight gradient for right node.
Definition: GNodeArray.hpp:133
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
std::string print(const GChatter &chatter=NORMAL) const
Print nodes.
Definition: GNodeArray.cpp:850
#define G_INTERPOLATE
Definition: GNodeArray.cpp:47
virtual double real(const int &row, const int &inx=0) const =0
GNodeArray * clone(void) const
Clone node array.
Definition: GNodeArray.cpp:253
void insert(const int &index, const double &node)
Insert one node into array.
Definition: GNodeArray.cpp:376
FITS binary table class.
#define G_REMOVE
Definition: GNodeArray.cpp:46
Exception handler interface definition.
GFitsHDU * append(const GFitsHDU &hdu)
Append HDU to FITS file.
Definition: GFits.cpp:678
FITS binary table class definition.
void read(const GFitsTable &table)
Read nodes from FITS table.
Definition: GNodeArray.cpp:775
const int & size(void) const
Return size of vector.
Definition: GVector.hpp:178
Vector class.
Definition: GVector.hpp:46
bool m_is_linear
Nodes form a linear array.
Definition: GNodeArray.hpp:122
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
void close(void)
Close FITS file.
Definition: GFits.cpp:1342
#define G_AT
Definition: GNodeArray.cpp:44
FITS table double column.
Filename class interface definition.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489
FITS table abstract base class interface definition.
void copy_members(const GNodeArray &array)
Copy class members.
Definition: GNodeArray.cpp:942