tsgIndexSets.hpp

/*
 * Copyright (c) 2017, Miroslav Stoyanov
 *
 * This file is part of
 * Toolkit for Adaptive Stochastic Modeling And Non-Intrusive ApproximatioN: TASMANIAN
 *
 * Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions
 *    and the following disclaimer in the documentation and/or other materials provided with the distribution.
 *
 * 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse
 *    or promote products derived from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
 * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
 * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * UT-BATTELLE, LLC AND THE UNITED STATES GOVERNMENT MAKE NO REPRESENTATIONS AND DISCLAIM ALL WARRANTIES, BOTH EXPRESSED AND IMPLIED.
 * THERE ARE NO EXPRESS OR IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, OR THAT THE USE OF THE SOFTWARE WILL NOT INFRINGE ANY PATENT,
 * COPYRIGHT, TRADEMARK, OR OTHER PROPRIETARY RIGHTS, OR THAT THE SOFTWARE WILL ACCOMPLISH THE INTENDED RESULTS OR THAT THE SOFTWARE OR ITS USE WILL NOT RESULT IN INJURY OR DAMAGE.
 * THE USER ASSUMES RESPONSIBILITY FOR ALL LIABILITIES, PENALTIES, FINES, CLAIMS, CAUSES OF ACTION, AND COSTS AND EXPENSES, CAUSED BY, RESULTING FROM OR ARISING OUT OF,
 * IN WHOLE OR IN PART THE USE, STORAGE OR DISPOSAL OF THE SOFTWARE.
 */
 
#ifndef __TASMANIAN_SPARSE_GRID_INDEX_SETS_HPP
#define __TASMANIAN_SPARSE_GRID_INDEX_SETS_HPP
 
#include "tsgIOHelpers.hpp"


 
namespace TasGrid{

template<typename T>
std::vector<T> spltVector2D(std::vector<T> const &x, size_t stride, int ibegin, int iend){
    size_t sbegin(ibegin), send(iend);
    size_t num_strips = x.size() / stride;
    size_t new_stride = send - sbegin;
    std::vector<T> result(num_strips * new_stride);
    auto ix = x.begin();
    auto ir = result.begin();
    for(size_t i=0; i<num_strips; i++){
        std::copy_n(ix + sbegin, new_stride, ir);
        std::advance(ix, stride);
        std::advance(ir, new_stride);
    }
    return result;
}

template<typename T>
class Data2D{
public:
    Data2D() = default;
    template<typename IntTypeA, typename IntTypeB>
    Data2D(IntTypeA new_stride, IntTypeB new_num_strips) : stride(static_cast<size_t>(new_stride)), num_strips(static_cast<size_t>(new_num_strips)),
                                                           vec(Utils::size_mult(stride, num_strips)){}
    template<typename IntTypeA, typename IntTypeB>
    Data2D(IntTypeA new_stride, IntTypeB new_num_strips, T val) : stride(static_cast<size_t>(new_stride)),
                                                                  num_strips(static_cast<size_t>(new_num_strips)),
                                                                  vec(Utils::size_mult(stride, num_strips), val){}
    template<typename IntTypeA, typename IntTypeB>
    Data2D(IntTypeA new_stride, IntTypeB new_num_strips, std::vector<T> &&data) : stride(static_cast<size_t>(new_stride)), num_strips(static_cast<size_t>(new_num_strips)),
        vec(std::forward<std::vector<T>>(data)){}

    template<bool iomode, IO::IOPad pad>
    void writeVector(std::ostream &os) const{ IO::writeVector<iomode, pad, T>(vec, os); }

    bool empty() const{ return (num_strips == 0); }

    Data2D<T> splitData(int ibegin, int iend) const{
        if (stride == 0) return Data2D<T>(); // empty data object
        Data2D<T> result(iend - ibegin, 0);
        result.num_strips = num_strips;
        result.vec = spltVector2D(vec, stride, ibegin, iend);
        return result;
    }

    T* getStrip(int i){ return &(vec[i*stride]); }
    T const* getStrip(int i) const{ return &(vec[i*stride]); }
    typename std::vector<T>::iterator getIStrip(int i){ return vec.begin() + Utils::size_mult(stride, i); }
    size_t getStride() const{ return stride; }
    int getNumStrips() const{ return (int) num_strips; }
    size_t getTotalEntries() const{ return vec.size(); }
    T* data(){ return vec.data(); }
    T const* data() const{ return vec.data(); }
    size_t totalSize() const{ return vec.size(); }
    void clear(){
        stride = 0;
        num_strips = 0;
        vec = std::vector<double>();
    }

    inline typename std::vector<T> release(){ return std::move(vec); }

    inline typename std::vector<T>::iterator begin(){ return vec.begin(); }
    inline typename std::vector<T>::const_iterator begin() const{ return vec.cbegin(); }
    inline typename std::vector<T>::iterator end(){ return vec.end(); }
    inline typename std::vector<T>::const_iterator end() const{ return vec.cend(); }

    inline typename std::vector<T>::reverse_iterator rbegin(){ return vec.rbegin(); }

    void appendStrip(typename std::vector<T>::const_iterator const &x){
        vec.insert(vec.end(), x, x + stride);
        num_strips++;
    }

    void appendStrip(const std::vector<T> &x){
        appendStrip(x.begin());
    }

    void appendStrip(int pos, const std::vector<T> &x){
        vec.insert(vec.begin() + static_cast<size_t>(pos) * stride, x.begin(), x.end());
        num_strips++;
    }

