forage_values.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2020 W. Traylor <wolfgang.traylor@senckenberg.de>
//
// SPDX-License-Identifier: LGPL-3.0-or-later
 
#ifndef FAUNA_FORAGE_VALUES_H
#define FAUNA_FORAGE_VALUES_H
 
#include <array>
#include <cmath>
#include <numeric>
#include <stdexcept>
#include "Fauna/average.h"
#include "Fauna/forage_types.h"
 
namespace Fauna {
class HerbivoreInterface;  // for ForageDistribution
 
enum class ForageValueTag {
  PositiveAndZero,
  ZeroToOne
};
 
 
template <ForageValueTag tag>
class ForageValues {
 public:
 
  ForageValues(const double init_value = 0.0) { set(init_value); }
 
 
  ForageValues<tag> divide_safely(const ForageValues<tag>& divisor,
                                  const double na_value) const {
    ForageValues<tag> result(*this);
    for (int ft = 0; ft < array.size(); ft++) {
      const double d = divisor.array[ft];
      if (d != 0.0)
        result.set((ForageType)ft, array[ft] / d);  // normal
      else
        result.set((ForageType)ft, na_value);  // division by zero
    }
    return result;
  }
 
 
  double get(const ForageType ft) const {
    if (ft == ForageType::Inedible)
      throw std::invalid_argument(
          "Fauna::ForageValues<>::get() "
          "The forage type `ForageType::Inedible` is not allowed.");
    assert((int)ft < array.size());
    assert((int)ft >= 0);
    return array[(int)ft];
  }
 
 
  ForageValues<tag>& merge(const ForageValues<tag>& other,
                           const double this_weight = 1.0,
                           const double other_weight = 1.0) {
    for (int ft = 0; ft < array.size(); ft++)
      set((ForageType)ft,
          average(array[ft], other.array[ft], this_weight, other_weight));
    return *this;
  }
 
 
  ForageValues<tag>& max(const ForageValues<tag>& other) {
    if (&other == this) return *this;
    for (int ft = 0; ft < array.size(); ft++)
      set((ForageType)ft, std::max(array[ft], other.array[ft]));
    return *this;
  }
 
 
  ForageValues<tag>& min(const ForageValues<tag>& other) {
    if (&other == this) return *this;
    for (int ft = 0; ft < array.size(); ft++)
      set((ForageType)ft, std::min(array[ft], other.array[ft]));
    return *this;
  }
 
  double operator[](const ForageType ft) const { return get(ft); }
 
  double& operator[](const ForageType ft) {
    if (ft == ForageType::Inedible)
      throw std::invalid_argument(
          "Fauna::ForageValues<>::get() "
          "The forage type `ForageType::Inedible` is not allowed.");
    assert((int)ft < array.size());
    assert((int)ft >= 0);
    return array[(int)ft];
  }
 
 
  void set(const ForageType forage_type, double value) {
    check_value(value);
    if (forage_type == ForageType::Inedible)
      throw std::invalid_argument((std::string)"ForageValues<> "
                            "Forage type `ForageType::Inedible` is not allowed."+
                            " ("+get_forage_type_name(forage_type)+")");
 
    // Change the value.
    assert((int)forage_type < array.size());
    assert((int)forage_type >= 0);
    array[(int)forage_type] = value;
  }
 
 
  void set(double value) {
    check_value(value);
    array.fill(value);
  }
 
  double sum() const {
    return std::accumulate(array.begin(), array.end(), 0.0);
  }
 
  ForageValues<tag>& operator+=(const double rhs) {
    for (int ft = 0; ft < array.size(); ft++)
      set((ForageType)ft, array[ft] + rhs);
    return *this;
  }
  ForageValues<tag>& operator-=(const double rhs) {
    for (int ft = 0; ft < array.size(); ft++)
      set((ForageType)ft, array[ft] - rhs);
    return *this;
  }
  ForageValues<tag>& operator*=(const double rhs) {
    for (int ft = 0; ft < array.size(); ft++)
      set((ForageType)ft, array[ft] * rhs);
    return *this;
  }
  ForageValues<tag>& operator/=(const double rhs) {
    if (rhs == 0)
      throw std::domain_error("Fauna::ForageValues<> Division by zero.");
    for (int ft = 0; ft < array.size(); ft++)
      set((ForageType)ft, array[ft] / rhs);
    return *this;
  }
 
