GBounds.cpp

Go to the documentation of this file.

/***************************************************************************
 *                      GBounds.cpp - Boundary class                       *
 * ----------------------------------------------------------------------- *
 *  copyright (C) 2026 by Juergen Knoedlseder                              *
 * ----------------------------------------------------------------------- *
 *                                                                         *
 *  This program is free software: you can redistribute it and/or modify   *
 *  it under the terms of the GNU General Public License as published by   *
 *  the Free Software Foundation, either version 3 of the License, or      *
 *  (at your option) any later version.                                    *
 *                                                                         *
 *  This program is distributed in the hope that it will be useful,        *
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
 *  GNU General Public License for more details.                           *
 *                                                                         *
 *  You should have received a copy of the GNU General Public License      *
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
 *                                                                         *
 ***************************************************************************/
/**
 * @file GBounds.cpp
 * @brief Boundary class implementation
 * @author Juergen Knoedlseder
 */
 
/* __ Includes ___________________________________________________________ */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <cmath>
#include "GException.hpp"
#include "GTools.hpp"
#include "GFilename.hpp"
#include "GBounds.hpp"
#include "GFits.hpp"
#include "GFitsTable.hpp"
#include "GFitsBinTable.hpp"
#include "GFitsTableDoubleCol.hpp"
 
/* __ Method name definitions ____________________________________________ */
#define G_REMOVE_INX                                  "GBounds::remove(int&)"
#define G_SET   "GBounds::set(int&, double&, double&, std::string&, double&)"
#define G_MIN_SET                               "GBounds::min(int&, double&)"
#define G_MAX_SET                               "GBounds::max(int&, double&)"
#define G_MIN_GET                                        "GBounds::min(int&)"
#define G_MAX_GET                                        "GBounds::max(int&)"
#define G_MEAN                                          "GBounds::mean(int&)"
#define G_LOGMEAN                                    "GBounds::logmean(int&)"
#define G_WIDTH                                        "GBounds::width(int&)"
#define G_INSERT_VALUE        "GBounds::insert_value(int&, double&, double&)"
 
/* __ Macros _____________________________________________________________ */
 
/* __ Coding definitions _________________________________________________ */
 
/* __ Debug definitions __________________________________________________ */
 
 
/*==========================================================================
 =                                                                         =
 =                         Constructors/destructors                        =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Void constructor
 *
 * Constructs empty boundaries.
 ***************************************************************************/
