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>::numPoints_[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>::numPoints_[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>::numPoints_[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::CharReaderBuilder rbuilder;
        Json::CharReader* reader = rbuilder.newCharReader();
 
        reader->parse(JsonSchema::grid.c_str(),
                      JsonSchema::grid.c_str() + JsonSchema::grid.size(),
                      &schema, nullptr);
        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_->numPoints_[i]) {
                    indices_.at(i) = 0;
                    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 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_->numPoints_[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>::numPoints_[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>::numPoints_[i] - 1;
        }
 
        iterator it(indices, this);
        return ++it;
    }
 
 
    void WriteToFile(const std::string& filename)
    {
        std::ofstream output(filename.c_str(), std::ofstream::out);
 
        // Print out header.
        output << "#! type grid\n";
        output << "#! dim  " << GridBase<T>::dimension_ << "\n"; 
        output << "#! count "; 
        for(auto& c : this->GetNumPoints()) 
            output << c << " "; 
        output << "\n";
        output << "#! lower "; 
        for(auto& l : this->GetLower())
            output << l << " "; 
        output << "\n"; 
        output << "#! upper "; 
        for(auto& u : this->GetUpper())
            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);
    if (!file)
        throw BuildException ({"Attempt to load grid " + filename + " unsuccessful."});
 
        std::string line, buff; 
 
        // Skip type for now. 
        std::getline(file, line); 
 
        // Get dimension.
        {
            size_t 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)
            {
                double spacing = 0;
                if(GridBase<T>::isPeriodic_[i])
                {
                    spacing = (GridBase<T>::edges_.second[i] - GridBase<T>::edges_.first[i]) / (GridBase<T>::numPoints_[i]);
                }               
                if(std::abs(lower-spacing/2 - lowers[i]) > 1e-8)
                    throw std::invalid_argument(filename + 
                    ": Expected lower " + std::to_string(lowers[i]+spacing/2) + 
                    " 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)
            {
                double spacing = 0;
                if(GridBase<T>::isPeriodic_[i])
                {
                    spacing = (GridBase<T>::edges_.second[i] - GridBase<T>::edges_.first[i]) / (GridBase<T>::numPoints_[i]);
                }
                if(std::abs((upper+spacing/2) - uppers[i]) > 1e-8)
                    throw std::invalid_argument(filename + 
                    ": Expected upper " + std::to_string(uppers[i]-spacing/2) + 
                    " 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(size_t i = 0; i < GridBase<T>::dimension_; ++i) 
                iss >> buff;
            
            iss >> d;
        }
    }
};
 
} // End namespace SSAGES