GammaLib  1.7.0.dev
 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-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 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::not_enough_nodes
529  * Not enough nodes for interpolation in node array.
530  * @exception GException::vector_mismatch
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) {
542  }
543 
544  // Throw exception if vectors have not the same size
545  if (m_node.size() != vector.size()) {
547  vector.size());
548  }
549 
550  // Set interpolation value (circumvent const correctness)
551  set_value(value);
552 
553  // Interpolate
554  double y = vector[inx_left()] * wgt_left() +
555  vector[inx_right()] * wgt_right();
556 
557  // Return
558  return y;
559 }
560 
561 
562 /***********************************************************************//**
563  * @brief Set indices and weighting factors for interpolation
564  *
565  * @param[in] value Value for which the interpolation should be done.
566  *
567  * @exception GException::invalid_value
568  * No nodes are available for interpolation.
569  *
570  * Set the indices that bound the specified value and the corresponding
571  * weighting factors for linear interpolation. If the array has a linear
572  * form (i.e. the nodes are equidistant), an analytic formula is used to
573  * determine the boundary indices. If the nodes are not equidistant the
574  * boundary indices are searched by bisection. If there is only a single
575  * node, no interpolation is done and the index of this node is returned.
576  *
577  * Note that this method needs to be called after changing the node array
578  * to setup the corrected interpolation indices and weights.
579  ***************************************************************************/
580 void GNodeArray::set_value(const double& value) const
581 {
582  // Get number of nodes
583  int nodes = m_node.size();
584 
585  /*
586  // Throw an exception if less than 2 nodes are available
587  if (nodes < 2) {
588  throw GException::not_enough_nodes(G_SET_VALUE, nodes);
589  }
590  */
591  // Throw an exception if there are no nodes
592  if (nodes < 1) {
593  std::string msg = "Attempting to set interpolating value without "
594  "having any nodes. Interpolation can only be "
595  "done if nodes are available.";
597  }
598 
599  // Initialize computation flag
600  bool compute = true;
601 
602  // Update cache if required. If cache was updated, computation is
603  // enforced. Otherwise we check if the value has changed. The value
604  // check is only done when a last value has been recorded.
605  if (m_need_setup) {
606  setup();
607  }
608  else {
609  if (m_has_last_value && value == m_last_value) {
610  compute = false;
611  }
612  }
613 
614  // Continue only if computation is required
615  if (compute) {
616 
617  // Handle special case of a single node
618  if (nodes == 1) {
619  m_inx_left = 0;
620  m_inx_right = 0;
621  m_wgt_left = 1.0;
622  m_wgt_right = 0.0;
623  }
624 
625  // Handle all other cases
626  else {
627 
628  // If array is linear then get left index from analytic formula
629  if (m_is_linear) {
630 
631  // Set left index
632  m_inx_left = int(m_linear_slope * value + m_linear_offset);
633 
634  // Keep index in valid range
635  if (m_inx_left < 0) {
636  m_inx_left = 0;
637  }
638  else if (m_inx_left >= nodes-1) {
639  m_inx_left = nodes - 2;
640  }
641 
642  } // endif: array is linear
643 
644  // ... otherwise search the relevant indices by bisection
645  else {
646 
647  // Set left index if value is before first node
648  if (value < m_node[0]) {
649  m_inx_left = 0;
650  }
651 
652  // Set left index if value is after last node
653  else if (value > m_node[nodes-1]) {
654  m_inx_left = nodes - 2;
655  }
656 
657  // Set left index by bisection
658  else {
659  int low = 0;
660  int high = nodes - 1;
661  while ((high - low) > 1) {
662  int mid = (low+high) / 2;
663  if (m_node[mid] > value) {
664  high = mid;
665  }
666  else {
667  low = mid;
668  }
669  }
670  m_inx_left = low;
671  } // endelse: did bisection
672  }
673 
674  // Set right index
675  m_inx_right = m_inx_left + 1;
676 
677  // Set weighting factors
679  m_wgt_left = 1.0 - m_wgt_right;
680 
681  } // endelse: more than one node was present
682 
683  } // endif: computation was required
684 
685  // Store last value and signal availability
686  m_last_value = value;
687  m_has_last_value = true;
688 
689  // Return
690  return;
691 }
692 
693 
694 /***********************************************************************//**
695  * @brief Load nodes from FITS file
696  *
697  * @param[in] filename FITS filename.
698  *
699  * Loads the node array from a FITS file.
700  *
701  * If no extension name is provided, the node array is loaded from the
702  * "NODES" extension.
703  ***************************************************************************/
704 void GNodeArray::load(const GFilename& filename)
705 {
706  // Open FITS file
707  GFits fits(filename);
708 
709  // Get nodes table
710  const GFitsTable& table = *fits.table(filename.extname("NODES"));
711 
712  // Read nodes from table
713  read(table);
714 
715  // Close FITS file
716  fits.close();
717 
718  // Return
719  return;
720 }
721 
722 
723 /***********************************************************************//**
724  * @brief Save node array into FITS file
725  *
726  * @param[in] filename FITS filename.
727  * @param[in] clobber Overwrite an existing node array extension?
728  *
729  * Saves node array into a FITS file. If a file with the given @p filename
730  * does not yet exist it will be created, otherwise the method opens the
731  * existing file. Node arrays can only be appended to an existing file if the
732  * @p clobber flag is set to `true` (otherwise an exception is thrown).
733  *
734  * The method will append a binary FITS table containing the node array to
735  * the FITS file. The extension name can be specified as part of the
736  * @p filename. For example the @p filename
737  *
738  * myfile.fits[NODE ARRAY]
739  *
740  * will save the node array in the `NODE ARRAY` extension of the
741  * `myfile.fits` file. If the extension exists already in the file it will be
742  * replaced, otherwise a new extension will be created. If no extension name
743  * is provided, the method will use `NODES` as the default extension name
744  * for the node array.
745  ***************************************************************************/
746 void GNodeArray::save(const GFilename& filename, const bool& clobber) const
747 {
748  // Open or create FITS file (without extension name since the requested
749  // extension may not yet exist in the file)
750  GFits fits(filename.url(), true);
751 
752  // Write node array to FITS object
753  write(fits, filename.extname("NODES"));
754 
755  // Save to file
756  fits.save(clobber);
757 
758  // Return
759  return;
760 }
761 
762 
763 /***********************************************************************//**
764  * @brief Read nodes from FITS table
765  *
766  * @param[in] table FITS table.
767  *
768  * Reads the nodes from a FITS @p table.
769  ***************************************************************************/
770 void GNodeArray::read(const GFitsTable& table)
771 {
772  // Clear nodes
773  clear();
774 
775  // Extract number of nodes in FITS file
776  int num = table.nrows();
777 
778  // Continue only if there are nodes bins
779  if (num > 0) {
780 
781  // Get the column with the name "Value"
782  const GFitsTableCol* col_nodes = table["Value"];
783 
784  // Set nodes
785  for (int i = 0; i < num; ++i) {
786  append(col_nodes->real(i));
787  }
788 
789  // Setup node distances and linear array handling
790  setup();
791 
792  } // endif: there were nodes
793 
794  // Return
795  return;
796 }
797 
798 
799 /***********************************************************************//**
800  * @brief Write nodes into FITS object
801  *
802  * @param[in] fits FITS file.
803  * @param[in] extname Nodes extension name (defaults to "NODES")
804  *
805  * Writes nodes into FITS object.
806  ***************************************************************************/
807 void GNodeArray::write(GFits& fits, const std::string& extname) const
808 {
809  // Set number of nodes
810  int num = m_node.size();
811 
812  // Create node column
813  GFitsTableDoubleCol col_nodes("Value", num);
814 
815  // Fill node column
816  for (int i = 0; i < num; ++i) {
817  col_nodes(i) = m_node[i];
818  }
819 
820  // Create nodes table
821  GFitsBinTable table(num);
822  table.append(col_nodes);
823  table.extname(extname);
824 
825  // If the FITS object contains already an extension with the same
826  // name then remove now this extension
827  if (fits.contains(extname)) {
828  fits.remove(extname);
829  }
830 
831  // Append NODES table to FITS file
832  fits.append(table);
833 
834  // Return
835  return;
836 }
837 
838 
839 /***********************************************************************//**
840  * @brief Print nodes
841  *
842  * @param[in] chatter Chattiness (defaults to NORMAL).
843  * @return String containing nodes information.
844  ***************************************************************************/
845 std::string GNodeArray::print(const GChatter& chatter) const
846 {
847  // Initialise result string
848  std::string result;
849 
850  // Continue only if chatter is not silent
851  if (chatter != SILENT) {
852 
853  // Append header
854  result.append("=== GNodeArray ===");
855 
856  // Append array type
857  result.append("\n"+gammalib::parformat("Number of nodes in array"));
858  result.append(gammalib::str(size()));
859  if (m_is_linear) {
860  result.append("\n"+gammalib::parformat("Array type")+"linear");
861  result.append("\n"+gammalib::parformat("Linear slope"));
862  result.append(gammalib::str(m_linear_slope));
863  result.append("\n"+gammalib::parformat("Linear offset"));
864  result.append(gammalib::str(m_linear_offset));
865  }
866  else {
867  result.append("\n"+gammalib::parformat("Array type")+"nonlinear");
868  }
869 
870  // EXPLICIT: Append indices and weights
871  if (chatter >= EXPLICIT) {
872  result.append("\n"+gammalib::parformat("Indices and weights"));
873  result.append("("+gammalib::str(m_inx_left)+",");
874  result.append(gammalib::str(m_inx_right)+")=");
875  result.append("("+gammalib::str(m_wgt_left)+",");
876  result.append(gammalib::str(m_wgt_right)+")");
877  }
878 
879  // VERBOSE: Append nodes
880  if (chatter >= VERBOSE) {
881  for (int i = 0; i < size(); ++i) {
882  result.append("\n"+gammalib::parformat("Node "+gammalib::str(i)));
883  result.append(gammalib::str(m_node[i]));
884  if (i < m_step.size()) {
885  result.append(" (delta="+gammalib::str(m_step[i])+")");
886  }
887  }
888  }
889 
890  } // endif: chatter was not silent
891 
892  // Return result
893  return result;
894 }
895 
896 
897 /*==========================================================================
898  = =
899  = Private methods =
900  = =
901  ==========================================================================*/
902 
903 /***********************************************************************//**
904  * @brief Initialise class members
905  ***************************************************************************/
907 {
908  // Initialise members
909  m_node.clear();
910  m_step.clear();
911  m_is_linear = false;
912  m_has_last_value = false;
913  m_last_value = 0.0;
914  m_linear_slope = 0.0;
915  m_linear_offset = 0.0;
916  m_inx_left = 0;
917  m_inx_right = 0;
918  m_wgt_left = 0.0;
919  m_wgt_right = 0.0;
920  m_need_setup = false;
921 
922  // Return
923  return;
924 }
925 
926 
927 /***********************************************************************//**
928  * @brief Copy class members
929  *
930  * @param[in] array Node array.
931  ***************************************************************************/
933 {
934  // Copy number of bins
935  m_node = array.m_node;
936  m_step = array.m_step;
937  m_is_linear = array.m_is_linear;
939  m_last_value = array.m_last_value;
942  m_inx_left = array.m_inx_left;
943  m_inx_right = array.m_inx_right;
944  m_wgt_left = array.m_wgt_left;
945  m_wgt_right = array.m_wgt_right;
946  m_need_setup = array.m_need_setup;
947 
948  // Return
949  return;
950 }
951 
952 
953 /***********************************************************************//**
954  * @brief Delete class members
955  ***************************************************************************/
957 {
958  // Return
959  return;
960 }
961 
962 
963 /***********************************************************************//**
964  * @brief Compute distance array and linear slope/offset
965  *
966  * Precomputes values for fast interpolation. The precomputation requires
967  * at least 2 nodes to be present in the node array. If less than two
968  * nodes are present, the distance vector m_step will be empty and no
969  * computation is done.
970  ***************************************************************************/
971 void GNodeArray::setup(void) const
972 {
973  // Reset distance vector
974  m_step.clear();
975 
976  // Get number of nodes
977  int nodes = m_node.size();
978 
979  // Continue only if there are at least 2 nodes
980  if (nodes > 1) {
981 
982  // Setup distance array between subsequent nodes
983  for (int i = 0; i < nodes-1; ++i) {
984  m_step.push_back(m_node[i+1] - m_node[i]);
985  }
986 
987  // Evaluate linear slope and offset
988  m_linear_slope = double(nodes-1) / (m_node[nodes-1] - m_node[0]);
990 
991  // Check if nodes form a linear array
992  m_is_linear = true;
993  for (int i = 0; i < nodes-1; ++i) {
994  double eps = m_linear_slope * m_node[i] + m_linear_offset - double(i);
995  if (std::abs(eps) > 1.0e-6) {
996  m_is_linear = false;
997  break;
998  }
999  }
1000 
1001  } // endif: there were at least two nodes
1002 
1003  // Signal that setup has been called
1004  m_need_setup = false;
1005 
1006  // Return
1007  return;
1008 }
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:188
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:472
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:1163
void free_members(void)
Delete class members.
Definition: GNodeArray.cpp:956
void load(const GFilename &filename)
Load nodes from FITS file.
Definition: GNodeArray.cpp:704
bool m_has_last_value
Last value is valid.
Definition: GNodeArray.hpp:121
double m_linear_offset
Offset for linear array.
Definition: GNodeArray.hpp:125
const double & wgt_left(void) const
Returns left node weight.
Definition: GNodeArray.hpp:260
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:216
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:580
double m_linear_slope
Slope for linear array.
Definition: GNodeArray.hpp:124
FITS table column abstract base class definition.
std::vector< double > m_node
Array of nodes.
Definition: GNodeArray.hpp:116
bool is_empty(void) const
Signals if there are no nodes in node array.
Definition: GNodeArray.hpp:202
void save(const GFilename &filename, const bool &clobber=false) const
Save node array into FITS file.
Definition: GNodeArray.cpp:746
void write(GFits &fits, const std::string &extname="NODES") const
Write nodes into FITS object.
Definition: GNodeArray.cpp:807
void init_members(void)
Initialise class members.
Definition: GNodeArray.cpp:906
double m_wgt_right
Weight for right node for linear interpolation.
Definition: GNodeArray.hpp:129
double m_last_value
Last requested value.
Definition: GNodeArray.hpp:123
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:128
void setup(void) const
Compute distance array and linear slope/offset.
Definition: GNodeArray.cpp:971
const double & wgt_right(void) const
Returns right node weight.
Definition: GNodeArray.hpp:275
void remove(const int &extno)
Remove HDU from FITS file.
Definition: GFits.cpp:848
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:122
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:231
bool m_need_setup
Call of setup is required.
Definition: GNodeArray.hpp:119
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:126
int m_inx_right
Index of right node for linear interpolation.
Definition: GNodeArray.hpp:127
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:245
const int & nrows(void) const
Return number of rows in table.
Definition: GFitsTable.hpp:119
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:845
#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:665
FITS binary table class definition.
void read(const GFitsTable &table)
Read nodes from FITS table.
Definition: GNodeArray.cpp:770
const int & size(void) const
Return size of vector.
Definition: GVector.hpp:180
Vector class.
Definition: GVector.hpp:46
bool m_is_linear
Nodes form a linear array.
Definition: GNodeArray.hpp:120
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
void close(void)
Close FITS file.
Definition: GFits.cpp:1314
#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:413
FITS table abstract base class interface definition.
void copy_members(const GNodeArray &array)
Copy class members.
Definition: GNodeArray.cpp:932