GridBase.h

#pragma once

#include <cmath>
#include <exception>
#include <vector>

namespace SSAGES
{


template<typename T>
class GridBase
{
protected:
    std::vector<T> data_;

    size_t dimension_;

    std::vector<int> numPoints_;

    std::pair< std::vector<double>,std::vector<double> > edges_;

    std::vector<bool> isPeriodic_;

    std::vector<int> wrapIndices(const std::vector<int> &indices) const
    {
        std::vector<int> newIndices(indices);
        for (size_t i=0; i<dimension_; ++i) {
            if (!GetPeriodic(i)) {
                continue;
            }

            int index = indices.at(i);
            while (index < 0) { index += GetNumPoints(i); }
            while (index >= GetNumPoints(i)) { index -= GetNumPoints(i); }
            newIndices.at(i) = index;
        }

        return newIndices;
    }


    virtual size_t mapTo1d(const std::vector<int> &indices) const = 0;

protected:

    GridBase(std::vector<int> numPoints,
             std::vector<double> lower,
             std::vector<double> upper,
             std::vector<bool> isPeriodic)
        : dimension_(numPoints.size()),
          numPoints_(numPoints),
          edges_(std::pair< std::vector<double>, std::vector<double> >(lower, upper)),
          isPeriodic_(isPeriodic)
    {
        // Check that vector sizes are correct
        if (edges_.first.size() != dimension_ ||
            edges_.second.size() != dimension_) {
            throw std::invalid_argument("Size of vector containing upper or "
                "lower edges, does not match size of vector containing "
                "number of grid points.");
        }
        if (isPeriodic_.size() == 0) {
            // Default: Non-periodic in all dimensions
            isPeriodic.resize(dimension_, false);
        } else if (isPeriodic_.size() != dimension_) {
            throw std::invalid_argument("Size of vector isPeriodic does not "
                    "match size of vector containing number of grid points.");
        }
        
    }

public:

    size_t GetDimension() const
    {
        return dimension_;
    }


    const std::vector<int>& GetNumPoints() const
    {
        return numPoints_;
    }


    int GetNumPoints(size_t dim) const
    {
        if (dim >= GetDimension()) {
            std::cerr << "Warning! Grid size requested for a dimension larger "
                         "than the grid dimensionality!\n";
            return 0;
        }

        return numPoints_.at(dim);
    }


    const std::vector<double>& GetLower() const
    {
        return edges_.first;
    }


    double GetLower(size_t dim) const
    {
        if (dim >= GetDimension()) {
            std::cerr << "Warning! Lower edge requested for a dimension larger "
                         "than the grid dimensionality!\n";
            return 0.0;
        }
        return GetLower().at(dim);
    }


    const std::vector<double>& GetUpper() const
    {
        return edges_.second;
    }



    double GetUpper(size_t dim) const
    {
        if (dim >= GetDimension()) {
            std::cerr << "Warning! Upper edge requested for a dimension larger "
                         "than the grid dimensionality!\n";
            return 0.0;
        }
        return GetUpper().at(dim);
    }


    const std::vector<bool>& GetPeriodic() const
    {
        return isPeriodic_;
    }


    bool GetPeriodic(size_t dim) const
    {
        if (dim >= GetDimension()) {
            std::cerr << "Warning! Periodicity requested for a dimension larger "
                         "than the grid dimensionality!\n";
            return false;
        }
        return GetPeriodic().at(dim);
    }


    size_t size() const
    {
        return data_.size();
    }


    T *data()
    {
        return data_.data();
    }


    T const *data() const
    {
        return data_.data();
    }


    std::vector<int> GetIndices(const std::vector<double> &x) const
    {
        // Check that input vector has the correct dimensionality
        if (x.size() != dimension_) {
            throw std::invalid_argument("Specified point has a larger "
                    "dimensionality than the grid.");
        }

        std::vector<int> indices(dimension_);
        for (size_t i = 0; i < dimension_; ++i) {
            double xpos = x.at(i);
            if (!GetPeriodic(i)) {
                if (xpos < GetLower(i)) {
                    indices.at(i) = -1;
                    continue;
                } else if (xpos > GetUpper(i)) {
                    indices.at(i) = GetNumPoints(i);
                    continue;
                }
            }

            // To make sure, the value is rounded in the correct direction.
            double spacing = (GetUpper(i) - GetLower(i)) / GetNumPoints(i);
            indices.at(i) = std::floor( (xpos - GetLower(i)) / spacing);
        }

        return wrapIndices(indices);
    }
    

