Grid.h

#pragma once

#include <exception>
#include <fstream>
#include <istream>

#include "Drivers/DriverException.h"
#include "GridBase.h"
#include "schema.h"
#include "Validator/ObjectRequirement.h"

namespace SSAGES
{


template<typename T>
class Grid : public GridBase<T>
{
private:

    size_t mapTo1d(const std::vector<int> &indices) const override
    {
        // Check if an index is out of bounds
        for (size_t i=0; i < GridBase<T>::GetDimension(); ++i) {
            int index = indices.at(i);
            int numpoints = GridBase<T>::GetNumPoints(i);
            if ( index < 0 || index >= numpoints )
            {
                throw std::out_of_range("Grid index out of range.");
            }
        }

        size_t idx = 0;
        size_t fac = 1;
        for (size_t i=0; i < GridBase<T>::GetDimension(); ++i) {
            idx += indices.at(i) * fac;
            fac *= GridBase<T>::GetNumPoints(i);
        }
        return idx;
    }

public:

    Grid(std::vector<int> numPoints,
              std::vector<double> lower,
              std::vector<double> upper,
              std::vector<bool> isPeriodic)
      : GridBase<T>(numPoints, lower, upper, isPeriodic)
    {
        size_t data_size = 1;
        for (size_t d = 0; d < GridBase<T>::GetDimension(); ++d) {
            size_t storage_size = GridBase<T>::GetNumPoints(d);
            data_size *= storage_size;
        }

        GridBase<T>::data_.resize(data_size);
    }


    static Grid<T>* BuildGrid(const Json::Value& json)
    {
        return BuildGrid(json, "#/Grid");
    }


    static Grid<T>* BuildGrid(const Json::Value& json, const std::string& path)
    {
        Json::ObjectRequirement validator;
        Json::Value schema;
        Json::Reader reader;

        reader.parse(JsonSchema::grid, schema);
        validator.Parse(schema, path);

        // Validate inputs.
        validator.Validate(json, path);
        if (validator.HasErrors()) {
            throw BuildException(validator.GetErrors());
        }

        // Read in Lower Grid edges.
        std::vector<double> lower;
        for (auto &lv : json["lower"]) {
           lower.push_back(lv.asDouble());
        }

        size_t dimension = lower.size();

        // Read in upper grid edges.
        std::vector<double> upper;
        for (auto &uv : json["upper"]) {
            upper.push_back(uv.asDouble());
        }

        if (upper.size() != dimension) {
            throw BuildException({"Number of upper values does not match "
                                  "number of lower values!"});
        }

        // Read in number of points.
        std::vector<int> number_points;
        for (auto &np : json["number_points"]) {
            number_points.push_back(np.asInt());
        }

        if (number_points.size() != dimension) {
            throw BuildException({"Arrays \"lower\" and \"number_points\" do "
                                  "not have the same size!"});
        }

        // Read in periodicity.
        std::vector<bool> isPeriodic;
        for (auto &periodic : json["periodic"]) {
            isPeriodic.push_back(periodic.asBool());
        }

        if (isPeriodic.size() == 0) {
            isPeriodic = std::vector<bool>(dimension, false);
        } else if (isPeriodic.size() != dimension) {
            throw BuildException({"Arrays \"lower\" and \"periodic\" do not "
                                  "have the same size!"});
        }

        // Construct the grid.
        Grid<T>* grid = new Grid(number_points, lower, upper, isPeriodic);

        return grid;
    }


    template<typename R>
    class GridIterator {
    public:
        typedef GridIterator self_type;

        typedef int difference_type;

        typedef R value_type;

        typedef R* pointer;

        typedef R& reference;

        typedef std::bidirectional_iterator_tag iterator_category;

        GridIterator() = default;


        GridIterator(const std::vector<int> &indices, Grid<T> *grid)
            : indices_(indices), grid_(grid)
        {
        }


        GridIterator(const std::vector<int> &indices, const Grid<T> *grid)
            : indices_(indices), grid_(grid)
        {
        }


        GridIterator(const GridIterator &other)
            : indices_(other.indices_), grid_(other.grid_)
        {
        }