  ForageValues<tag> operator+(const double rhs) const {
    ForageValues<tag> result(*this);
    return result.operator+=(rhs);
  }
  ForageValues<tag> operator-(const double rhs) const {
    ForageValues<tag> result(*this);
    return result.operator-=(rhs);
  }
  ForageValues<tag> operator*(const double rhs) const {
    ForageValues<tag> result(*this);
    return result.operator*=(rhs);
  }
  ForageValues<tag> operator/(const double rhs) const {
    ForageValues<tag> result(*this);
    return result.operator/=(rhs);
  }
 
  ForageValues<tag> operator+(const ForageValues<tag>& rhs) const {
    ForageValues<tag> result(*this);
    return result.operator+=(rhs);
  }
  ForageValues<tag> operator-(const ForageValues<tag>& rhs) const {
    ForageValues<tag> result(*this);
    return result.operator-=(rhs);
  }
  ForageValues<tag> operator*(const ForageValues<tag>& rhs) const {
    ForageValues<tag> result(*this);
    return result.operator*=(rhs);
  }
  ForageValues<tag> operator/(const ForageValues<tag>& rhs) const {
    ForageValues<tag> result(*this);
    return result.operator/=(rhs);
  }
 
  ForageValues<tag>& operator+=(const ForageValues<tag>& rhs) {
    for (int ft = 0; ft < array.size(); ft++)
      set((ForageType)ft, array[ft] + rhs.array[ft]);
    return *this;
  }
  ForageValues<tag>& operator-=(const ForageValues<tag>& rhs) {
    for (int ft = 0; ft < array.size(); ft++)
      set((ForageType)ft, array[ft] - rhs.array[ft]);
    return *this;
  }
  ForageValues<tag>& operator*=(const ForageValues<tag>& rhs) {
    for (int ft = 0; ft < array.size(); ft++)
      set((ForageType)ft, array[ft] * rhs.array[ft]);
    return *this;
  }
  ForageValues<tag>& operator/=(const ForageValues<tag>& rhs) {
    for (int ft = 0; ft < array.size(); ft++) {
      if (rhs.array[ft] == 0)
        throw std::domain_error(
            (std::string) "Fauna::ForageValues<> Division by zero." + " (" +
            get_forage_type_name((ForageType)ft) + ")");
      set((ForageType)ft, array[ft] / rhs.array[ft]);
    }
    return *this;
  }
 
  bool operator==(const ForageValues<tag>& rhs) const {
    for (int ft = 0; ft < array.size(); ft++)
      if (array[ft] != rhs.array[ft]) return false;
    return true;
  }
 
  // TODO: Is this function necessary and logical?
  bool operator!=(const ForageValues<tag>& rhs) const {
    for (int ft = 0; ft < array.size(); ft++)
      if (array[ft] == rhs.array[ft]) return false;
    return true;
  }
 
  bool operator<(const ForageValues<tag>& rhs) const {
    for (int ft = 0; ft < array.size(); ft++)
      if (array[ft] >= rhs.array[ft]) return false;
    return true;
  }
  bool operator<=(const ForageValues<tag>& rhs) const {
    for (int ft = 0; ft < array.size(); ft++)
      if (array[ft] > rhs.array[ft]) return false;
    return true;
  }
  bool operator>(const ForageValues<tag>& rhs) const {
    for (int ft = 0; ft < array.size(); ft++)
      if (array[ft] <= rhs.array[ft]) return false;
    return true;
  }
  bool operator>=(const ForageValues<tag>& rhs) const {
    for (int ft = 0; ft < array.size(); ft++)
      if (array[ft] < rhs.array[ft]) return false;
    return true;
  }  // Operator overload
 
 
  constexpr static const double IMPRECISION_TOLERANCE = 1e-3;
 