    double GetInterpolated(const std::vector<double> &x)
    {
        double interpvalue = 0;
        int tempindex;
        for (size_t i = 0 ; i < pow(2,dimension_) ; ++i)
        {
            
            std::vector<double> shiftedvector = x;
            tempindex = i;
                
            double accumulatedweight = 1;
            for(size_t j = 0; j < dimension_ ; ++j)
            {
                double spacing = (GetUpper(j) - GetLower(j)) / GetNumPoints(j);
                shiftedvector[j] += ((tempindex%2)-0.5)*spacing;
                tempindex = tempindex/2;
            }
            
            std::vector<int> shiftedindices = GetIndices(shiftedvector);
            std::vector<double> shiftedcenters = GetCoordinates(shiftedindices);
            
            for(size_t j = 0; j < dimension_ ; ++j)
            {   
                double spacing = (GetUpper(j) - GetLower(j)) / GetNumPoints(j);
                //Handle Edges
                if(shiftedindices[j] == -1)
                    {
                    accumulatedweight *= ((std::abs(x[j]-GetCoordinates(GetIndices(x))[j])/spacing));
                    shiftedvector[j] += spacing/2;
                    }
                else if (shiftedindices[j] == GetNumPoints(j))
                    {
                    accumulatedweight *= ((std::abs(x[j]-GetCoordinates(GetIndices(x))[j])/spacing));
                    shiftedvector[j] -= spacing/2;
                    }
                else
                    {
                    // Handle Periodicity
                    if(std::abs(x[j]-shiftedcenters[j]) > spacing)
                        accumulatedweight *= (1-(std::abs(std::abs(x[j]-shiftedcenters[j]) - (GetUpper(j) - GetLower(j)))/spacing));
                    else
                        accumulatedweight *= (1-(std::abs(x[j]-shiftedcenters[j])/spacing));
                    }
            }           
            interpvalue += accumulatedweight*at(shiftedvector);
        }
        return interpvalue;
    }


    int GetIndex(double x) const
    {
        if (dimension_ != 1) {
            throw std::invalid_argument("1d Grid index can only be accessed for "
                   "1d-Grids can be accessed with a.");
        }
        return GetIndices({x}).at(0);
    }


    std::vector<double> GetCoordinates(const std::vector<int> &indices)
    {
        if (indices.size() != dimension_) {
            throw std::invalid_argument(
                    "Grid indices specified for GetCoordinates() do not have "
                    "the same dimensionality as the grid.");
        }

        std::vector<double> v(dimension_);

        for (size_t i = 0; i < dimension_; ++i) {
            double spacing = (GetUpper(i) - GetLower(i)) / GetNumPoints(i);
            v.at(i) = GetLower(i) + (indices[i] + 0.5)*spacing;
        }

        return v;
    }


    double GetCoordinate(int index)
    {
        return GetCoordinates({index}).at(0);
    }


    const T& at(const std::vector<int> &indices) const
    {
        // Check that indices are in bound.
        if (indices.size() != GetDimension()) {
            throw std::invalid_argument("Dimension of indices does not match "
                    "dimension of the grid.");
        }

        return data_.at(mapTo1d(indices));
    }


    T& at(const std::vector<int> &indices)
    {
        return const_cast<T&>(static_cast<const GridBase<T>* >(this)->at(indices));
    }


    template<typename R>
    const T& at(std::initializer_list<R>&& x) const
    {
        return at(static_cast<std::vector<R> >(x));
    }


    template<typename R>
    T& at(std::initializer_list<R>&& x)
    {
        return at(static_cast<std::vector<R> >(x));
    }


    const T& at(int index) const
    {
        if (dimension_ != 1) {
            throw std::invalid_argument("Only 1d-Grids can be accessed with a "
                    "single integer as the index.");
        }
        return at({index});
    }


    T& at(int index)
    {
        return const_cast<T&>(static_cast<const GridBase<T>* >(this)->at(index));
    }


    const T& at(const std::vector<double> &x) const
    {
        return at(GetIndices(x));
    }


    T& at(const std::vector<double> &x)
    {
        return at(GetIndices(x));
    }


    const T& at(double x) const
    {
        if (dimension_ != 1) {
            throw std::invalid_argument("Only 1d-Grids can be accessed with a "
                    "single float as the specified point.");
        }
        return at({x});
    }


    T& at(double x)
    {
        return const_cast<T&>(static_cast<const GridBase<T>* >(this)->at(x));
    }


    const T& operator[](const std::vector<int> &indices) const
    {
        return at(indices);
    }


    T& operator[](const std::vector<int> &indices)
    {
        return at(indices);
    }


    template<typename R>
    const T& operator[](std::initializer_list<R>&& x) const
    {
        return at(static_cast<std::vector<R> >(x));
    }


    template<typename R>
    T& operator[](std::initializer_list<R>&& x)
    {
        return at(static_cast<std::vector<R> >(x));
    }


    const T& operator[](int index) const
    {
        return at(index);
    }


    T& operator[](int index)
    {
        return at(index);
    }


    const T& operator[](const std::vector<double> &x) const
    {
        return at(x);
    }


    T& operator[](const std::vector<double> &x)
    {
        return at(x);
    }


    const T& operator[](double x) const
    {
        return at(x);
    }


    T& operator[](double x)
    {
        return at(x);
    }
};

} // End namespace SSAGES