size_bin_mapping.hpp

Go to the documentation of this file.

#ifndef FIMS_POPULATION_DYNAMICS_SIZE_BIN_MAPPING_HPP
#define FIMS_POPULATION_DYNAMICS_SIZE_BIN_MAPPING_HPP
 
#include <algorithm>
#include <cstddef>
#include <stdexcept>
 
#include "../size_grid.hpp"
#include "../../../common/fims_vector.hpp"
 
namespace fims_popdy {

struct SizeBinMapping {
  static double OverlapWidth(double left_a, double right_a, double left_b,
                             double right_b) {
    return (std::max)(
        0.0, (std::min)(right_a, right_b) - (std::max)(left_a, left_b));
  }

  static bool HasStrictlyIncreasingEdges(const fims::Vector<double>& edges) {
    if (edges.size() < 2) {
      return false;
    }
 
    for (std::size_t i = 0; i + 1 < edges.size(); ++i) {
      if (!(edges[i + 1] > edges[i])) {
        return false;
      }
    }
 
    return true;
  }

  static fims::Vector<fims::Vector<double>> BuildRebinWeights(
      const fims::Vector<double>& source_edges,
      const fims::Vector<double>& destination_edges) {
    if (!HasStrictlyIncreasingEdges(source_edges) ||
        !HasStrictlyIncreasingEdges(destination_edges)) {
      throw std::runtime_error(
          "SizeBinMapping requires strictly increasing source and "
          "destination edges");
    }
 
    if (!DestinationCoversSourceRange(source_edges, destination_edges)) {
      throw std::runtime_error(
          "SizeBinMapping requires destination edges to cover the full "
          "source-bin range");
    }
 
    const std::size_t n_source_bins = source_edges.size() - 1;
    const std::size_t n_destination_bins = destination_edges.size() - 1;
 
    fims::Vector<fims::Vector<double>> weights;
    weights.resize(n_destination_bins);
 
    for (std::size_t i = 0; i < n_destination_bins; ++i) {
      weights[i].resize(n_source_bins);
 
      for (std::size_t j = 0; j < n_source_bins; ++j) {
        const double source_left = source_edges[j];
        const double source_right = source_edges[j + 1];
        const double destination_left = destination_edges[i];
        const double destination_right = destination_edges[i + 1];
        const double source_width = source_right - source_left;
        const double overlap = OverlapWidth(
            source_left, source_right, destination_left, destination_right);
 
        weights[i][j] = source_width > 0.0 ? overlap / source_width : 0.0;
      }
    }
 
    return weights;
  }

  static bool DestinationCoversSourceRange(
      const fims::Vector<double>& source_edges,
      const fims::Vector<double>& destination_edges) {
    if (!HasStrictlyIncreasingEdges(source_edges) ||
        !HasStrictlyIncreasingEdges(destination_edges)) {
      return false;
    }
 
    return destination_edges[0] <= source_edges[0] &&
           destination_edges[destination_edges.size() - 1] >=
               source_edges[source_edges.size() - 1];
  }

  static fims::Vector<double> ExpandDestinationEdgesToCoverSourceRange(
      const fims::Vector<double>& source_edges,
      const fims::Vector<double>& destination_edges) {
    if (!HasStrictlyIncreasingEdges(source_edges) ||
        !HasStrictlyIncreasingEdges(destination_edges)) {
      throw std::runtime_error(
          "SizeBinMapping requires strictly increasing source and "
          "destination edges");
    }
 
    fims::Vector<double> expanded_edges = destination_edges;
 
    if (expanded_edges[0] > source_edges[0]) {
      expanded_edges[0] = source_edges[0];
    }
 
    const std::size_t last_destination_edge = expanded_edges.size() - 1;
    const std::size_t last_source_edge = source_edges.size() - 1;
 
    if (expanded_edges[last_destination_edge] <
        source_edges[last_source_edge]) {
      expanded_edges[last_destination_edge] = source_edges[last_source_edge];
    }
 
    return expanded_edges;
  }

  template <typename Type>
  static fims::Vector<Type> ApplyRebinWeights(
      const fims::Vector<fims::Vector<double>>& weights,
      const fims::Vector<Type>& source_mass) {
    fims::Vector<Type> destination_mass;
    destination_mass.resize(weights.size());
 
    for (std::size_t i = 0; i < weights.size(); ++i) {
      Type total = Type(0.0);
 
      if (weights[i].size() != source_mass.size()) {
        throw std::runtime_error(
            "SizeBinMapping weight rows must match source mass length");
      }
 
      for (std::size_t j = 0; j < source_mass.size(); ++j) {
        total += source_mass[j] * weights[i][j];
      }
 
      destination_mass[i] = total;
    }
 
    return destination_mass;
  }
};
 
}  // namespace fims_popdy
 
#endif /* FIMS_POPULATION_DYNAMICS_SIZE_BIN_MAPPING_HPP */