CompoundOperators.h

/*
 * ########################################################################
 * The contents of this file is free and unencumbered software released into the
 * public domain. For more information, please refer to <http://unlicense.org/>
 * ########################################################################
 */
 
#ifndef BSPLINE_OPERATORS_COMPOUNDOPERATORS_H
#define BSPLINE_OPERATORS_COMPOUNDOPERATORS_H
 
#include <bspline/operators/GenericOperators.h>
 
namespace bspline::operators {
template <typename O1, typename O2,
          std::enable_if_t<are_operators_v<O1, O2>, bool> = true>
class OperatorProduct final : public Operator {
 private:
  O1 _o1;
  O2 _o2;
 
 public:
  OperatorProduct(O1 o1, O2 o2) : _o1(std::move(o1)), _o2(std::move(o2)){};
 
  static constexpr size_t outputOrder(size_t inputOrder) {
    return O1::outputOrder(O2::outputOrder(inputOrder));
  }
 
  template <typename T, size_t size>
  std::array<T, outputOrder(size - 1) + 1> transform(
      const std::array<T, size> &input, const support::Grid<T> &grid,
      size_t intervalIndex) const {
    return _o1.transform(_o2.transform(input, grid, intervalIndex), grid,
                         intervalIndex);
  }
};
 
template <typename O1, typename O2,
          std::enable_if_t<are_operators_v<O1, O2>, bool> = true>
OperatorProduct<O1, O2> operator*(O1 &&o1, O2 &&o2) {
  return OperatorProduct(std::forward<O1>(o1), std::forward<O2>(o2));
}
 
// ######################### OperatorProduct #############################
// #######################################################################
 
// #######################################################################
// ########################### OperatorSum ###############################
 
enum class AdditionOperation {
  ADDITION,
  SUBTRACTION,
};
 
template <AdditionOperation operation>
inline constexpr bool is_valid_operation_v = (operation ==
                                              AdditionOperation::ADDITION) ||
                                             (operation ==
                                              AdditionOperation::SUBTRACTION);
 
template <
    typename O1, typename O2, AdditionOperation operation,
    std::enable_if_t<are_operators_v<O1, O2> && is_valid_operation_v<operation>,
                     bool> = true>
class OperatorSum final : public Operator {
 private:
  O1 _o1;
  O2 _o2;
 
  template <typename T, size_t sizea, size_t sizeb>
  static std::array<T, std::max(sizea, sizeb)> &add(std::array<T, sizea> &a,
                                                    std::array<T, sizeb> &b) {
    if constexpr (sizeb > sizea) {
      return add(b, a);
    } else {
      for (size_t i = 0; i < sizeb; i++) {
        a[i] += b[i];
      }
      return a;
    }
  }
 
 public:
  OperatorSum(O1 o1, O2 o2) : _o1(std::move(o1)), _o2(std::move(o2)){};
 
  static constexpr size_t outputOrder(size_t inputOrder) {
    return std::max(O1::outputOrder(inputOrder), O2::outputOrder(inputOrder));
  }
 
  template <typename T, size_t size>
  std::array<T, outputOrder(size - 1) + 1> transform(
      const std::array<T, size> &input, const support::Grid<T> &grid,
      size_t intervalIndex) const {
    auto a = _o1.transform(input, grid, intervalIndex);
    auto b = _o2.transform(input, grid, intervalIndex);
 
    // Negate b if subtraction is requested.
    if constexpr (operation == AdditionOperation::SUBTRACTION) {
      for (T &el : b) {
        el *= static_cast<T>(-1);
      }
    }
 
    return add(a, b);
  }
};
 
template <typename O1, typename O2,
          std::enable_if_t<are_operators_v<O1, O2>, bool> = true>
auto operator+(O1 &&o1, O2 &&o2) {
  return OperatorSum<O1, O2, AdditionOperation::ADDITION>(std::forward<O1>(o1),
                                                          std::forward<O2>(o2));
}
 
template <typename O1, typename O2,
          std::enable_if_t<are_operators_v<O1, O2>, bool> = true>
auto operator-(O1 &&o1, O2 &&o2) {
  return OperatorSum<O1, O2, AdditionOperation::SUBTRACTION>(
      std::forward<O1>(o1), std::forward<O2>(o2));
}
 
}  // namespace bspline::operators
#endif  // BSPLINE_OPERATORS_COMPOUNDOPERATORS_H