        reference operator*() { return grid_->at(indices_); }


        self_type &operator++()
        {
            indices_.at(0) += 1;
            for (size_t i = 0; i < grid_->GetDimension() - 1; ++i) {
                if (indices_.at(i) >= grid_->GetNumPoints(i)) {
                    indices_.at(i) = 0;
                    indices_.at(i+1) += 1;
                }
            }

            return *this;
        }


        self_type operator++(int)
        {
            GridIterator it(*this);
            ++(*this);
            return it;
        }

        /*
         * \param shift Vector of shifts in each dimension.
         * \return Reference to itself.
         *
         * This operator shifts the current position of the iterator by the
         * given amount in each dimension.
         *
         * Example:
         *
         * \code{.cpp}
         * if += {1,1,1};
         * \endcode
         *
         * In this example the current position of the iterator is shifted
         * diagonally.
         */
        self_type &operator+=(std::vector<int> shift)
        {
            if (shift.size() != grid_->GetDimension()) {
                throw std::invalid_argument("Vector to shift iterator does not "
                                            "match grid dimension.");
            }

            for (size_t i = 0; i < grid_->GetDimension(); ++i) {
                indices_.at(i) += shift.at(i);
            }

            return *this;
        }


        const self_type operator+(std::vector<int> shift)
        {
            return GridIterator(*this) += shift;
        }


        self_type &operator--()
        {
            indices_.at(0) -= 1;
            for (size_t i = 0; i < grid_->GetDimension() - 1; ++i) {
                if (indices_.at(i) < 0) {
                    indices_.at(i) = grid_->GetNumPoints(i)-1;
                    indices_.at(i+1) -= 1;
                }
            }

            return *this;
        }


        self_type operator--(int)
        {
            GridIterator it(*this);
            --(*this);
            return it;
        }


        self_type &operator-=(std::vector<int> shift)
        {
            if (shift.size() != grid_->GetDimension()) {
                throw std::invalid_argument("Vector to shift iterator does not "
                                            "match histogram dimension.");
            }

            for (size_t i = 0; i < grid_->GetDimension(); ++i) {
                indices_.at(i) -= shift.at(i);
            }

            return *this;
        }


        const self_type operator-(std::vector<int> shift)
        {
            return GridIterator(*this) -= shift;
        }


        bool operator==(const self_type &rhs) const
        {
            return indices_ == rhs.indices_ && grid_ == rhs.grid_;
        }


        bool operator!=(const self_type &rhs) const
        {
            return !( (*this) == rhs );
        }


        std::vector<int> &indices()
        {
            return indices_;
        }


        int &index(size_t d)
        {
            return indices()[d];
        }


        std::vector<double> coordinates() const
        {
            return grid_->GetCoordinates(indices_);
        }


        double coordinate(size_t d) const
        {
            return coordinates()[d];
        }

    private:
        std::vector<int> indices_;

        Grid<T> *grid_;
    };

    typedef GridIterator<T> iterator;

    typedef GridIterator<const T> const_iterator;


    iterator begin()
    {
        std::vector<int> indices(GridBase<T>::GetDimension(), 0);

        return iterator(indices, this);
    }


    iterator end()
    {
        std::vector<int> indices(GridBase<T>::GetDimension());
        for (size_t i = 0; i < indices.size(); ++i) {
            indices.at(i) = GridBase<T>::GetNumPoints(i) - 1;
        }

        iterator it(indices, this);
        return ++it;
    }


    const_iterator begin() const
    {
        std::vector<int> indices(GridBase<T>::GetDimension(), 0);

        return iterator(indices, this);
    }


    const_iterator end() const
    {
        std::vector<int> indices(GridBase<T>::GetDimension());
        for (size_t i = 0; i < indices.size(); ++i) {
            indices.at(i) = GridBase<T>::GetNumPoints(i) - 1;
        }

        iterator it(indices, this);
        return ++it;
    }


