GVectorSparse.cpp

Go to the documentation of this file.

/***************************************************************************
 *                  GVectorSparse.cpp - Sparse vector class                *
 * ----------------------------------------------------------------------- *
 *  copyright (C) 2024 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 GVectorSparse.cpp
 * @brief Sparse vector class implementation
 * @author Juergen Knoedlseder
 */
 
/* __ Includes ___________________________________________________________ */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include "GVector.hpp"
#include "GVectorSparse.hpp"
#include "GTools.hpp"
 
 
/* __ Method name definitions ____________________________________________ */
#define G_CONSTRUCTOR                    "GVectorSparse::GVectorSparse(int&)"
#define G_OPERATOR                          "GVectorSparse::operator[](int&)"
 
 
/*==========================================================================
 =                                                                         =
 =                         Constructors/destructors                        =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Void sparse vector constructor
 ***************************************************************************/
GVectorSparse::GVectorSparse(void)
{
    // Initialise class members
    init_members();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Sparse vector constructor
 *
 * @param[in] num Vector size.
 * @param[in] alloc Number of elements for allocation.
 *
 * Constructs a sparse vector, specifying the size of the expanded vector
 * @p num and the number of elements @p alloc for which memory should be
 * allocated. Using pre-allocated memory reduces the dynamic reallocation
 * of memory.
 ***************************************************************************/
GVectorSparse::GVectorSparse(const int& num, const int& alloc)
{
    // Initialise class members
    init_members();
 
    // Store vector size
    m_num = num;
 
    // Allocate memory
    alloc_members(alloc);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Vector conversion constructor
 *
 * @param[in] vector Vector.
 *
 * Converts a regular vector into a sparse vector by stripping all zero
 * elements.
 ***************************************************************************/
GVectorSparse::GVectorSparse(const GVector& vector)
{
    // Initialise class members
    init_members();
 
    // Set attributes
    m_num      = vector.size();
    m_elements = vector.non_zeros();
 
    // Allocate memory
    alloc_members(m_elements);
 
    // Copy non-zero elements
    for (int i = 0, idst = 0; i < vector.size(); ++i) {
        if (vector[i] != 0.0) {
            m_inx[idst]  = i;
            m_data[idst] = vector[i];
            idst++;
        }
    }
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Copy constructor
 *
 * @param[in] vector Sparse vector.
 ***************************************************************************/
GVectorSparse::GVectorSparse(const GVectorSparse& vector)
{
    // Initialise class members
    init_members();
 
    // Copy members
    copy_members(vector);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Destructor
 ***************************************************************************/
GVectorSparse::~GVectorSparse(void)
{
    // Free members
    free_members();
 
    // Return
    return;
}
 
 
/*==========================================================================
 =                                                                         =
 =                               Operators                                 =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Assignment operator
 *
 * @param[in] vector Sparse vector.
 * @return Sparse vector.
 ***************************************************************************/
GVectorSparse& GVectorSparse::operator=(const GVectorSparse& vector)
{
    // Execute only if object is not identical
    if (this != &vector) {
 
        // Free members
        free_members();
 
        // Initialise private members
        init_members();
 
        // Copy members
        copy_members(vector);
 
    } // endif: object was not identical
 
    // Return this object
    return *this;
}
 
 
/***********************************************************************//**
 * @brief Equality operator
 *
 * @param[in] vector Sparse vector.
 * @return True if sparse vectors are identical.
 *
 * Returns true if both sparse vectors are identical. Sparse vectors are
 * considered identical if they have the same size and if all their elements
 * are identical.
 ***************************************************************************/
bool GVectorSparse::operator==(const GVectorSparse& vector) const
{
    // Initalise result depending on sparse vector size and number of
    // elements. If the test succeeds we can simply check all elements.
    bool result = ((m_num      == vector.m_num) &&
                   (m_elements == vector.m_elements));
 
    // Test for difference. Break at first difference
    if (result) {
        for (int i = 0; i < m_elements; ++i) {
            if ((m_inx[i]  != vector.m_inx[i]) ||
                (m_data[i] != vector.m_data[i])) {
                result = false;
                break;
            }
        }
    }
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Non-equality operator
 *
 * @param[in] vector Sparse vector.
 * @return True if both sparse vectors are different.
 ***************************************************************************/
bool GVectorSparse::operator!=(const GVectorSparse& vector) const
{
    // Get negated result of equality operation
    bool result = !(this->operator==(vector));
 
    // Return result
    return result;
}
 
 
/*==========================================================================
 =                                                                         =
 =                             Public methods                              =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Clear sparse vector
 ***************************************************************************/
void GVectorSparse::clear(void)
{
    // Free members
    free_members();
 
    // Initialise private members
    init_members();
 
    // Return
    return; 
}
 
 
/***********************************************************************//**
 * @brief Clone sparse vector
 *
 * @return Pointer to deep copy of sparse vector.
 ***************************************************************************/
GVectorSparse* GVectorSparse::clone(void) const
{
    // Clone sparse vector
    return new GVectorSparse(*this);
}
 
 
/***********************************************************************//**
 * @brief Sparse vector element access with range checking
 *
 * @param[in] index Element index [0,...,size()[.
 * @return Reference to vector element.
 *
 * @exception GException::out_of_range
 *            Element index is out of range.
 *
 * Returns the sparse vector element on the basis of the @p index of the
 * full vector. If no element exists yet for this @p index, a element with
 * value of zero will be added to the vector. Note that using the method
 * will be very slow, and if no values are set, the sparse vector will be
 * populated with zero values.
 ***************************************************************************/
double& GVectorSparse::operator[](const int& index)
{
    // Throw an exception if index is out of range
    #if defined(G_RANGE_CHECK)
    if (index < 0 || index >= size()) {
        throw GException::out_of_range(G_OPERATOR, "Vector element index",
                                       index, size());
    }
    #endif
 
    // Allocate pointer to result
    double* result = NULL;
 
    // Find first position after the index
    for (int i = 0; i < m_elements; ++i) {
 
        // If element with requested index is already in vector
        // then return pointer to this element.
        if (m_inx[i] == index) {
            result = &(m_data[i]);
            break;
        }
 
        // ... otherwise search first inx after requested index
        else if (m_inx[i] > index) {
 
            // Insert index and value of zero
            insert(i, index, 0.0);
 
            // Set pointer to inserted element
            result = &(m_data[i]);
 
            // Exit loop
            break;
 
        } // endelse: searched first inx after requested index
 
    } // endfor: looped over inx array
 
    // If we have no valid pointer then append element to sparse
    // vector
    if (result == NULL) {
 
        // Save
        int num = m_elements;
 
        // Append index and value of zero
        insert(m_elements, index, 0.0);
 
        // Set pointer to appended element
        result = &(m_data[num]);
 
    }
 
    // Return reference
    return (*result);
}
 
 
/***********************************************************************//**
 * @brief Sparse vector element access with range checking
 *
 * @param[in] index Element index [0,...,size()[.
 * @return Reference to vector element.
 *
 * @exception GException::out_of_range
 *            Element index is out of range.
 *
 * Returns the sparse vector element on the basis of the @p index of the
 * full vector. If no element exists for this @p index, a static zero value
 * is returned.
 ***************************************************************************/
const double& GVectorSparse::operator[](const int& index) const
{
    // Set static zero value
    static double zero = 0.0;
 
    // Throw an exception if index is out of range
    #if defined(G_RANGE_CHECK)
    if (index < 0 || index >= size()) {
        throw GException::out_of_range(G_OPERATOR, "Vector element index",
                                       index, size());
    }
    #endif
 
    // Get sparse index
    int inx = index2inx(index);
 
    // Set result pointer
    const double* result = (inx != -1) ? &(m_data[inx]) : &zero;
 
    // Return reference
    return (*result);
}
 
 
/***********************************************************************//**
 * @brief Print sparse vector information
 *
 * @param[in] chatter Chattiness.
 * @return String containing vector information.
 ***************************************************************************/
std::string GVectorSparse::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("=== GVectorSparse ===");
 
        // Append information
        result.append("\n"+gammalib::parformat("Size of vector"));
        result.append(gammalib::str(m_num));
        result.append("\n"+gammalib::parformat("Number of elements"));
        result.append(gammalib::str(m_elements));
        result.append("\n"+gammalib::parformat("Allocated elements"));
        result.append(gammalib::str(m_alloc));
        if (m_colinx != -1) {
            result.append("\n"+gammalib::parformat("Column index"));
            result.append(gammalib::str(m_colinx));
        }
 
    } // endif: chatter was not silent
 
    // Return result
    return result;
}
 
 
/*==========================================================================
 =                                                                         =
 =                             Private methods                             =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Initialise class members
 ***************************************************************************/
void GVectorSparse::init_members(void)
{
    // Initialise members
    m_num      = 0;
    m_elements = 0;
    m_alloc    = 0;
    m_colinx   = -1;
    m_inx      = NULL;
    m_data     = NULL;
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Copy class members
 *
 * @param[in] vector Sparse vector.
 ***************************************************************************/
void GVectorSparse::copy_members(const GVectorSparse& vector)
{
    // Copy attributes
    m_num      = vector.m_num;
    m_elements = vector.m_elements;
    m_alloc    = vector.m_alloc;
    m_colinx   = vector.m_colinx;
 
    // Copy index array
    if (vector.m_inx != NULL) {
        m_inx = new int[m_alloc];
        for (int i = 0; i < m_elements; ++i) {
            m_inx[i] = vector.m_inx[i];
        }
    }
 
    // Copy data array
    if (vector.m_data != NULL) {
        m_data = new double[m_alloc];
        for (int i = 0; i < m_elements; ++i) {
            m_data[i] = vector.m_data[i];
        }
    }
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Delete class members
 ***************************************************************************/
void GVectorSparse::free_members(void)
{
    // Delete arrays
    if (m_inx  != NULL) delete [] m_inx;
    if (m_data != NULL) delete [] m_data;
 
    // Initialise pointers
    m_inx  = NULL;
    m_data = NULL;
 
    // Set allocated memory to zero
    m_alloc = 0;
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Allocate memory for elements
 *
 * @param[in] alloc Number of elements for allocation.
 *
 * Allocates memory for @p alloc elements. Memory for existing elements will
 * be deleted.
 ***************************************************************************/
void GVectorSparse::alloc_members(const int& alloc)
{
    // Delete existing elements
    if (m_inx  != NULL) delete [] m_inx;
    if (m_data != NULL) delete [] m_data;
 
    // If number of elements is positive then allocate memory
    if (alloc > 0) {
        m_inx  = new int[alloc];
        m_data = new double[alloc];
    }
 
    // Set allocated memory
    m_alloc = alloc;
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Return inx for element index
 *
 * @param[in] index Element index [0,...,size()[.
 * @return Sparse index (-1 if not found).
 *
 * Returns sparse index for a given element index. If no sparse index is
 * found the method returns -1.
 ***************************************************************************/
int GVectorSparse::index2inx(const int& index) const
{
    // Initialise inx to "not found"
    int inx = -1;
 
    // Search for index
    for (int i = 0; i < m_elements; ++i) {
        if (m_inx[i] == index) {
            inx = i;
            break;
        }
    }
 
    // Return sparse index
    return inx;
}
 
 
/***********************************************************************//**
 * @brief Insert one element into sparse vector
 *
 * @param[in] index Element index [0,...,elements()+1[.
 * @param[in] inx Vector index [0,...,size()[.
 * @param[in] data Vector data.
 *
 * Insert one element into the sparse vector. An element consists of an
 * index @p inx and a value @p data. The element is inserted before the
 * element with @p index. If @p index is equal to elements() the element is
 * appended.
 *
 * The method does only allocate fresh memory in case that the existing
 * memory is exhausted. Fresh memory is allocated in blocks of 512 elements
 * until the total vector size is filled.
 ***************************************************************************/
void GVectorSparse::insert(const int& index, const int& inx, const double& data)
{
    // Initialise pointers to elements
    int*    new_inx  = m_inx;
    double* new_data = m_data;
 
    // Determine if new memory allocation is needed
    bool allocate = (m_elements + 1 > m_alloc);
 
    // If there is not enough memory allocated then allocate new memory
    if (allocate) {
        m_alloc  = m_elements + G_SPARSE_VECTOR_DEFAULT_MEM_BLOCK;
        if (m_alloc > m_num) {
            m_alloc = m_num;
        }
        new_inx  = new int[m_alloc];
        new_data = new double[m_alloc];
    }
 
    // Copy information before inserted element. This is only needed in
    // case that new memory has been allocated.
    if (allocate) {
        for (int i = 0; i < index; ++i) {
            new_inx[i]  = m_inx[i];
            new_data[i] = m_data[i];
        }
    }
 
    // Copy information after inserted element. This needs to be done in
    // inverse order to preserve content in case that no new memory was
    // allocated.
    for (int i = m_elements; i > index; --i) {
        new_inx[i]  = m_inx[i-1];
        new_data[i] = m_data[i-1];
    }
 
    // Insert element
    new_inx[index]  = inx;
    new_data[index] = data;
 
    // If new memory was allocated then replace old memory
    if (allocate) {
 
        // Delete old memory
        delete [] m_inx;
        delete [] m_data;
 
        // Store pointers to new memory
        m_inx  = new_inx;
        m_data = new_data;
 
    } // endif: replaced old memory
 
    // Increment number of elements
    m_elements++;
 
    // Return
    return;
}