Histogram.h

#pragma once

#include <exception>

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

namespace SSAGES
{


template<typename T>
class Histogram : 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);
            bool periodic = GridBase<T>::GetPeriodic(i);
            int numpoints = GridBase<T>::GetNumPoints(i);
            if ( (periodic && (index < 0 || index >= numpoints)) ||
                 (periodic && (index < -1 || index > numpoints)) )
            {
                throw std::out_of_range("Bin index out of range.");
            }
        }

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

public:

    Histogram(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);
            if (!GridBase<T>::GetPeriodic(d)) { storage_size += 2; }
            data_size *= storage_size;
        }

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


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


    static Histogram<T>* BuildHistogram(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.
        Histogram<T>* hist = new Histogram(number_points, lower, upper, isPeriodic);

        return hist;
    }


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

        typedef int difference_type;

        typedef R value_type;

        typedef R* pointer;

        typedef R& reference;

        typedef std::bidirectional_iterator_tag iterator_category;

        HistIterator() = default;


        HistIterator(const std::vector<int> &indices, Histogram<T> *hist)
            : indices_(indices), hist_(hist)
        {
        }


        HistIterator(const HistIterator &other)
            : indices_(other.indices_), hist_(other.hist_)
        {
        }


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


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

                    indices_.at(i) = 0;
                    indices_.at(i+1) += 1;
                } else if (!hist_->GetPeriodic(i) &&
                           indices_.at(i) > hist_->GetNumPoints(i)) {

                    indices_.at(i) = -1;
                    indices_.at(i+1) += 1;
                }
            }

            return *this;
        }


        self_type operator++(int)
        {
            HistIterator 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:
         *
         * if += {1,1,1};
         *
         * In this example the current position of the iterator is shifted
         * diagonally.
         */
        self_type &operator+=(std::vector<int> shift)
        {
            if (shift.size() != hist_->GetDimension()) {
                throw std::invalid_argument("Vector to shift iterator does not "
                                            "match histogram dimension.");
            }

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

            return *this;
        }


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


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

            return *this;
        }


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


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

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

            return *this;
        }


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


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


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


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


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


        bool isUnderOverflowBin() const
        {
            for (size_t i = 0; i < indices_.size(); ++i) {
                if (indices_.at(i) == -1 || indices_.at(i) == hist_->GetNumPoints(i)) {
                    return true;
                }
            }
            return false;
        }


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


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

    private:
        std::vector<int> indices_;

        Histogram<T> *hist_;
    };

    typedef HistIterator<T> iterator;

    typedef HistIterator<const T> const_iterator;


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

        return iterator(indices, this);
    }


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

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


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

        return iterator(indices, this);
    }


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

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

} // End namespace SSAGES