    void WriteToFile(const std::string& file)
    {
        std::ofstream output(file.c_str(), std::ofstream::out);

        // Print out header.
        output << "#! type grid\n";
        output << "#! dim  " << GridBase<T>::dimension_ << "\n"; 
        output << "#! count "; 
        for(auto& c : GridBase<T>::numPoints_) 
            output << c << " "; 
        output << "\n";
        output << "#! lower "; 
        for(auto& l : GridBase<T>::edges_.first)
            output << l << " "; 
        output << "\n"; 
        output << "#! upper "; 
        for(auto& u : GridBase<T>::edges_.second)
            output << u << " "; 
        output << "\n"; 
        output << "#! periodic "; 
        for(auto p : GridBase<T>::isPeriodic_)
            output << p << " "; 
        output << "\n";

        for(auto it = this->begin(); it != this->end(); ++it)
        {
            auto coords = it.coordinates(); 
            for(auto& c : coords)
                output << std::setprecision(8) << std::fixed << c << " "; 
            output.unsetf(std::ios_base::fixed);
            output << std::setprecision(16) << *it << "\n";
        }

        output.close();
    }


    void LoadFromFile(const std::string& filename)
    {
        std::ifstream file(filename);
        std::string line, buff; 

        // Skip type for now. 
        std::getline(file, line); 

        // Get dimension.
        {
            int dim = 0;
            std::getline(file, line);
            std::istringstream iss(line);
            iss >> buff >> buff >> dim; 
            if(dim != GridBase<T>::dimension_)
                throw std::invalid_argument(filename + 
                    ": Expected dimension " + std::to_string(GridBase<T>::dimension_) + 
                    " but got " + std::to_string(dim) + " instead.");
        }

        // Get size. 
        {
            std::getline(file, line); 
            std::istringstream iss(line);
            int count = 0;
            auto counts = GridBase<T>::numPoints_;
            iss >> buff >> buff;
            for(int i = 0; iss >> count; ++i)
            {
                if(count != counts[i])
                    throw std::invalid_argument(filename + 
                    ": Expected count " + std::to_string(counts[i]) + 
                    " on dimension " + std::to_string(i) + " but got " + 
                    std::to_string(count) + " instead.");
            }
        }

        // Get lower bounds. 
        {
            std::getline(file, line); 
            std::istringstream iss(line);
            double lower = 0; 
            auto lowers = GridBase<T>::edges_.first; 
            iss >> buff >> buff;
            for(int i = 0; iss >> lower; ++i)
            {
                if(std::abs(lower - lowers[i]) > 1e-8)
                    throw std::invalid_argument(filename + 
                    ": Expected lower " + std::to_string(lowers[i]) + 
                    " on dimension " + std::to_string(i) + " but got " + 
                    std::to_string(lower) + " instead.");
            }
        }

        // Get upper bounds. 
        {
            std::getline(file, line); 
            std::istringstream iss(line);
            double upper = 0; 
            auto uppers = GridBase<T>::edges_.second; 
            iss >> buff >> buff;
            for(int i = 0; iss >> upper; ++i)
            {
                if(std::abs(upper - uppers[i]) > 1e-8)
                    throw std::invalid_argument(filename + 
                    ": Expected upper " + std::to_string(uppers[i]) + 
                    " on dimension " + std::to_string(i) + " but got " + 
                    std::to_string(upper) + " instead.");
            }
        }

        // Get periodic. 
        {
            std::getline(file, line); 
            std::istringstream iss(line);
            bool periodic; 
            auto periodics = GridBase<T>::isPeriodic_; 
            iss >> buff >> buff;
            for(int i = 0; iss >> periodic; ++i)
            {
                if(periodic != periodics[i])
                    throw std::invalid_argument(filename + 
                    ": Expected periodic " + std::to_string(periodics[i]) + 
                    " on dimension " + std::to_string(i) + " but got " + 
                    std::to_string(periodic) + " instead.");
            }
        }

        // Finally load data.
        for(auto& d : GridBase<T>::data_)
        {
            std::getline(file, line);
            std::istringstream iss(line);
            
            // Skip coordinates.
            for(int i = 0; i < GridBase<T>::dimension_; ++i) 
                iss >> buff;
            
            iss >> d;
        }
    }
};

} // End namespace SSAGES