GBounds::GBounds(void)
{
    // Initialise members
    init_members();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief FITS file constructor
 *
 * @param[in] filename FITS file name.
 *
 * Constructs boundaries from a FITS file.
 ***************************************************************************/
GBounds::GBounds(const GFilename& filename)
{
    // Initialise members
    init_members();
 
    // Load FITS file
    load(filename);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Array constructor
 *
 * @param[in] array Array.
 * @param[in] unit Boundary units.
 *
 * Constructs boundaries from an array of values. Optionally boundary units
 * may be specified for the array values.
 ***************************************************************************/
GBounds::GBounds(const std::vector<double>& array, const std::string& unit)
{
    // Initialise members
    init_members();
 
    // Set boundaries from array
    set(array);
 
    // Set units
    this->unit(unit);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Copy constructor
 *
 * @param[in] bounds Boundaries.
 ***************************************************************************/
GBounds::GBounds(const GBounds& bounds)
{
    // Initialise class members for clean destruction
    init_members();
 
    // Copy members
    copy_members(bounds);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Interval constructor
 *
 * @param[in] min Minimum of the interval.
 * @param[in] max Maximum of the interval.
 * @param[in] unit Boundary units.
 *
 * Constructs boundaries for one [@p min, @p max[ interval. Optionally
 * boundary units may be specified for the interval values.
 ***************************************************************************/
GBounds::GBounds(const double& min, const double& max, const std::string& unit)
{
    // Initialise members
    init_members();
 
    // Append interval
    append(min, max);
 
    // Set units
    this->unit(unit);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Boundary constructor
 *
 * @param[in] num Number of intervals.
 * @param[in] min Minimum value of first interval.
 * @param[in] max Maximum value of last interval.
 * @param[in] method Spacing method (one of "LIN", "LOG" or "POW").
 * @param[in] unit Boundary units.
 * @param[in] gamma Power law index for @p POW method.
 *
 * Constructs boundaries by defining @p num successive intervals between
 * @p min and @p max. The @p method parameter controls the spacing of the
 * boundaries. Optionally boundary units may be specified for the interval
 * values. For the @p POW method, the parameter @p gamma controls the slope
 * of the power law. See the set() method for more information.
 ***************************************************************************/
GBounds::GBounds(const int&         num,
                 const double&      min,
                 const double&      max,
                 const std::string& method,
                 const std::string& unit,
                 const double&      gamma)
{
    // Initialise members
    init_members();
 
    // Set intervals
    set(num, min, max, method, gamma);
 
    // Set units
    this->unit(unit);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Destructor
 ***************************************************************************/
GBounds::~GBounds(void)
{
    // Free members
    free_members();
 
    // Return
    return;
}
 
 
/*==========================================================================
 =                                                                         =
 =                              Operators                                  =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Assignment operator
 *
 * @param[in] bounds Boundaries to be assigned.
 * @return Boundaries.
 ***************************************************************************/
GBounds& GBounds::operator=(const GBounds& bounds)
{
    // Execute only if object is not identical
    if (this != &bounds) {
 
        // Free members
        free_members();
 
        // Initialise private members for clean destruction
        init_members();
 
        // Copy members
        copy_members(bounds);
 
    } // endif: object was not identical
 
    // Return this object
    return *this;
}
 
 
/*==========================================================================
 =                                                                         =
 =                             Public methods                              =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Clear boundaries
 ***************************************************************************/
void GBounds::clear(void)
{
    // Free members
    free_members();
 
    // Initialise private members
    init_members();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Clone boundaries
 *
 * @return Pointer to deep copy of boundaries.
 ***************************************************************************/
GBounds* GBounds::clone(void) const
{
    return new GBounds(*this);
}
 
 
/***********************************************************************//**
 * @brief Append interval
 *
 * @param[in] min Minimum value of interval.
 * @param[in] max Maximum value of interval.
 *
 * Appends an interval to the end of the boundaries container.
 ***************************************************************************/
void GBounds::append(const double& min, const double& max)
{
    // Append interval at the end of the container
    insert_interval(size(), min, max);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Insert interval
 *
 * @param[in] min Minimum value of interval.
 * @param[in] max Maximum value of interval.
 *
 * Inserts an interval into the boundaries after the first boundary that has
 * a minimum value smaller than @p min. The method implicitely assumes that
 * the intervals are ordered by increasing minimum value.
 ***************************************************************************/
void GBounds::insert(const double& min, const double& max)
{
    // Determine index at which interval should be inserted
    int inx = 0;
    for (; inx < size(); ++inx) {
        if (min < m_min[inx])
            break;
    }
 
    // Insert interval
    insert_interval(inx, min, max);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Remove interval
 *
 * @param[in] index Interval index [0,...,size()[.
 *
 * Removes interval at @p index from the boundaries container. All intervals
 * after the specified @p index are moved forward by one position.
 ***************************************************************************/
void GBounds::remove(const int& index)
{
    #if defined(G_RANGE_CHECK)
    // If index is outside boundary then throw an error
    if (index < 0 || index >= size()) {
        throw GException::out_of_range(G_REMOVE_INX, "Interval index",
                                       index, size());
    }
    #endif
 
    // Erase boundary from container
    m_min.erase(m_min.begin() + index);
    m_max.erase(m_max.begin() + index);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Reserve space for intervals
 *
 * @param[in] num Number of elements.
 ***************************************************************************/
void GBounds::reserve(const int& num)
{
    // Reserve space
    m_min.reserve(num);
    m_max.reserve(num);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Append boundaries
 *
 * @param[in] ebds Energy boundaries.
 *
 * Append boundaries to the container.
 ***************************************************************************/
void GBounds::extend(const GBounds& bounds)
{
    // Do nothing if boundaries are empty
    if (!bounds.is_empty()) {
 
        // Get size. Note that we extract the size first to avoid an
        // endless loop that arises when a container is appended to
        // itself.
        int num = bounds.size();
 
        // Reserve enough space
        reserve(size() + num);
 
        // Loop over all elements and append them to container
        for (int i = 0; i < num; ++i) {
            m_min.push_back(m_min[i]);
            m_max.push_back(m_max[i]);
        }
 
    } // endif: boundaries were not empty
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Set boundaries from array
 *
 * @param[in] array Array.
 *
 * Sets boundaries from an array. Each two subsequent values in the array
 * will form a boundary. This means that n values will lead to n-1 boundaries
 * with the following mapping:
 *
 *      [array[0], array[1]]
 *      [array[1], array[2]]
 *      ...
 *      [array[n-2], array[n-1]]
 *
 * If there is only one value in the array the following empty boundary will
 * be appended:
 *
 *      [array[0], array[0]]
 *
 * Note that the method conserves the units of the boundaries.
 ***************************************************************************/
void GBounds::set(const std::vector<double>& array)
{
    // Converse boundary units
    std::string unit = m_unit;
 
    // Initialise members
    clear();
 
    // Put back boundary units
    m_unit = unit;
 
    // Get number of values in container
    int num = array.size();
 
    // If there is only one value in the container then append an empty
    // boundary
    if (num == 1) {
        append(array[0], array[0]);
    }
 
    // ... otherwise if there is more than one value in the container
    // then append subsequent values as boundaries
    else if (num > 1) {
        for (int i = 0; i < num-1; ++i) {
            append(array[i], array[i+1]);
        }
    }
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Set intervals
 *
 * @param[in] num Number of intervals.
 * @param[in] min Minimum value of first interval.
 * @param[in] max Maximum value of last interval.
 * @param[in] method Energy spacing method (one of "LIN", "LOG" or "POW").
 * @param[in] gamma Power law index for @p POW method.
 *
 * @exception GException::invalid_argument
 *            @p num is negative.
 *            @p min value is not smaller than @p max value.
 *            Non-positive @p min value specified for LOG or POW method.
 *            Invalid spacing @p method specified.
 *
 * Sets boundaries by defining @p num successive intervals between @p min
 * and @p max. The @p method parameter controls the spacing of the
 * boundaries. Note that the @p max value needs to be larger than the @p min
 * value. For the "LOG" and "POW" methods the @p min value needs to be
 * positive.
 *
 * Note that the method conserves the units of the boundaries.
 ***************************************************************************/
void GBounds::set(const int&         num,
                  const double&      min,
                  const double&      max,
                  const std::string& method,
                  const double&      gamma)
{
    // Converse boundary units
    std::string unit = m_unit;
 
    // Initialise members
    clear();
 
    // Put back boundary units
    m_unit = unit;
 
    // Convert method to upper-case string
    std::string umethod = gammalib::toupper(method);
 
    // Check validity of number of values
    if (num < 0) {
        std::string msg = "Negative number of values "+gammalib::str(num)+
                          " specified. Please specify a non-negative number "
                          "of values.";
        throw GException::invalid_argument(G_SET, msg);
    }
 
    // Check validity of minimum value for logarithmic and powerlaw boundaries
    if ((umethod == "LOG") || (umethod == "POW")) {
        if (min <= 0.0) {
            std::string msg = "Minimum value "+gammalib::str(min)+" is not "
                              "positive. \""+method+"\" method requires a "
                              "positive minimum value. Please specify a "
                              "positive minimum value.";
            throw GException::invalid_argument(G_SET, msg);
        }
    }
 
    // Check validity of minimum and maximum values
    if (min > max) {
        std::string msg = "Minimum value "+gammalib::str(min)+" is larger than "
                          "maximum value "+gammalib::str(max)+". Please specify "
                          "a minimum value that does not exceed the "
                          "maximum value.";
        throw GException::invalid_argument(G_SET, msg);
    }
 
    // Reserve elements
    reserve(num);
 
    // Set linear boundaries
    if (umethod == "LIN") {
        double bin   = (max - min)/double(num);
        double value = min;
        for (int i = 0; i < num; ++i) {
            append(value, value+bin);
            value += bin;
        }
    }
 
    // ... otherwise set logarithmic boundaries
    else if (umethod == "LOG") {
        double logmin = std::log10(min);
        double logmax = std::log10(max);
        double logbin = (logmax - logmin)/double(num);
        double vmin   = min;
        for (int i = 0; i < num; ++i) {
            double vmax = std::pow(10.0, double(i+1)*logbin + logmin);
            append(vmin, vmax);
            vmin = vmax;
        }
    }
 
    // ... otherwise set power law boundaries
    else if (umethod == "POW") {
        double a    = 1.0 - gamma;
        double c    = (a == 0.0) ? 1.0 / (std::log(max)   - std::log(min))
                                 : a   / (std::pow(max,a) - std::pow(min,a));
        double b    = double(c*(num));
        double vmin = min;
        for (int i = 0; i < num; ++i) {
            double logvmax = (a == 0.0) ? 1.0/b + std::log(vmin)
                                        : std::log(a/b + std::pow(vmin,a)) / a;
            double vmax    = std::exp(logvmax);
            append(vmin, vmax);
            vmin = vmax;
        }
    }
 
    // ... otherwise throw an exception
    else {
        std::string msg = "Invalid spacing method \""+umethod+"\" "
                          "specified. Please provide one of \"LIN\", \"LOG\""
                          ", or \"POW\".";
        throw GException::invalid_argument(G_SET, msg);
    }
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Load boundaries from FITS file
 *
 * @param[in] filename FITS file name.
 *
 * Loads boundaries from a FITS file.
 *
 * If no extension name is provided, the boundaries are loaded from the
 * `BOUNDS` extension.
 ***************************************************************************/
void GBounds::load(const GFilename& filename)
{
    // Open FITS file
    GFits fits(filename);
 
    // Get boundary table
    const GFitsTable& table = *fits.table(filename.extname(gammalib::extname_bounds));
 
    // Read boundaries from table
    read(table);
 
    // Close FITS file
    fits.close();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Save boundaries into FITS file
 *
 * @param[in] filename FITS file name.
 * @param[in] clobber Overwrite an existing boundaries extension?
 *
 * Saves boundaries into a FITS file. If a file with the given @p filename
 * does not yet exist it will be created, otherwise the method opens the
 * existing file. Boundaries can only be appended to an existing file if the
 * @p clobber flag is set to "true" (otherwise an exception is thrown).
 *
 * The method will append a binary FITS table containing the boundaries to
 * the FITS file. The extension name can be specified as part of the
 * @p filename. For example the @p filename
 *
 *      myfile.fits[BOUNDARIES]
 *
 * will save the boundaries in the `BOUNDARIES` extension of the
 * "myfile.fits" file. If the extension exists already in the file it will
 * be replaced, otherwise a new extension will be created. If no extension
 * name is provided, the method will use `BOUNDS` as the default extension
 * name for boundaries.
 ***************************************************************************/
void GBounds::save(const GFilename& filename,
                   const bool&      clobber) const
{
    // Open or create FITS file (without extension name since the requested
    // extension may not yet exist in the file)
    GFits fits(filename.url(), true);
 
    // Write boundaries to FITS file
    write(fits, filename.extname(gammalib::extname_bounds));
 
    // Save to file
    fits.save(clobber);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read boundaries from FITS table
 *
 * @param[in] table FITS table.
 *
 * Reads the boundaries from a FITS table.
 ***************************************************************************/
void GBounds::read(const GFitsTable& table)
{
    // Clear members
    clear();
 
    // Get number of boundaries in file
    int num = table.nrows();
 
    // Continue only if there are boundaries
    if (num > 0) {
 
        // Get columns
        const GFitsTableCol* col_min = table["MIN"];
        const GFitsTableCol* col_max = table["MAX"];
 
        // Extract unit
        m_unit = col_min->unit();
 
        // Reserve space for boundaries
        reserve(num);
 
        // Read boundaries
        for (int i = 0; i < num; ++i) {
            append(col_min->real(i), col_max->real(i));
        }
 
    }
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Write boundaries into FITS object
 *
 * @param[in] fits FITS file.
 * @param[in] extname Boundary extension name.
 *
 * Writes the boundaries into a FITS object.
 ***************************************************************************/
void GBounds::write(GFits& fits, const std::string& extname) const
{
    // Create boundary columns
    GFitsTableDoubleCol col_min("MIN", size());
    GFitsTableDoubleCol col_max("MAX", size());
 
    // Fill boundary columns
    for (int i = 0; i < size(); ++i) {
        col_min(i) = m_min[i];
        col_max(i) = m_max[i];
    }
 
    // Set units
    col_min.unit(m_unit);
    col_max.unit(m_unit);
 
    // Create binary table
    GFitsBinTable table;
 
    // Append columns
    table.append(col_min);
    table.append(col_max);
 
    // Set extension name
    table.extname(extname);
 
    // If the FITS object contains already an extension with the same
    // name then remove now this extension
    if (fits.contains(extname)) {
        fits.remove(extname);
    }
 
    // Append boundary table to FITS file
    fits.append(table);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Returns bin index for a value
 *
 * @param[in] value Value.
 * @return Bin index.
 *
 * Returns the boundary bin index for a given value. By convention, the limits
 * for a bin are defined as
 *
 *      min <= value < max
 *
 * i.e. a value equaling to max falls above the boundary.
 *
 * If the value falls outside all boundaries, -1 is returned.
 ***************************************************************************/
int GBounds::index(const double& value) const
{
    // Initialise index with 'not found'
    int index = -1;
 
    // Search all boundaries for containment
    for (int i = 0; i < size(); ++i) {
        if (value >= m_min[i] && value < m_max[i]) {
            index = i;
            break;
        }
    }
 
    // Return index
    return index;
}
 
 
/***********************************************************************//**
 * @brief Return minimum value of all intervals
 *
 * @return Minimum value of all intervals.
 ***************************************************************************/
double GBounds::min(void) const
{
    // Initialise minimum value
    double min = 0.0;
 
    // If there are intervals then determine their minimum value
    if (size() > 0) {
        min = m_min[0];
        for (int i = 1; i < size(); ++i) {
            if (m_min[i] < min) {
                min = m_min[i];
            }
        }
    }
 
    // Return minimum value
    return min;
}
 
 
/***********************************************************************//**
 * @brief Return maximum value of all intervals
 *
 * @return Maximum value of all intervals.
 ***************************************************************************/
double GBounds::max(void) const
{
    // Initialise maximum value
    double max = 0.0;
 
    // If there are intervals then determine their maximum value
    if (size() > 0) {
        max = m_max[0];
        for (int i = 1; i < size(); ++i) {
            if (m_max[i] > max) {
                max = m_max[i];
            }
        }
    }
 
    // Return maximum value
    return max;
}
 
 
/***********************************************************************//**
 * @brief Set minimum value for a given interval
 *
 * @param[in] index Energy interval index (0,...,size()-1).
 * @param[in] value Minimum value of interval.
 *
 * @exception GException::out_of_range
 *            Specified index is out of range.
 *
 * Sets the minimum value for the interval @p index.
 ***************************************************************************/
void GBounds::min(const int& index, const double& value)
{
    #if defined(G_RANGE_CHECK)
    // Throw an exception if index is outside valid range
    if (index < 0 || index >= size()) {
        throw GException::out_of_range(G_MIN_SET, "Boundary index",
                                       index, size());
    }
    #endif
 
    // Set minimum value
    m_min[index] = value;
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Set maximum value for a given interval
 *
 * @param[in] index Interval index (0,...,size()-1).
 * @param[in] value Maximum value of interval.
 *
 * @exception GException::out_of_range
 *            Specified index is out of range.
 *
 * Sets the maximum value for the interval @p index.
 ***************************************************************************/
void GBounds::max(const int& index, const double& value)
{
    #if defined(G_RANGE_CHECK)
    // Throw an exception if index is outside valid range
    if (index < 0 || index >= size()) {
        throw GException::out_of_range(G_MAX_SET, "boundary index",
                                       index, size());
    }
    #endif
 
    // Set maximum value
    m_max[index] = value;
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Returns minimum value for a given interval
 *
 * @param[in] index Interval index (0,...,size()-1).
 * @return Minimum value of interval.
 *
 * @exception GException::out_of_range
 *            Specified index is out of range.
 ***************************************************************************/
double GBounds::min(const int& index) const
{
    #if defined(G_RANGE_CHECK)
    // Throw an exception if index is outside valid range
    if (index < 0 || index >= size()) {
        throw GException::out_of_range(G_MIN_GET, "Boundary index",
                                       index, size());
    }
    #endif
 
    // Return minimum value
    return (m_min[index]);
}
 
 
/***********************************************************************//**
 * @brief Returns maximum value for a given interval
 *
 * @param[in] index Interval index (0,...,size()-1).
 * @return Maximum value of interval.
 *
 * @exception GException::out_of_range
 *            Specified index is out of range.
 ***************************************************************************/
double GBounds::max(const int& index) const
{
    #if defined(G_RANGE_CHECK)
    // Throw an exception if index is outside valid range
    if (index < 0 || index >= size()) {
        throw GException::out_of_range(G_MAX_GET, "Boundary index",
                                       index, size());
    }
    #endif
 
    // Return maximum value
    return (m_max[index]);
}
 
 
/***********************************************************************//**
 * @brief Returns mean value for a given interval
 *
 * @param[in] index Interval index (0,...,size()-1).
 * @return Mean value of interval.
 *
 * @exception GException::out_of_range
 *            Specified index is out of range.
 *
 * Computes the mean value for the interval @p index.
 ***************************************************************************/
double GBounds::mean(const int& index) const
{
    #if defined(G_RANGE_CHECK)
    // If index is outside boundary then throw an error
    if (index < 0 || index >= size()) {
        throw GException::out_of_range(G_MEAN, "Boundary index",
                                       index, size());
    }
    #endif
 
    // Compute mean value
    double mean = 0.5 * (m_min[index] + m_max[index]);
 
    // Return mean value
    return mean;
}
 
 
/***********************************************************************//**
 * @brief Returns logarithmic mean value for a given interval
 *
 * @param[in] index Interval index (0,...,size()-1).
 * @return Logarithmic mean value of interval.
 *
 * @exception GException::out_of_range
 *            Specified index is out of range.
 *
 * Computes the logarithmic mean value for the interval @p index.
 ***************************************************************************/
double GBounds::logmean(const int& index) const
{
    #if defined(G_RANGE_CHECK)
    // If index is outside boundary then throw an error
    if (index < 0 || index >= size()) {
        throw GException::out_of_range(G_LOGMEAN, "Boundary index",
                                       index, size());
    }
    #endif
 
    // Compute logarithmic mean value
    double logmean = std::sqrt(m_min[index] * m_max[index]);
 
    // Return logarithmic mean value
    return logmean;
}
 
 
/***********************************************************************//**
 * @brief Returns interval width
 *
 * @param[in] index Interval index (0,...,size()-1).
 * @return Interval width.
 *
 * @exception GException::out_of_range
 *            Specified index is out of range.
 ***************************************************************************/
double GBounds::width(const int& index) const
{
    #if defined(G_RANGE_CHECK)
    // If index is outside boundary then throw an error
    if (index < 0 || index >= size()) {
        throw GException::out_of_range(G_WIDTH, "Boundary index",
                                       index, size());
    }
    #endif
 
    // Return
    return (m_max[index]-m_min[index]);
}
 
 
/***********************************************************************//**
 * @brief Checks whether boundaries contain value
 *
 * @param[in] value Value to be checked.
 * @return True if value falls in at least one interval, false otherwise.
 *
 * Checks if the @p value falls in at least one of the intervals. The method
 * exits when the first matching interval has been found.
 ***************************************************************************/
bool GBounds::contains(const double& value) const
{
    // Initialise result
    bool contained = false;
 
    // Test all boundaries
    for (int i = 0; i < size(); ++i) {
        if (value >= m_min[i] && value <= m_max[i]) {
            contained = true;
            break;
        }
    }
 
    // Return result
    return contained;
}
 
 
/***********************************************************************//**
 * @brief Checks whether boundaries contain an interval
 *
 * @param[in] min Minimum value of interval to be checked.
 * @param[in] max Maximum value of interval to be checked.
 * @return True if interval [min, max] exists within boundariesi.
 *
 * Checks if the interval [@p min, @p max] exists within the boundary
 * contained.
 ***************************************************************************/
bool GBounds::contains(const double& min, const double& max) const
{
    // Initialise result
    bool contained = false;
 
    // Test all boundaries
    for (int i = 0; i < size(); ++i) {
        if (min == m_min[i] && max == m_max[i]) {
            contained = true;
            break;
        }
    }
 
    // Return result
    return contained;
}
 
 
/***********************************************************************//**
 * @brief Print boundaries
 *
 * @param[in] chatter Chattiness.
 * @return String containing boundary information.
 ***************************************************************************/
std::string GBounds::print(const GChatter& chatter) const
{
    // Initialise result string
    std::string result;
 
    // Continue only if chatter is not silent
    if (chatter != SILENT) {
 
        // Append Header
        result.append("=== GBounds ===");
 
        // Append information
        result.append("\n"+gammalib::parformat("Number of intervals"));
        result.append(gammalib::str(size()));
        result.append("\n"+gammalib::parformat("Value range"));
        result.append(gammalib::str(min()));
        result.append(" - ");
        result.append(gammalib::str(max()));
        if (unit().length() > 0) {
            result.append(" ");
            result.append(unit());
        }
 
        // If there are multiple bins then append them
        if (chatter >= EXPLICIT) {
            if (size() > 1) {
                for (int i = 0; i < size(); ++i) {
                    result.append("\n");
                    result.append(gammalib::parformat("Interval "+
                                  gammalib::str(i)));
                    result.append(gammalib::str(min(i)));
                    result.append(" - ");
                    result.append(gammalib::str(max(i)));
                }
            }
        } // endif: chatter was less than explicit
 
    } // endif: chatter was not silent
 
    // Return
    return result;
}
 
 
/*==========================================================================
 =                                                                         =
 =                            Private methods                              =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Initialise class members
 ***************************************************************************/
void GBounds::init_members(void)
{
    // Initialise members
    m_unit.clear();
    m_min.clear();
    m_max.clear();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Copy class members
 *
 * @param[in] ebds Energy boundaries.
 ***************************************************************************/
void GBounds::copy_members(const GBounds& bounds)
{
    // Copy members
    m_unit = bounds.m_unit;
    m_min  = bounds.m_min;
    m_max  = bounds.m_max;
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Delete class members
 ***************************************************************************/
void GBounds::free_members(void)
{
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Insert interval
 *
 * @param[in] index Index after with interval is inserted.
 * @param[in] min Minimum value of interval.
 * @param[in] max Maximum value of interval.
 *
 * @exception GException::invalid_argument
 *            Minimum value larger than maximum value
 *
 * Inserts an interval after the specified @p index in the  boundaries. The
 * method does not reorder the intervals, instead the client needs to
 * determine the approriate @p index.
 *
 * Invalid parameters do not produce any exception, but are handled
 * transparently. If the interval is invalid (i.e. @p min > @p max) an
 * exception is thrown. If the @p index is out of the valid range, the
 * index will be adjusted to either the first or the last element.
 ***************************************************************************/
void GBounds::insert_interval(const int&    index,
                              const double& min,
                              const double& max)
{
    // Throw an exception if interval is invalid
    if (min > max) {
        std::string msg = "Invalid interval specified. Minimum"
                          " value "+gammalib::str(min)+" can not be"
                          " larger than maximum value "+
                          gammalib::str(max)+".";
        throw GException::invalid_argument(G_INSERT_VALUE, msg);
    }
 
    // Inserts interval
    m_min.insert(m_min.begin()+index, min);
    m_max.insert(m_max.begin()+index, max);
 
    // Return
    return;
}
 
 
/*==========================================================================
 =                                                                         =
 =                                Friends                                  =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Boundaries equality operator friend
 *
 * @param[in] a First boundaries.
 * @param[in] b Second boundaries.
 * @return True if both boundaries are identical.
 ***************************************************************************/
bool operator==(const GBounds& a, const GBounds& b)
{
    // Initialise identify flag
    bool identity = true;
 
    // Check that both boundaries have the same size
    if (a.size() != b.size()) {
        identity = false;
    }
 
    // Check all boundaries
    else {
        for (int i = 0; i < a.size(); ++i) {
            if ((a.min(i) != b.min(i)) || (a.max(i) != b.max(i))) {
                identity = false;
                break;
            }
        }
    }
 
    // Return identity flag
    return identity;
}