    void append(Data2D<T> const &other) {
        vec.insert(vec.end(), other.vec.begin(), other.vec.end());
        num_strips += other.num_strips;
    }
 
private:
    size_t stride = 0, num_strips = 0;
    std::vector<T> vec;
};
 
namespace IO{
    template<typename iomode, typename DataType, typename IndexStride, typename IndexNumStrips>
    Data2D<DataType> readData2D(std::istream &is, IndexStride stride, IndexNumStrips num_strips){
        return Data2D<DataType>(stride, num_strips, readVector<iomode, DataType>(is, Utils::size_mult(stride, num_strips)));
    }
}

class MultiIndexSet{
public:
    MultiIndexSet() = default;
    MultiIndexSet(size_t cnum_dimensions, std::vector<int> &&new_indexes) :
        num_dimensions(cnum_dimensions), cache_num_indexes((int)(new_indexes.size() / cnum_dimensions)),
        indexes(std::move(new_indexes)){}
    MultiIndexSet(Data2D<int> const &data);
    template<typename iomode> MultiIndexSet(std::istream &is, iomode) :
        num_dimensions((size_t) IO::readNumber<iomode, int>(is)),
        cache_num_indexes(IO::readNumber<iomode, int>(is)),
        indexes(IO::readVector<iomode, int>(is, Utils::size_mult(num_dimensions, cache_num_indexes)))
        {}


    template<bool useAscii> void write(std::ostream &os) const;

    inline bool empty() const{ return indexes.empty(); }

    inline size_t getNumDimensions() const{ return num_dimensions; }
    inline int getNumIndexes() const{ return cache_num_indexes; }

    void addSortedIndexes(std::vector<int> const &addition);
    inline MultiIndexSet& operator += (MultiIndexSet const &addition){
        num_dimensions = addition.getNumDimensions();
        addSortedIndexes(addition.indexes);
        return *this;
    }

    inline std::vector<int>::const_iterator begin() const{ return indexes.cbegin(); }
    inline std::vector<int>::const_iterator end() const{ return indexes.cend(); }
    int* data(){ return indexes.data(); }
    int const* data() const{ return indexes.data(); }
    inline size_t totalSize() const{ return indexes.size(); }

    inline std::vector<int> release(){ return std::move(indexes); }

    int getSlot(const int *p) const;
    inline int getSlot(const std::vector<int> &p) const{ return getSlot(p.data()); }
    inline bool missing(const std::vector<int> &p) const{ return (getSlot(p.data()) == -1); }

    inline const int *getIndex(int i) const{ return &(indexes[static_cast<size_t>(i) * num_dimensions]); }

    inline std::vector<int> copyIndex(int i) const{
        return std::vector<int>(&indexes[static_cast<size_t>(i) * num_dimensions], &indexes[static_cast<size_t>(i) * num_dimensions] + num_dimensions);
    }

    MultiIndexSet operator -(const MultiIndexSet &substract) const;

    void removeIndex(const std::vector<int> &p);

    int getMaxIndex() const{ return (empty()) ? 0 : *std::max_element(indexes.begin(), indexes.end()); }
 
private:
    size_t num_dimensions = 0;
    int cache_num_indexes = 0;
    std::vector<int> indexes;
};

class StorageSet{
public:
    StorageSet() = default;
    StorageSet(int cnum_outputs, int cnum_values, std::vector<double> &&vals) :
        num_outputs(cnum_outputs), num_values(cnum_values), values(std::move(vals)){}
    template<typename iomode> StorageSet(std::istream &is, iomode) :
        num_outputs((size_t) IO::readNumber<iomode, int>(is)),
        num_values((size_t) IO::readNumber<iomode, int>(is)),
        values((IO::readFlag<iomode>(is)) ? IO::readVector<iomode, double>(is, Utils::size_mult(num_outputs, num_values)) : std::vector<double>())
    {}

    template<bool useAscii> void write(std::ostream &os) const;

    void resize(int cnum_outputs, int cnum_values){
        num_outputs = (size_t) cnum_outputs;
        num_values  = (size_t) cnum_values;
        values = std::vector<double>();
    }

    size_t getNumOutputs() const{ return num_outputs; }

    double const* getValues(int i) const{ return &(values[i*num_outputs]); }
    double* getValues(int i){ return &(values[i*num_outputs]); }

    void setValues(const double vals[]){ values = std::vector<double>(vals, vals + num_outputs * num_values); }
    void setValues(std::vector<double> &&vals){
        num_values = vals.size() / num_outputs;
        values = std::move(vals); // move assignment
    }

    StorageSet splitValues(int ibegin, int iend) const{
        return {iend - ibegin, (int) num_values, spltVector2D(values, num_outputs, ibegin, iend)};
    }

    inline std::vector<double>::const_iterator begin() const{ return values.cbegin(); }
    inline std::vector<double>::const_iterator end() const{ return values.cend(); }

    inline std::vector<double> release(){ return std::move(values); }

    void addValues(const MultiIndexSet &old_set, const MultiIndexSet &new_set, const double new_vals[]);

    double* data(){ return values.data(); }
    double const* data() const{ return values.data(); }
    size_t totalSize() const{ return values.size(); }
 
private:
    size_t num_outputs = 0, num_values = 0; // kept as size_t to avoid conversions in products, but each one is small individually
    std::vector<double> values;
};
 
}
 
#endif