VariantMatrix.hpp

#ifndef SEQUENCE_VARIANT_MATRIX_HPP__
#define SEQUENCE_VARIANT_MATRIX_HPP__

#include <algorithm>
#include <cstdint>
#include <cstddef>
#include <utility>
#include <vector>
#include <limits>
#include <stdexcept>
#include <type_traits>

static_assert(sizeof(std::int8_t) == sizeof(char),
              "sizeof(char) is not 8 bits");

namespace Sequence
{

    class VariantMatrix
    {
      private:
        std::int8_t
        set_max_allele(const std::int8_t max_allele_value)
        {
            if (max_allele_value < 0 && !data.empty())
                {
                    auto itr = std::max_element(data.begin(), data.end());
                    return *itr;
                }
            //special case to allow construction of empty data sets
            //without throwing an exception
            else if (data.empty())

                {
                    return 0;
                }
            return max_allele_value;
        }

      public:
        std::vector<std::int8_t> data;
        std::vector<double> positions;
        std::size_t nsites;
        std::size_t nsam;
        static const std::int8_t mask;
        using value_type = std::int8_t;
        const std::int8_t max_allele;
        template <typename data_input, typename positions_input>
        VariantMatrix(data_input&& data_, positions_input&& positions_,
                      const std::int8_t max_allele_value = -1)
            : data(std::forward<data_input>(data_)),
              positions(std::forward<positions_input>(positions_)),
              nsites(positions.size()),
              nsam((nsites > 0) ? data.size() / positions.size() : 0),
              max_allele{ set_max_allele(max_allele_value) }
        {
            if (max_allele < 0)
                {
                    throw std::invalid_argument("max allele must be >= 0");
                }
            if ((!data.empty() && !positions.empty())
                && data.size() % positions.size() != 0.0)
                {
                    throw std::invalid_argument("incorrect dimensions");
                }
        }

        // Non range-checked access

        std::int8_t& get(const std::size_t site, const std::size_t haplotype);
        const std::int8_t& get(const std::size_t site,
                               const std::size_t haplotype) const;

        // Ranged-checked access after std::vector<T>::at.
        std::int8_t& at(const std::size_t site, const std::size_t haplotype);
        const std::int8_t& at(const std::size_t site,
                              const std::size_t haplotype) const;
    };

} // namespace Sequence

#endif