 private:
  std::array<double, 1> array;
 
 
  void check_value(double& value) const {
    switch (tag) {
      case ForageValueTag::PositiveAndZero:
        if (value < 0.0) {
          // Correct floating point rounding errors.
          if (value >= -IMPRECISION_TOLERANCE)
            value = 0.0;
          else
            throw std::invalid_argument(
                "ForageValues<PositiveAndZero> Value < 0 not allowed. "
                "(value == " +
                std::to_string(value) + ")");
        }
        break;
      case ForageValueTag::ZeroToOne:
        if (value < 0.0 || value > 1.0) {
          // Correct floating point rounding errors.
          if (value < 0.0 && value >= 0.0 - IMPRECISION_TOLERANCE)
            value = 0.0;
          else if (value > 1.0 && value <= 1.0 + IMPRECISION_TOLERANCE)
            value = 1.0;
          else
            throw std::invalid_argument(
                "ForageValues<ZeroToOne> Value is not in interval [0,1]. "
                "(value == " +
                std::to_string(value) + ")");
        }
        break;
      default:
        throw std::logic_error(
            "ForageValues<> "
            "ForageValueTag not implemented.");
    }
    if (std::isnan(value))
      throw std::invalid_argument(
          "ForageValues<> NAN is not allowed as a value.");
    if (std::isinf(value))
      throw std::invalid_argument(
          "ForageValues<> INFINITY is not allowed as a value.");
  }
};
 
typedef ForageValues<ForageValueTag::ZeroToOne> Digestibility;
 
typedef ForageValues<ForageValueTag::PositiveAndZero> ForageEnergy;
 
typedef ForageValues<ForageValueTag::PositiveAndZero> ForageEnergyContent;
 
typedef ForageValues<ForageValueTag::ZeroToOne> ForageFraction;
 
typedef ForageValues<ForageValueTag::PositiveAndZero> ForageMass;
 
typedef std::vector<std::pair<HerbivoreInterface*, ForageMass>>
    ForageDistribution;
 
 
inline ForageValues<ForageValueTag::PositiveAndZero> operator*(
    const double lhs, const ForageFraction& rhs) {
  ForageValues<ForageValueTag::PositiveAndZero> result;
  for (const auto ft : FORAGE_TYPES) result.set(ft, rhs[ft] * lhs);
  return result;
}
 
inline ForageValues<ForageValueTag::PositiveAndZero> operator*(
    const ForageFraction& lhs,
    const ForageValues<ForageValueTag::PositiveAndZero>& rhs) {
  ForageValues<ForageValueTag::PositiveAndZero> result;
  for (const auto ft : FORAGE_TYPES) result.set(ft, rhs[ft] * lhs[ft]);
  return result;
}
 
inline ForageValues<ForageValueTag::PositiveAndZero> operator*(
    const ForageValues<ForageValueTag::PositiveAndZero>& lhs,
    const ForageFraction& rhs) {
  return operator*(rhs, lhs);
}
 
// Note that we don’t use the typedef "ForageFraction" as the type of parameter
// "fractions" here because Doxygen (1.8.16) does not recognize it and throws a
// warning.
ForageValues<ForageValueTag::PositiveAndZero> foragefractions_to_foragevalues(
    const ForageValues<ForageValueTag::ZeroToOne>& fractions);
 
 
ForageFraction foragevalues_to_foragefractions(
    const ForageValues<ForageValueTag::PositiveAndZero> values,
    const double tolerance);
 
 
// Note that we don’t use the typedefs "ForageFraction" and "ForageEnergy" as
// the parameter types here because Doxygen (1.8.16) does not recognize it and
// throws a warning.
ForageFraction convert_mj_to_kg_proportionally(
    const ForageValues<ForageValueTag::PositiveAndZero>& mj_per_kg,
    const ForageValues<ForageValueTag::ZeroToOne>& mj_proportions);
 
}  // namespace Fauna
#endif  // FAUNA_FORAGE_VALUES_H