diff --git a/algorithm/include/gnuradio-4.0/algorithm/dataset/DataSetEstimators.hpp b/algorithm/include/gnuradio-4.0/algorithm/dataset/DataSetEstimators.hpp new file mode 100644 index 00000000..ef393708 --- /dev/null +++ b/algorithm/include/gnuradio-4.0/algorithm/dataset/DataSetEstimators.hpp @@ -0,0 +1,511 @@ +#ifndef DATASETESTIMATORS_HPP +#define DATASETESTIMATORS_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include "DataSetHelper.hpp" + +namespace gr::dataset { + +template +[[nodiscard]] constexpr T inverseDecibel(T x) noexcept { + return gr::math::pow(T(10), x / T(20)); // Inverse decibel => 10^(value / 20) +} + +namespace estimators { +template +[[nodiscard]] constexpr T computeCentreOfMass(const DataSet& ds, std::size_t minIndex = 0UZ, std::size_t maxIndex = 0UZ, std::size_t signalIndex = 0) { + if (maxIndex == 0UZ) { // renormalise default range + maxIndex = ds.axisValues(dim::X).size(); + } + detail::checkRangeIndex(minIndex, maxIndex, ds.axisValues(dim::X).size()); + T com = T(0); + T mass = T(0); + + for (std::size_t i = minIndex; i < maxIndex; i++) { + const T x = getIndexValue(ds, dim::X, i); + const T value = getIndexValue(ds, dim::Y, i, signalIndex); + if (gr::math::isfinite(x) && gr::math::isfinite(value)) { + com += x * value; + mass += value; + } + } + if (gr::value(mass) == gr::value(T(0))) { + return std::numeric_limits::quiet_NaN(); + } + return com / mass; +} + +template> +[[nodiscard]] constexpr TValue computeFWHM(const T& data, std::size_t index) { + using value_t = gr::meta::fundamental_base_value_type_t; // innermost value + if (!(index > 0UZ && index < data.size() - 1UZ)) { + return std::numeric_limits::quiet_NaN(); + } + TValue maxHalf = value_t(0.5) * gr::value(data[index]); + std::size_t lowerLimit; + std::size_t upperLimit; + for (upperLimit = index; upperLimit < data.size() && data[upperLimit] > maxHalf; upperLimit++) { + // done in condition + } + for (lowerLimit = index; lowerLimit >= 0 && data[lowerLimit] > maxHalf; lowerLimit--) { + // done in condition + } + if (upperLimit >= data.size() || lowerLimit < 0) { + return std::numeric_limits::quiet_NaN(); + } + return static_cast(upperLimit - lowerLimit); +} + +template +[[nodiscard]] constexpr TValue computeInterpolatedFWHM(const T& data, std::size_t index) { + using value_t = gr::meta::fundamental_base_value_type_t; // innermost value + if (!(index > 0 && index < data.size() - 1)) { + return std::numeric_limits::quiet_NaN(); + } + TValue maxHalf = value_t(0.5) * data[index]; + std::size_t lowerLimit; + std::size_t upperLimit; + for (upperLimit = index; upperLimit < data.size() && data[upperLimit] > maxHalf; upperLimit++) { + // done in condition + } + for (lowerLimit = index; lowerLimit >= 0 && data[lowerLimit] > maxHalf; lowerLimit--) { + // done in condition + } + if (upperLimit >= data.size() || lowerLimit < 0) { + return std::numeric_limits::quiet_NaN(); + } + TValue lowerRefined = detail::linearInterpolate(value_t(lowerLimit), value_t(lowerLimit + 1), data[lowerLimit], data[lowerLimit + 1], maxHalf); + TValue upperRefined = detail::linearInterpolate(value_t(upperLimit - 1), value_t(upperLimit), data[upperLimit - 1], data[upperLimit], maxHalf); + return upperRefined - lowerRefined; +} + +template> +[[nodiscard]] constexpr int getLocationMinimum(const DataSet& dataSet, std::size_t indexMin = 0UZ, std::size_t indexMax = 0UZ, std::size_t signalIndex = 0UZ) { + if (indexMax == 0UZ) { // renormalise default range + indexMax = dataSet.axisValues(dim::X).size(); + } + detail::checkRangeIndex(indexMin, indexMax, dataSet.axisValues(dim::X).size()); + int locMin = -1; + T minVal = std::numeric_limits::infinity(); + for (std::size_t i = indexMin; i < indexMax; i++) { + if (T actual = getIndexValue(dataSet, dim::Y, i, signalIndex); gr::math::isfinite(actual) && actual < minVal) { + minVal = actual; + locMin = static_cast(i); + } + } + return locMin; +} + +template> +[[nodiscard]] constexpr int getLocationMaximum(const DataSet& dataSet, std::size_t indexMin = 0UZ, std::size_t indexMax = 0UZ, std::size_t signalIndex = 0UZ) { + if (indexMax == 0UZ) { // renormalise default range + indexMax = dataSet.axisValues(dim::X).size(); + } + detail::checkRangeIndex(indexMin, indexMax, dataSet.axisValues(dim::X).size()); + int locMax = -1; + T maxVal = -std::numeric_limits::infinity(); + for (std::size_t i = indexMin; i < indexMax; i++) { + if (T actual = getIndexValue(dataSet, dim::Y, i, signalIndex); gr::math::isfinite(actual) && actual > maxVal) { + maxVal = actual; + locMax = static_cast(i); + } + } + return locMax; +} + +template> +[[nodiscard]] constexpr T getMinimum(const DataSet& dataSet, std::size_t indexMin = 0UZ, std::size_t indexMax = 0UZ, std::size_t signalIndex = 0UZ) { + if (indexMax == 0UZ) { // renormalise default range + indexMax = dataSet.axisValues(dim::X).size(); + } + detail::checkRangeIndex(indexMin, indexMax, dataSet.axisValues(dim::X).size()); + + T val = std::numeric_limits::max(); + bool foundAny = false; + for (std::size_t i = indexMin; i < indexMax; i++) { + if (T actual = getIndexValue(dataSet, dim::Y, i, signalIndex); gr::math::isfinite(actual)) { + foundAny = true; + val = std::min(val, actual); + } + } + return foundAny ? val : std::numeric_limits::quiet_NaN(); +} + +template> +[[nodiscard]] constexpr T getMaximum(const DataSet& dataSet, std::size_t indexMin = 0UZ, std::size_t indexMax = 0UZ, std::size_t signalIndex = 0UZ) { + if (indexMax == 0UZ) { // renormalise default range + indexMax = dataSet.axisValues(dim::X).size(); + } + detail::checkRangeIndex(indexMin, indexMax, dataSet.axisValues(dim::X).size()); + T val = -std::numeric_limits::min(); + bool foundAny = false; + for (std::size_t i = indexMin; i < indexMax; i++) { + if (T actual = getIndexValue(dataSet, dim::Y, i, signalIndex); gr::math::isfinite(actual)) { + foundAny = true; + val = std::max(val, actual); + } + } + return foundAny ? val : std::numeric_limits::quiet_NaN(); +} + +template> +[[nodiscard]] constexpr T getMean(const DataSet& dataSet, std::size_t indexMin = 0UZ, std::size_t indexMax = 0UZ, std::size_t signalIndex = 0UZ) { + if (indexMax == 0UZ) { + indexMax = dataSet.axisValues(dim::X).size(); + } + detail::checkRangeIndex(indexMin, indexMax, dataSet.axisValues(dim::X).size()); + + auto signalRange = dataSet.signalValues(signalIndex) | std::views::drop(indexMin) | std::views::take(indexMax - indexMin); + auto finiteValues = signalRange | std::views::filter(gr::math::isfinite); + + T sum = 0; + T count = 0; + for (const auto& val : finiteValues) { + sum += val; + count += T(1); + } + + return count > T(0) ? sum / count : std::numeric_limits::quiet_NaN(); +} + +template> +[[nodiscard]] constexpr T getMedian(const DataSet& dataSet, std::size_t indexMin = 0UZ, std::size_t indexMax = 0UZ, std::size_t signalIndex = 0UZ) { + if (indexMax == 0UZ) { // renormalise default range + indexMax = dataSet.axisValues(dim::X).size(); + } + std::span values = dataSet.signalValues(signalIndex).subspan(indexMin, indexMax - indexMin); + + if (values.empty()) { + return std::numeric_limits::quiet_NaN(); + } + + std::vector data(values.begin(), values.end()); // Temporary mutable copy for in-place partitioning + auto mid = data.begin() + data.size() / 2; + std::ranges::nth_element(data, mid); + + if ((data.size() & 1UZ) == 0UZ) { + // even-sized data, calculate the mean of the two middle elements + auto midPrev = std::ranges::max_element(data.begin(), mid); + return static_cast(0.5) * (*midPrev + *mid); + } + + return static_cast(*mid); // odd-sized data, return the middle element +} + +template +[[nodiscard]] constexpr T getRange(const DataSet& dataSet, std::size_t indexMin = 0UZ, std::size_t indexMax = 0UZ, std::size_t signalIndex = 0UZ) { + if (indexMax == 0UZ) { // Default range + indexMax = dataSet.axisValues(dim::X).size(); + } + detail::checkRangeIndex(indexMin, indexMax, dataSet.axisValues(dim::X).size()); + + auto signalRange = dataSet.signalValues(signalIndex) | std::views::drop(indexMin) | std::views::take(indexMax - indexMin); + auto finiteValues = signalRange | std::views::filter(gr::math::isfinite); + + auto [minIt, maxIt] = std::ranges::minmax_element(finiteValues); + if (minIt == finiteValues.end() || maxIt == finiteValues.end()) { + return std::numeric_limits::quiet_NaN(); + } + + return *maxIt - *minIt; +} + +template +[[nodiscard]] constexpr T getRms(const DataSet& dataSet, std::size_t indexMin = 0UZ, std::size_t indexMax = 0UZ, std::size_t signalIndex = 0UZ) { + if (indexMax == 0UZ) { + indexMax = dataSet.axisValues(dim::X).size(); + } + detail::checkRangeIndex(indexMin, indexMax, dataSet.axisValues(dim::X).size()); + + auto signalRange = dataSet.signalValues(signalIndex) | std::views::drop(indexMin) | std::views::take(indexMax - indexMin); + auto finiteValues = signalRange | std::views::filter(gr::math::isfinite); + + if (finiteValues.empty()) { + return T(0); + } + + T sum = 0, sum2 = 0; + T count = 0; + for (const auto& val : finiteValues) { + sum += val; + sum2 += val * val; + count += T(1); + } + T mean1 = (sum / count) * (sum / count); + T mean2 = sum2 / count; + return gr::math::sqrt(mean2 > mean1 ? mean2 - mean1 : mean1 - mean2); // abs for safety +} + +template +[[nodiscard]] constexpr T getIntegral(const DataSet& dataSet, std::size_t indexMin = 0UZ, std::size_t indexMax = 0UZ, std::size_t signalIndex = 0UZ) { + if (indexMax == 0UZ) { // renormalise default range + indexMax = dataSet.axisValues(dim::X).size(); + } + detail::checkRangeIndex(indexMin, indexMax, dataSet.axisValues(dim::X).size()); + + // sign if reversed + T sign_ = detail::sign(T(1), int(indexMax) - int(indexMin)); + std::size_t start = std::min(indexMin, indexMax); + std::size_t stop = std::max(indexMin, indexMax); + T integral = 0; + for (std::size_t i = start; i < stop - 1; i++) { // compute integral via triangulation + T x0 = getIndexValue(dataSet, dim::X, i, signalIndex); + T x1 = getIndexValue(dataSet, dim::X, i + 1, signalIndex); + T y0 = getIndexValue(dataSet, dim::Y, i, signalIndex); + T y1 = getIndexValue(dataSet, dim::Y, i + 1, signalIndex); + T localIntegral = (x1 - x0) * T(0.5) * (y0 + y1); + if (!gr::math::isfinite(localIntegral)) { + // skip + } else { + integral += localIntegral; + } + } + return static_cast(sign_ * integral); +} + +template> +[[nodiscard]] constexpr T getEdgeDetect(const DataSet& dataSet, std::size_t indexMin = 0UZ, std::size_t indexMax = 0UZ, std::size_t signalIndex = 0UZ) { + if (indexMax == 0UZ) { // renormalise default range + indexMax = dataSet.axisValues(dim::X).size(); + } + if (dataSet.axisCount() == 0 || dataSet.axisValues(0).empty() || indexMin >= indexMax) { + return std::numeric_limits::quiet_NaN(); + } + detail::checkRangeIndex(indexMin, indexMax, dataSet.axisValues(0).size()); + + T minVal = getMinimum(dataSet, indexMin, indexMax, signalIndex); + T maxVal = getMaximum(dataSet, indexMin, indexMax, signalIndex); + T range = maxVal > minVal ? maxVal - minVal : minVal - maxVal; // abs + if (!gr::math::isfinite(range) || range == TValue(0)) { + return std::numeric_limits::quiet_NaN(); + } + // check if falling or rising + T startVal = getIndexValue(dataSet, dim::Y, indexMin, signalIndex); + T endVal = getIndexValue(dataSet, dim::Y, indexMax - 1, signalIndex); + bool inverted = (startVal > endVal); + T startTime = getIndexValue(dataSet, dim::X, indexMin, signalIndex); + + T halfCross = inverted ? (maxVal - TValue(0.5) * range) : (minVal + TValue(0.5) * range); + for (std::size_t i = indexMin; i < indexMax; i++) { + T actual = getIndexValue(dataSet, dim::Y, i, signalIndex); + if (!gr::math::isfinite(actual)) { + continue; + } + if ((!inverted && actual > halfCross) || (inverted && actual < halfCross)) { + T xNow = getIndexValue(dataSet, dim::X, i, signalIndex); + return xNow - startTime; + } + } + return std::numeric_limits::quiet_NaN(); +} + +template> +[[nodiscard]] constexpr T getDutyCycle(const DataSet& dataSet, std::size_t indexMin = 0UZ, std::size_t indexMax = 0UZ, std::size_t signalIndex = 0UZ) { + if (indexMax == 0UZ) { + indexMax = dataSet.axisValues(dim::X).size(); + } + detail::checkRangeIndex(indexMin, indexMax, dataSet.axisValues(dim::X).size()); + + T minVal = getMinimum(dataSet, indexMin, indexMax, signalIndex); + T maxVal = getMaximum(dataSet, indexMin, indexMax, signalIndex); + + if (!gr::math::isfinite(minVal) || !gr::math::isfinite(maxVal) || minVal == maxVal) { + return std::numeric_limits::quiet_NaN(); + } + + auto signalRange = dataSet.signalValues(signalIndex) | std::views::drop(indexMin) | std::views::take(indexMax - indexMin); + auto finiteValues = signalRange | std::views::filter(gr::math::isfinite); + + T range = maxVal > minVal ? maxVal - minVal : minVal - maxVal; + const T thresholdMin = minVal + TValue(0.45) * range; + const T thresholdMax = minVal + TValue(0.55) * range; + int countLow = 0, countHigh = 0; + for (const auto& val : finiteValues) { + if (val < thresholdMin) { + countLow++; + } else if (val > thresholdMax) { + countHigh++; + } + } + + const int totalCount = countLow + countHigh; + return totalCount > 0 ? static_cast(countHigh) / static_cast(totalCount) : std::numeric_limits::quiet_NaN(); +} + +template> +[[nodiscard]] constexpr T getFrequencyEstimate(const DataSet& dataSet, std::size_t indexMin = 0UZ, std::size_t indexMax = 0UZ, std::size_t signalIndex = 0UZ) { + // getFrequencyEstimate => naive/simple approach counting edges + if (indexMax == 0UZ) { + indexMax = dataSet.axisValues(dim::X).size(); + } + detail::checkRangeIndex(indexMin, indexMax, dataSet.axisValues(dim::X).size()); + + T minVal = getMinimum(dataSet, indexMin, indexMax, signalIndex); + T maxVal = getMaximum(dataSet, indexMin, indexMax, signalIndex); + if (!gr::math::isfinite(minVal) || !gr::math::isfinite(maxVal) || maxVal == minVal) { + return std::numeric_limits::quiet_NaN(); + } + T range = maxVal > minVal ? maxVal - minVal : minVal - maxVal; // abs + T thresholdMin = minVal + TValue(0.45) * range; + T thresholdMax = minVal + TValue(0.55) * range; + + T startRisingEdge = std::numeric_limits::quiet_NaN(); + T startFallingEdge = std::numeric_limits::quiet_NaN(); + T avgPeriod = 0; + int avgPeriodCount = 0; + + // 0 => below threshold, 1 => above threshold + float actualState = 0.f; + for (std::size_t i = indexMin; i < indexMax; i++) { + T actual = getIndexValue(dataSet, dim::Y, i, signalIndex); + if (!gr::math::isfinite(actual)) { + continue; + } + T x = getIndexValue(dataSet, dim::X, i, signalIndex); + + if (actualState < 0.5f) { + // was low, check if we cross thresholdMax => rising edge + if (actual > thresholdMax) { + actualState = 1.0f; + if (gr::math::isfinite(startRisingEdge)) { + T period = x - startRisingEdge; + startRisingEdge = x; + avgPeriod += period; + avgPeriodCount++; + } else { + startRisingEdge = x; + } + } + } else { + // was high, check if we cross thresholdMin => falling edge + if (actual < thresholdMin) { + actualState = 0.0f; + if (gr::math::isfinite(startFallingEdge)) { + T period = x - startFallingEdge; + startFallingEdge = x; + avgPeriod += period; + avgPeriodCount++; + } else { + startFallingEdge = x; + } + } + } + } + if (avgPeriodCount == 0) { + return std::numeric_limits::quiet_NaN(); + } + return static_cast(avgPeriodCount) / avgPeriod; // freq = # of periods / total time +} + +template +[[nodiscard]] constexpr TValue interpolateGaussian(const T& data, std::size_t index) { + using value_t = gr::meta::fundamental_base_value_type_t; + if (!(index > 0 && index < data.size() - 1)) { + return static_cast(index); // fallback + } + const TValue left = data[index - 1]; + const TValue center = data[index]; + const TValue right = data[index + 1]; + + if (!gr::math::isfinite(left) || !gr::math::isfinite(right) || !gr::math::isfinite(center)) { + return static_cast(index); + } + + if (gr::value(left) <= value_t(0) || gr::value(right) <= value_t(0) || gr::value(center) <= value_t(0)) { + return static_cast(index); + } + TValue val = static_cast(index); + TValue numerator = value_t(0.5f) * gr::math::log(right / left); + TValue denominator = gr::math::log((center * center) / (left * right)); + if (gr::value(denominator) == TValue(0)) { + return val; + } + return val + numerator / denominator; +} + +// getLocationMaximumGaussInterpolated => returns X-value +template> +[[nodiscard]] constexpr T getLocationMaximumGaussInterpolated(const DataSet& dataSet, std::size_t indexMin = 0UZ, std::size_t indexMax = 0UZ, std::size_t signalIndex = 0UZ) { + // getFrequencyEstimate => naive/simple approach counting edges + if (indexMax == 0UZ) { + indexMax = dataSet.axisValues(dim::X).size(); + } + detail::checkRangeIndex(indexMin, indexMax, dataSet.axisValues(dim::X).size()); + + const int locMax_ = getLocationMaximum(dataSet, indexMin, indexMax, signalIndex); + if (locMax_ <= static_cast(indexMin + 1) || locMax_ >= static_cast(indexMax - 1) || locMax_ < 0 || indexMax == indexMin) { + return std::numeric_limits::quiet_NaN(); + } + const std::size_t locMax = static_cast(locMax_); + + T refinedIndex = interpolateGaussian(dataSet.signalValues(signalIndex), locMax); + T x0 = getIndexValue(dataSet, dim::X, locMax, signalIndex); + // approximate step + // we’ll do a naive approach: X_{locMax+1} - X_{locMax} + if (locMax + 1UZ >= static_cast(dataSet.axisValues(dim::X).size())) { + return std::numeric_limits::quiet_NaN(); + } + T x1 = getIndexValue(dataSet, dim::X, locMax + 1UZ, signalIndex); + T diff = x1 - x0; + T deltaBin = refinedIndex - static_cast(locMax); + return x0 + deltaBin * diff; +} + +template> +[[nodiscard]] constexpr T getZeroCrossing(const DataSet& dataSet, TValue threshold, std::size_t signalIndex = 0) { + if (dataSet.axisCount() == 0 || dataSet.axis_values.empty() || dataSet.axis_values[0].empty()) { + return std::numeric_limits::quiet_NaN(); + } + std::size_t nLength = dataSet.axis_values[0].size(); + T initial = getIndexValue(dataSet, dim::Y, 0, signalIndex); + + if (initial < threshold) { // rising edge + for (std::size_t i = 1; i < nLength; ++i) { + T yPrev = getIndexValue(dataSet, dim::Y, i - 1, signalIndex); + T yCurr = getIndexValue(dataSet, dim::Y, i, signalIndex); + if (gr::math::isfinite(yPrev) && gr::math::isfinite(yCurr) && yCurr >= threshold) { + T xPrev = getIndexValue(dataSet, dim::X, i - 1, signalIndex); + T xCurr = getIndexValue(dataSet, dim::X, i, signalIndex); + return gr::dataset::detail::linearInterpolate(xPrev, xCurr, yPrev, yCurr, T(threshold)); + } + } + } else if (initial > threshold) { // falling edge + for (std::size_t i = 1; i < nLength; ++i) { + T yPrev = getIndexValue(dataSet, dim::Y, i - 1, signalIndex); + T yCurr = getIndexValue(dataSet, dim::Y, i, signalIndex); + if (gr::math::isfinite(yPrev) && gr::math::isfinite(yCurr) && yCurr <= threshold) { + T xPrev = getIndexValue(dataSet, dim::X, i - 1, signalIndex); + T xCurr = getIndexValue(dataSet, dim::X, i, signalIndex); + return gr::dataset::detail::linearInterpolate(xPrev, xCurr, yPrev, yCurr, T(threshold)); + } + } + } else { // exactly at the threshold + return getIndexValue(dataSet, dim::X, 0, signalIndex); + } + return std::numeric_limits::quiet_NaN(); // No crossing found +} + +} // namespace estimators + +} // namespace gr::dataset + +#endif // DATASETESTIMATORS_HPP diff --git a/algorithm/include/gnuradio-4.0/algorithm/dataset/DataSetHelper.hpp b/algorithm/include/gnuradio-4.0/algorithm/dataset/DataSetHelper.hpp new file mode 100644 index 00000000..c14191ba --- /dev/null +++ b/algorithm/include/gnuradio-4.0/algorithm/dataset/DataSetHelper.hpp @@ -0,0 +1,224 @@ +#ifndef DATASETHELPER_HPP +#define DATASETHELPER_HPP + +namespace gr::dataset { +namespace dim { +// N.B. explicitly and individually defined indices, rather than enum class +// since dimensionality is a priori open +inline static constexpr std::size_t X = 0UZ; /// X-axis index +inline static constexpr std::size_t Y = 1UZ; /// Y-axis index +inline static constexpr std::size_t Z = 2UZ; /// Z-axis index +} // namespace dim + +namespace detail { +inline void checkRangeIndex(std::size_t minIndex, std::size_t maxIndex, std::size_t dataCount, std::source_location location = std::source_location::current()) { + if (minIndex > maxIndex || maxIndex > dataCount) { + throw gr::exception(fmt::format("{} invalid range: [{},{}), dataCount={}", location, minIndex, maxIndex, dataCount)); + } +} + +template +[[nodiscard]] constexpr U linearInterpolate(T x0, T x1, U y0, U y1, U y) noexcept { + if (y1 == y0) { + return U(x0); + } + return x0 + (y - y0) * (x1 - x0) / (y1 - y0); +} + +template +[[nodiscard]] constexpr T sign(T positiveFactor, int val) noexcept { + return (val >= 0) ? positiveFactor : -positiveFactor; +} + +} // namespace detail + +template> +[[nodiscard]] constexpr T getIndexValue(const DataSet& dataSet, std::size_t dimIndex, std::size_t sampleIndex, std::size_t signalIndex = 0UZ) { + if (dimIndex == dim::X) { + if (dataSet.axisCount() == 0UZ || dataSet.axisValues(0UZ).size() <= sampleIndex) { + return std::numeric_limits::quiet_NaN(); + } + return static_cast(dataSet.axisValues(dim::X)[sampleIndex]); + } + + if (dimIndex == dim::Y) { + if (signalIndex >= dataSet.size()) { + return std::numeric_limits::quiet_NaN(); + } + std::span vals = dataSet.signalValues(signalIndex); + if (vals.size() <= sampleIndex) { + return std::numeric_limits::quiet_NaN(); + } + return static_cast(vals[sampleIndex]); + } + + throw gr::exception(fmt::format("axis dimIndex {} is out of range [0, {}] and", dimIndex, dataSet.axisCount())); +} + +template +T getDistance(const DataSet& dataSet, std::size_t dimIndex, std::size_t indexMin = 0UZ, std::size_t indexMax = 0UZ, std::size_t signalIndex = 0UZ) { + if (indexMax == 0UZ) { // renormalise default range + indexMax = dataSet.axisValues(dim::X).size() - 1UZ; + } + detail::checkRangeIndex(indexMin, indexMax, dataSet.axisValues(dim::X).size()); + + if (dimIndex == dim::X || dimIndex == dim::Y) { + return getIndexValue(dataSet, dimIndex, indexMax, signalIndex) - getIndexValue(dataSet, dimIndex, indexMin, signalIndex); + } + + throw gr::exception(fmt::format("axis dimIndex {} is out of range [0, {}] and", dimIndex, dataSet.axisCount())); +} + +/*! + * \brief Interpolate in Y dimension at a given X coordinate xValue. + * If out of range, we clamp or return NaN. + */ +template, std::convertible_to U> +[[nodiscard]] constexpr T getValue(const DataSet& dataSet, std::size_t dimIndex, U xValue, std::size_t signalIndex = 0) { + if (dimIndex == dim::X) { + return xValue; // if user requests X dimension at xValue, that’s basically identity + } + + if (dimIndex == dim::Y) { // linear interpolation in X-axis=0 + if (dataSet.axisCount() == 0) { + return std::numeric_limits::quiet_NaN(); + } + const auto& xAxis = dataSet.axisValues(dim::X); + const auto ySpan = dataSet.signalValues(signalIndex); + + if (xAxis.empty() || ySpan.empty()) { + return std::numeric_limits::quiet_NaN(); + } + + // clamp + if (xValue <= static_cast(xAxis.front())) { + return static_cast(ySpan.front()); + } else if (xValue >= static_cast(xAxis.back())) { + return static_cast(ySpan.back()); + } + + // binary-search for interval + auto it = std::lower_bound(xAxis.begin(), xAxis.end(), T(xValue)); + if (it == xAxis.end()) { + return static_cast(ySpan.back()); + } + auto idxHigh = std::distance(xAxis.begin(), it); + if (idxHigh == 0) { + return static_cast(ySpan.front()); + } + std::size_t iHigh = static_cast(idxHigh); + std::size_t iLow = iHigh - 1; + + T xLow = static_cast(xAxis[iLow]); + T xHigh = static_cast(xAxis[iHigh]); + T yLow = static_cast(ySpan[iLow]); + T yHigh = static_cast(ySpan[iHigh]); + + // linear interpolation + T dx = xHigh - xLow; + if (std::abs(gr::value(dx)) == T(0)) { + return yLow; + } + const T t = (xValue - xLow) / dx; + return yLow + t * (yHigh - yLow); + } + + throw gr::exception(fmt::format("axis dimIndex {} is out of range [0, {}] and", dimIndex, dataSet.axisCount())); +} + +template +std::vector getSubArrayCopy(const DataSet& ds, std::size_t indexMin, std::size_t indexMax, std::size_t signalIndex = 0) { + if (indexMax <= indexMin) { + return {}; + } + detail::checkRangeIndex(indexMin, indexMax, ds.axisValues(0).size()); + std::span values = ds.signalValues(signalIndex); + return std::vector(values.begin() + static_cast(indexMin), values.begin() + static_cast(indexMax)); +} + +template> +[[nodiscard]] constexpr T tenLog10(T x) noexcept { + if (std::abs(gr::value(x)) <= T(0)) { + return -std::numeric_limits::infinity(); + } + return T(10) * gr::math::log10(x); // 10 * log10(x) +} + +template> +[[nodiscard]] constexpr T decibel(T x) noexcept { + if (std::abs(gr::value(x)) <= T(0)) { + return -std::numeric_limits::infinity(); + } + return T(20) * gr::math::log10(x); // 20 * log10(x) +} + +template +[[maybe_unused]] constexpr bool verify(const DataSet& dataSet, std::source_location location = std::source_location::current()) { + auto handleFailure = [&](const std::string& message) -> bool { + if constexpr (ThrowOnFailure) { + throw gr::exception(message, location); + } + return false; + }; + + // axes checks + if (dataSet.axisCount() == 0UZ) { + return handleFailure("DataSet has zero axes."); + } + if (dataSet.axis_names.size() != dataSet.axisCount()) { + return handleFailure("axis_names size does not match axisCount()."); + } + if (dataSet.axis_units.size() != dataSet.axisCount()) { + return handleFailure("axis_units size does not match axisCount()."); + } + if (dataSet.axis_values.size() != dataSet.axisCount()) { + return handleFailure("axis_values size does not match axisCount()."); + } + for (std::size_t axis = 0; axis < dataSet.axisCount(); ++axis) { + if (dataSet.axis_values[axis].empty()) { + return handleFailure(fmt::format("axis_values[{}] is empty.", axis)); + } + } + + // verify extends + if (dataSet.extents.empty()) { + return handleFailure("DataSet has no extents defined."); + } + for (std::size_t dim = 0; dim < dataSet.extents.size(); ++dim) { + if (dataSet.extents[dim] <= 0) { + return handleFailure(fmt::format("Extent at dimension {} is non-positive.", dim)); + } + } + + // verify signal_names, signal_quantities, and signal_units sizes match extents[0] + std::size_t num_signals = static_cast(dataSet.extents[0]); + if (dataSet.signal_names.size() != num_signals) { + return handleFailure("signal_names size does not match extents[0]."); + } + if (dataSet.signal_quantities.size() != num_signals) { + return handleFailure("signal_quantities size does not match extents[0]."); + } + if (dataSet.signal_units.size() != num_signals) { + return handleFailure("signal_units size does not match extents[0]."); + } + + std::size_t expected_signal_size = std::accumulate(dataSet.extents.begin(), dataSet.extents.end(), static_cast(1), std::multiplies<>()); + if (dataSet.signal_values.size() != expected_signal_size) { + return handleFailure("signal_values size does not match the product of extents."); + } + if (dataSet.signal_ranges.size() != num_signals) { + return handleFailure("signal_ranges size does not match the number of signals."); + } + if (dataSet.meta_information.size() != num_signals) { // Assuming meta_information per signal value + return handleFailure("meta_information size does not match the number of signals."); + } + if (dataSet.timing_events.size() != num_signals) { // Assuming timing_events per signal value + return handleFailure("timing_events size does not match the number of signals."); + } + + return true; // all checks passed +} + +} // namespace gr::dataset + +#endif // DATASETHELPER_HPP diff --git a/algorithm/include/gnuradio-4.0/algorithm/dataset/DataSetMath.hpp b/algorithm/include/gnuradio-4.0/algorithm/dataset/DataSetMath.hpp new file mode 100644 index 00000000..99ad2224 --- /dev/null +++ b/algorithm/include/gnuradio-4.0/algorithm/dataset/DataSetMath.hpp @@ -0,0 +1,130 @@ +#ifndef DATASETMATH_HPP +#define DATASETMATH_HPP + +#include + +#include +#include +#include + +#include "DataSetHelper.hpp" + +namespace gr::dataset { + +enum class MathOp { ADD = 0, SUBTRACT, MULTIPLY, DIVIDE, SQR, SQRT, LOG10, DB, INV_DB, IDENTITY }; + +template +[[nodiscard]] constexpr bool sameHorizontalBase(const DataSet& ds1, const DataSet& ds2) { + if (ds1.axisCount() == 0 || ds2.axisCount() == 0) { + return false; + } + const auto& x1 = ds1.axisValues(0); + const auto& x2 = ds2.axisValues(0); + if (x1.size() != x2.size()) { + return false; + } + for (std::size_t i = 0; i < x1.size(); i++) { + if (x1[i] != x2[i]) { + return false; + } + } + return true; +} + +namespace detail { +template> +[[nodiscard]] constexpr T applyMathOperation(MathOp op, T y1, T y2) { + switch (op) { + case MathOp::ADD: return y1 + y2; + case MathOp::SUBTRACT: return y1 - y2; + case MathOp::MULTIPLY: return y1 * y2; + case MathOp::DIVIDE: return (y2 == TValue(0)) ? std::numeric_limits::quiet_NaN() : (y1 / y2); + case MathOp::SQR: return (y1 + y2) * (y1 + y2); + case MathOp::SQRT: return (y1 + y2) > TValue(0) ? gr::math::sqrt(y1 + y2) : std::numeric_limits::quiet_NaN(); + case MathOp::LOG10: return tenLog10((y1 + y2)); + case MathOp::DB: return decibel((y1 + y2)); + case MathOp::INV_DB: return inverseDecibel(y1); + case MathOp::IDENTITY: + default: return (y1 + y2); + } +} +} // namespace detail + +/*! + * \brief mathFunction(DataSet, DataSet, MathOp) => merges or interpolates ds2 onto ds1’s base + */ +template> +[[nodiscard]] constexpr DataSet mathFunction(const DataSet& ds1, const DataSet& ds2, MathOp op) { + + // Create new DataSet + DataSet ret; + ret.axis_names = ds1.axis_names; + ret.axis_units = ds1.axis_units; + ret.axis_values = ds1.axis_values; // or we do an interpolation base + ret.layout = ds1.layout; + ret.extents = ds1.extents; + ret.signal_names = {"mathOp"}; + ret.signal_quantities = {ds1.signal_quantities.empty() ? "quantity" : ds1.signal_quantities[0]}; + ret.signal_units = {ds1.signal_units.empty() ? "" : ds1.signal_units[0]}; + ret.signal_ranges = {{T(0.0), T(1.0)}}; // placeholder + ret.meta_information = ds1.meta_information; // or merge + ret.timing_events = ds1.timing_events; + + std::size_t dataCount = 0U; + if (!ret.axis_values.empty()) { + dataCount = ret.axis_values[0].size(); + } + ret.signal_values.resize(dataCount); + + bool needsInterpolation = !sameHorizontalBase(ds1, ds2); + + for (std::size_t i = 0; i < dataCount; i++) { + TValue x = gr::value(getIndexValue(ds1, dim::X, i)); + T Y1 = getIndexValue(ds1, dim::Y, i); + T Y2 = needsInterpolation ? getValue(ds2, dim::Y, x) : getIndexValue(ds2, dim::Y, i); + T result = detail::applyMathOperation(op, Y1, Y2); + ret.signal_values[i] = static_cast(result); + } + return ret; +} + +template U, typename TValue = gr::meta::fundamental_base_value_type_t> +[[nodiscard]] constexpr DataSet mathFunction(const DataSet& ds, U value, MathOp op) { + + DataSet ret = ds; // copy + auto dataCount = ds.axis_values.empty() ? 0UZ : ds.axis_values[0].size(); + for (std::size_t i = 0; i < dataCount; i++) { + T Y1 = getIndexValue(ds, dim::Y, i); + T result = 0; + switch (op) { + case MathOp::ADD: result = Y1 + value; break; + case MathOp::SUBTRACT: result = Y1 - value; break; + case MathOp::MULTIPLY: result = Y1 * value; break; + case MathOp::DIVIDE: result = (value == TValue(0)) ? std::numeric_limits::quiet_NaN() : (Y1 / value); break; + case MathOp::SQR: result = (Y1 + value) * (Y1 + value); break; + case MathOp::SQRT: result = (Y1 + value) > TValue(0) ? gr::math::sqrt(Y1 + value) : std::numeric_limits::quiet_NaN(); break; + case MathOp::LOG10: result = tenLog10(Y1 + value); break; + case MathOp::DB: result = decibel(Y1 + value); break; + case MathOp::INV_DB: result = inverseDecibel(Y1); break; + case MathOp::IDENTITY: + default: result = Y1; break; + } + ret.signal_values[i] = static_cast(result); + } + return ret; +} + +template U> +[[nodiscard]] constexpr DataSet addFunction(const DataSet& ds, U value) { + return mathFunction(ds, value, MathOp::ADD); +} +template +[[nodiscard]] constexpr DataSet addFunction(const DataSet& ds1, const DataSet& ds2) { + return mathFunction(ds1, ds2, MathOp::ADD); +} + +// WIP complete other convenience and missing math functions + +} // namespace gr::dataset + +#endif // DATASETMATH_HPP diff --git a/algorithm/include/gnuradio-4.0/algorithm/dataset/DataSetTestFunctions.hpp b/algorithm/include/gnuradio-4.0/algorithm/dataset/DataSetTestFunctions.hpp new file mode 100644 index 00000000..7381318c --- /dev/null +++ b/algorithm/include/gnuradio-4.0/algorithm/dataset/DataSetTestFunctions.hpp @@ -0,0 +1,162 @@ +#ifndef DATASETTESTFUNCTIONS_HPP +#define DATASETTESTFUNCTIONS_HPP + +#include +#include +#include + +#include "DataSetHelper.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace gr::dataset::generate { + +namespace detail { +template +constexpr auto initialize = [](typename gr::meta::fundamental_base_value_type_t value, typename gr::meta::fundamental_base_value_type_t uncertainty = {}) -> T { + if constexpr (gr::UncertainValueLike) { + return T{value, uncertainty}; + } else { + return value; + } +}; +} // namespace detail + +template> +[[nodiscard]] constexpr gr::DataSet triangular(std::string name, std::size_t count, TValue offset = 0, TValue amplitude = 1) { + assert(count > 2UZ); + gr::DataSet ds; + ds.signal_names = {name}; + ds.signal_quantities = {"Amplitude"}; + ds.signal_units = {""}; + ds.axis_names = {"Time"}; + ds.axis_units = {"s"}; + ds.axis_values.resize(1UZ); // one X-axis + ds.axis_values[gr::dataset::dim::X].resize(count); + ds.signal_values.resize(count); + ds.meta_information.resize(1UZ); + ds.timing_events.resize(1UZ); + ds.extents = {1, static_cast(count)}; + + for (std::size_t i = 0UZ; i < count; ++i) { + ds.axis_values[0][i] = static_cast(i); + } + + std::size_t midpointLeft = count / 2; + for (std::size_t i = 0UZ; i < midpointLeft; ++i) { + TValue factor = static_cast(i) / static_cast(midpointLeft - (count & 1 ? 0UZ : 1UZ)); + ds.signal_values[i] = detail::initialize(offset + amplitude * factor, amplitude / TValue(10)); + ds.signal_values[count - i - 1] = ds.signal_values[i]; // ensures symmetry + } + if (count & 1) { // centre point + ds.signal_values[midpointLeft] = detail::initialize(offset + amplitude, amplitude / TValue(10)); + } + + ds.signal_ranges.push_back({T(offset), T(amplitude)}); + return ds; +} + +template> +[[nodiscard]] constexpr gr::DataSet ramp(std::string name, std::size_t count, TValue offset = TValue(0), TValue amplitude = TValue(1)) { + gr::DataSet ds; + ds.signal_names = {name}; + ds.signal_quantities = {"Amplitude"}; + ds.signal_units = {""}; + ds.axis_names = {"Time"}; + ds.axis_units = {"s"}; + ds.axis_values.resize(1); + ds.axis_values[0].resize(count); + ds.signal_values.resize(count); + ds.meta_information.resize(1UZ); + ds.timing_events.resize(1UZ); + ds.extents = {1, static_cast(count)}; + + for (std::size_t i = 0; i < count; i++) { + ds.axis_values[0][i] = static_cast(i); + ds.signal_values[i] = detail::initialize(TValue(offset) + TValue(amplitude) * TValue(i) / TValue(count), TValue(amplitude) / TValue(10)); + } + ds.signal_ranges.push_back({T(offset), T(amplitude)}); + return ds; +} + +template, std::convertible_to U> +[[nodiscard]] constexpr gr::DataSet gaussFunction(std::string name, std::size_t count, U mean = TValue(0), U sigma = U(3), U offset = U(0), U amplitude = U(1)) { + if (count <= 2) { + throw std::invalid_argument("Count must be greater than 2"); + } + if (sigma <= 0) { + throw std::invalid_argument("Sigma must be positive"); + } + + gr::DataSet ds; + ds.signal_names = {name}; + ds.signal_quantities = {"Amplitude"}; + ds.signal_units = {""}; + ds.axis_names = {"Time"}; + ds.axis_units = {"s"}; + ds.axis_values.resize(1); + ds.axis_values[0].resize(count); + ds.signal_values.resize(count); + ds.meta_information.resize(1UZ); + ds.timing_events.resize(1UZ); + ds.extents = {1, static_cast(count)}; + + auto gauss = [](U x, U mu, U sigma) -> TValue { return std::exp(-std::pow((TValue(x) - TValue(mu)) / TValue(sigma), 2) / U(2)) / (TValue(sigma) * std::sqrt(TValue(2) * std::numbers::pi_v)); }; + + for (std::size_t i = 0; i < count; ++i) { + ds.axis_values[0][i] = static_cast(i); + const TValue val = gauss(static_cast(i), U(mean), U(sigma)) * TValue(amplitude) + TValue(offset); + ds.signal_values[i] = detail::initialize(val, TValue(amplitude) / TValue(10)); + } + + auto minmax = std::minmax_element(ds.signal_values.begin(), ds.signal_values.end()); + ds.signal_ranges.push_back(Range{*minmax.first, *minmax.second}); + + return ds; +} + +template> +[[nodiscard]] gr::DataSet randomStepFunction(std::string name, std::size_t count, std::uint64_t seed = 0U) { + if (count == 0) { + throw std::invalid_argument("Count must be greater than 0"); + } + + gr::DataSet ds; + ds.signal_names = {name}; + ds.signal_quantities = {"Amplitude"}; + ds.signal_units = {""}; + ds.axis_names = {"Time"}; + ds.axis_units = {"s"}; + ds.axis_values.resize(1); + ds.axis_values[0].resize(count); + ds.signal_values.resize(count); + ds.meta_information.resize(1); + ds.timing_events.resize(1); + ds.extents = {1, static_cast(count)}; + + std::random_device rd; + std::mt19937_64 rng(seed == 0 ? rd() : seed); + std::uniform_int_distribution dist(0, count - 1); + std::size_t step = dist(rng); + + for (std::size_t i = 0; i < count; ++i) { + ds.axis_values[0][i] = static_cast(i); + TValue val = static_cast(i < step ? 0.0 : 1.0); + ds.signal_values[i] = detail::initialize(static_cast(val), static_cast(1) / static_cast(10)); + } + + ds.signal_ranges.push_back(Range{static_cast(0), static_cast(1)}); + + return ds; +} + +} // namespace gr::dataset::generate + +#endif // DATASETTESTFUNCTIONS_HPP diff --git a/algorithm/include/gnuradio-4.0/algorithm/dataset/DataSetUtils.hpp b/algorithm/include/gnuradio-4.0/algorithm/dataset/DataSetUtils.hpp index 574148a6..50da41a6 100644 --- a/algorithm/include/gnuradio-4.0/algorithm/dataset/DataSetUtils.hpp +++ b/algorithm/include/gnuradio-4.0/algorithm/dataset/DataSetUtils.hpp @@ -17,9 +17,9 @@ template [[maybe_unused]] bool draw(const TDataSet& dataSet, const DefaultChartConfig config = {}) { using TValueType = typename TDataSet::value_type; - if (dataSet.signal_values.empty() // check for empty data - || dataSet.axis_values.empty() || dataSet.axis_values[0].empty() // empty axis definition - || dataSet.signal_ranges.empty() || dataSet.signal_ranges[0].empty() // empty min/max definition + if (dataSet.signal_values.empty() // check for empty data + || dataSet.axis_values.empty() || dataSet.axis_values[0].empty() // empty axis definition + || dataSet.signal_ranges.empty() // empty min/max definition ) { return false; } @@ -38,8 +38,8 @@ template assert(!dataSet.signal_ranges.empty()); const TValueType xMin = dataSet.axis_values[0].front(); const TValueType xMax = dataSet.axis_values[0].back(); - TValueType yMin = dataSet.signal_ranges[0].front(); - TValueType yMax = dataSet.signal_ranges[0].back(); + TValueType yMin = dataSet.signal_ranges[0].min; + TValueType yMax = dataSet.signal_ranges[0].max; if constexpr (std::is_arithmetic_v) { const auto [min, max] = std::ranges::minmax_element(dataSet.signal_values); @@ -92,15 +92,15 @@ template template void updateMinMax(DataSet& dataSet) { if constexpr (std::is_arithmetic_v) { - const auto [min, max] = std::ranges::minmax_element(dataSet.signal_values); - dataSet.signal_ranges[0][0] = *min; - dataSet.signal_ranges[0][1] = *max; + const auto [min, max] = std::ranges::minmax_element(dataSet.signal_values); + dataSet.signal_ranges[0].min = *min; + dataSet.signal_ranges[0].max = *max; } else if constexpr (gr::meta::complex_like) { const auto [min, max] = std::ranges::minmax_element(dataSet.signal_values, // [](const T& a, const T& b) { return std::abs(a) < std::abs(b); }); - dataSet.signal_ranges[0][0] = *min; - dataSet.signal_ranges[0][1] = *max; + dataSet.signal_ranges[0].min = *min; + dataSet.signal_ranges[0].max = *max; } else { static_assert(std::is_arithmetic_v || std::is_same_v>, "Unsupported type for DataSet"); } diff --git a/algorithm/test/CMakeLists.txt b/algorithm/test/CMakeLists.txt index 685b8b3a..60eea202 100644 --- a/algorithm/test/CMakeLists.txt +++ b/algorithm/test/CMakeLists.txt @@ -1,8 +1,10 @@ add_ut_test(qa_algorithm_fourier) +add_ut_test(qa_DataSetEstimators) add_ut_test(qa_FilterTool) add_ut_test(qa_ImChart) add_ut_test(qa_SchmittTrigger) target_link_libraries(qa_algorithm_fourier PRIVATE gnuradio-algorithm) +target_link_libraries(qa_DataSetEstimators PRIVATE gnuradio-algorithm) target_link_libraries(qa_FilterTool PRIVATE gnuradio-algorithm) target_link_libraries(qa_ImChart PRIVATE gnuradio-algorithm) target_link_libraries(qa_SchmittTrigger PRIVATE gnuradio-algorithm) diff --git a/algorithm/test/qa_DataSetEstimators.cpp b/algorithm/test/qa_DataSetEstimators.cpp new file mode 100644 index 00000000..8beccd18 --- /dev/null +++ b/algorithm/test/qa_DataSetEstimators.cpp @@ -0,0 +1,380 @@ +#include +#include +#include +#include // for draw(...) +#include + +#include + +#include +#include +#include + +namespace test::detail { + +template +[[nodiscard]] constexpr auto approx(const TLhs& lhs, const TRhs& rhs, const TEpsilon& epsilon) { + if constexpr (gr::meta::complex_like) { + return boost::ut::detail::and_{boost::ut::detail::approx_{lhs.real(), rhs.real(), epsilon}, boost::ut::detail::approx_{lhs.imag(), rhs.imag(), epsilon}}; + } else { + return boost::ut::detail::approx_{lhs, rhs, epsilon}; + } +} +} // namespace test::detail + +const boost::ut::suite<"DataSet visual test functions"> _DataSetTestFcuntions = [] { + using namespace boost::ut; + using namespace gr::dataset; + constexpr static std::size_t nSamples = 201UZ; + + "triangular DataSet"_test = [] { + gr::DataSet ds = generate::triangular("triagonal", nSamples); + gr::dataset::draw(ds); + + expect(gr::dataset::verify(ds)); + + gr::DataSet ds1 = generate::triangular("triagonal - odd", 11); + fmt::println("\"{:20}\": {}", ds1.signalName(0UZ), ds1.signal_values); + expect(eq(ds1.signalValues().front(), ds1.signalValues().back())); + expect(eq(ds1.signalValues()[5UZ], 1.0)); + + gr::DataSet ds2 = generate::triangular("triagonal - even", 10); + fmt::println("\"{:20}\": {}", ds2.signalName(0UZ), ds2.signal_values); + expect(eq(ds2.signalValues().front(), ds2.signalValues().back())); + expect(eq(ds2.signalValues()[4UZ], ds2.signalValues()[5UZ])); + }; + + "ramp DataSet"_test = [] { + gr::DataSet ds = generate::ramp("ramp", nSamples); + expect(gr::dataset::verify(ds)); + gr::dataset::draw(ds); + }; + + "gaussFunction DataSet"_test = [] { + constexpr T mean = T(nSamples) / T(2); + constexpr T sigma = T(nSamples) / T(10); + const T normalisation = sigma * gr::math::sqrt(T(2) * std::numbers::pi_v); + + gr::DataSet ds = generate::gaussFunction("gaussFunction", nSamples, mean, sigma, T(0), normalisation); + expect(gr::dataset::verify(ds)); + gr::dataset::draw(ds); + }; + + "ramp + gauss DataSet"_test = [] { + constexpr T mean = T(nSamples) / T(2); + constexpr T sigma = T(nSamples) / T(10); + const T normalisation = sigma * gr::math::sqrt(T(2) * std::numbers::pi_v); + + gr::DataSet ds1 = generate::ramp("ramp", nSamples, 0.0, 0.2); + gr::DataSet ds2 = generate::gaussFunction("gaussFunction", nSamples, mean, sigma, T(0), normalisation); + gr::DataSet ds = addFunction(ds1, ds2); + + expect(gr::dataset::verify(ds)); + gr::dataset::draw(ds); + }; + + "randomStepFunction DataSet"_test = [] { + gr::DataSet ds = generate::randomStepFunction("randomStepFunction", nSamples); + expect(gr::dataset::verify(ds)); + gr::dataset::draw(ds); + }; +}; + +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wdouble-promotion" +#pragma GCC diagnostic ignored "-Wimplicit-float-conversion" // GCC & Clang: Loss or promotion of floating-point precision, disabled only for unit-tests + +const boost::ut::suite<"DataSet element-wise accessor"> _dataSetAccessors = [] { + using namespace boost::ut; + using namespace gr::dataset; + using test::detail::approx; + + "basic element-wise access API "_test = [] { + using value_t = gr::meta::fundamental_base_value_type_t; + auto ds = generate::triangular("triag", 3UZ); // data ~ [0, 1, 0] + + expect(approx(getIndexValue(ds, dim::Y, 0), value_t(0), value_t(1e-3f))); + expect(approx(getIndexValue(ds, dim::Y, 1), value_t(1), value_t(1e-3f))); + expect(approx(getIndexValue(ds, dim::Y, 0), value_t(0), value_t(1e-3f))); + + expect(!gr::math::isfinite(getIndexValue(ds, dim::X, 3))) << fmt::format("element is not NaN: {}", getIndexValue(ds, dim::X, 3)); + expect(!gr::math::isfinite(getIndexValue(ds, dim::Y, 3))) << fmt::format("element is not NaN: {}", getIndexValue(ds, dim::Y, 3)); + + expect(throws([&] { std::ignore = getIndexValue(ds, dim::Z, 0); })); + + expect(approx(getDistance(ds, dim::X, 0UZ, 2UZ), value_t(2), value_t(1e-3f))); + expect(approx(getDistance(ds, dim::X), value_t(2), value_t(1e-3f))); + + expect(approx(getDistance(ds, dim::Y, 0UZ, 1UZ), value_t(1), value_t(1e-3f))); // Y-distance of first to middle + expect(approx(getDistance(ds, dim::Y), T(0), value_t(1e-3f))); // Y-distance of first to last + + expect(approx(getValue(ds, dim::X, value_t(0.123f)), static_cast(0.123f), value_t(1e-3f))); // identity + expect(approx(getValue(ds, dim::Y, value_t(0.5)), static_cast(0.5f), value_t(1e-3f))); + + std::vector copyX = getSubArrayCopy(ds, dim::X, 0UZ, 2UZ); + for (std::size_t i = 0; i < copyX.size(); i++) { + expect(eq(copyX[i], getIndexValue(ds, dim::X, i))) << fmt::format("X-index {} mismatch", i); + } + + std::vector copyY = getSubArrayCopy(ds, dim::X, 0UZ, 2UZ); + for (std::size_t i = 0; i < copyY.size(); i++) { + expect(eq(copyY[i], getIndexValue(ds, dim::Y, i))) << fmt::format("X-index {} mismatch", i); + } + } | std::tuple, gr::UncertainValue>{}; +}; + +const boost::ut::suite<"DSP helper"> _dspHelper = [] { + using namespace boost::ut; + using namespace gr::dataset; + using test::detail::approx; + + "common dsp helper functions"_test = [] { + using value_t = T; + const static std::string typeName = gr::meta::type_name(); + + "tenLog10"_test = [] { + expect(approx(tenLog10(value_t(10)), value_t(10), value_t(1e-3f))) << "10 * log10(10) should be 10"; + expect(approx(tenLog10(value_t(1)), value_t(0), value_t(1e-3f))) << "10 * log10(1) should be 0"; + expect(approx(tenLog10(value_t(0.1)), value_t(-10), value_t(1e-3f))) << "10 * log10(0.1) should be -10"; + + // edge cases + expect(!gr::math::isfinite(tenLog10(value_t(0)))) << fmt::format("tenLog10<{}>(0) = {} should be -inf", typeName, tenLog10(value_t(0))); + expect(!gr::math::isfinite(tenLog10(value_t(-1)))) << fmt::format("tenLog10<{}>(-1) = {} should be -inf", typeName, tenLog10(value_t(-1))); + }; + + "decibel"_test = [] { + expect(approx(decibel(value_t(10)), value_t(20), value_t(1e-3f))) << "20 * log10(10) should be 20"; + expect(approx(decibel(value_t(1)), value_t(0), value_t(1e-3f))) << "20 * log10(1) should be 0"; + expect(approx(decibel(value_t(0.1)), value_t(-20), value_t(1e-3f))) << "20 * log10(0.1) should be -20"; + + // edge cases + expect(!gr::math::isfinite(decibel(value_t(0)))) << "decibel(0) should be -inf"; + expect(!gr::math::isfinite(decibel(value_t(-1)))) << "20 * log10(-1) should return -inf"; + }; + + "inverseDecibel"_test = [] { + expect(approx(inverseDecibel(value_t(20)), value_t(10), value_t(1e-3f))) << "10^(20 / 20) should be 10"; + expect(approx(inverseDecibel(value_t(0)), value_t(1), value_t(1e-3f))) << "10^(0 / 20) should be 1"; + expect(approx(inverseDecibel(value_t(-20)), value_t(0.1), value_t(1e-3f))) << "10^(-20 / 20) should be 0.1"; + + // edge cases + expect(approx(inverseDecibel(value_t(100)), gr::math::pow(value_t(10), value_t(5)), value_t(1e-3f))) << "10^(100 / 20) should be 10^5"; + expect(approx(inverseDecibel(value_t(-100)), gr::math::pow(value_t(10), value_t(-5)), value_t(1e-3f))) << "10^(-100 / 20) should be 10^-5"; + }; + } | std::tuple, gr::UncertainValue>{}; +}; + +const boost::ut::suite<"DataSet estimator"> _qaDataSetEstimators = [] { + using namespace boost::ut; + using namespace gr::dataset; + using test::detail::approx; + + constexpr static size_t nSamples = 11; + + "basic estimators"_test = [] { + using value_t = gr::meta::fundamental_base_value_type_t; + // 10-sample triangle => ascending up to i=4..5 => data ~ [0, 0.2, 0.4, 0.6, 0.8, 1, 0.8, 0.6, 0.4, 0.2, 0] + auto ds = generate::triangular("triag", nSamples); + + expect(approx(estimators::computeCentreOfMass(ds), T(5), T(1e-3))); + expect(gr::math::isfinite(estimators::computeCentreOfMass(ds, 0UZ, nSamples))) << "should be finite for partial range"; + + std::vector data{1, 2, 3, 2, 1}; + expect(approx(estimators::computeFWHM(data, 2), T(4), T(1e-5))); + expect(approx(estimators::computeInterpolatedFWHM(data, 2), T(3), T(1e-5))); + + expect(eq(estimators::getLocationMaximum(ds, 0UZ, nSamples), int(5))); + expect(eq(estimators::getLocationMinimum(ds, 0UZ, nSamples), int(0))); + expect(eq(estimators::getLocationMaximum(ds), int(5))); + expect(eq(estimators::getLocationMinimum(ds), int(0))); + + expect(eq(gr::value(estimators::getMaximum(ds, 0UZ, nSamples)), value_t(1))); + expect(eq(gr::value(estimators::getMaximum(ds)), value_t(1))); + expect(eq(gr::value(estimators::getMinimum(ds, 0UZ, nSamples)), value_t(0))); + expect(eq(gr::value(estimators::getMinimum(ds)), value_t(0))); + + expect(approx(estimators::getMean(ds, 0UZ, nSamples), T(0.454545), T(1e-3f))); + expect(approx(estimators::getMean(ds), T(0.454545), T(1e-3f))); + + expect(approx(estimators::getMedian(ds, 0UZ, nSamples), T(0.4), T(1e-3f))); + expect(approx(estimators::getMedian(ds), T(0.4), T(1e-3f))); + + expect(eq(gr::value(estimators::getRange(ds, 0UZ, nSamples)), value_t(1))); + expect(eq(gr::value(estimators::getRange(ds)), value_t(1))); + + expect(approx(estimators::getRms(ds, 0UZ, nSamples), T(0.320124), T(1e-3))); + expect(approx(estimators::getRms(ds), T(0.320124), T(1e-3))); + if constexpr (gr::UncertainValueLike) { + expect(neq(gr::uncertainty(estimators::getRms(ds, 0UZ, nSamples)), value_t(0))); + expect(neq(gr::uncertainty(estimators::getRms(ds)), value_t(0))); + } + + expect(approx(estimators::getIntegral(ds, 0UZ, nSamples), T(5.0), T(1e-3))); + expect(approx(estimators::getIntegral(ds), T(5.0), T(1e-3))); + if constexpr (gr::UncertainValueLike) { + expect(neq(gr::uncertainty(estimators::getIntegral(ds, 0UZ, nSamples)), value_t(0))); + expect(neq(gr::uncertainty(estimators::getIntegral(ds)), value_t(0))); + } + + expect(approx(estimators::getEdgeDetect(ds, 0UZ, nSamples), T(3), T(0.5))) << "50% is ~ 0.5 => crossing near i=3"; + expect(approx(estimators::getEdgeDetect(ds), T(3), T(0.5))) << "50% is ~ 0.5 => crossing near i=3"; + } | std::tuple, gr::UncertainValue>{}; + ; + + "getDutyCycle"_test = [] { + using value_t = gr::meta::fundamental_base_value_type_t; + std::vector localY = {0, 0, 0, 1, 1, 1}; // simple data set with 0,0,0,1,1,1 => 50% high + gr::DataSet ds; + ds.axis_names = {"Time"}; + ds.axis_units = {"s"}; + ds.axis_values.resize(1); + ds.axis_values[0].resize(localY.size()); + ds.signal_names = {"simple step"}; + ds.signal_values.resize(localY.size()); + ds.extents = {1, static_cast(localY.size())}; + for (std::size_t i = 0; i < localY.size(); i++) { + ds.axis_values[0][i] = T(i); + ds.signal_values[i] = localY[i]; + } + ds.signal_ranges.push_back({T(0), T(1)}); + + expect(eq(estimators::getDutyCycle(ds, 0UZ, localY.size()), value_t(0.5))) << "simple 3-high/3-low => duty=0.5"; + expect(eq(estimators::getDutyCycle(ds), value_t(0.5))) << "simple 3-high/3-low => duty=0.5"; + } | std::tuple, gr::UncertainValue>{}; + + "getFrequencyEstimate"_test = [] { + using value_t = gr::meta::fundamental_base_value_type_t; + std::vector localY = {0, 1, 0, 1, 0, 1}; + gr::DataSet ds; + ds.axis_names = {"Time"}; + ds.axis_units = {"s"}; + ds.axis_values.resize(1); + ds.axis_values[0].resize(localY.size()); + ds.signal_names = {"oscillator"}; + ds.signal_values.resize(localY.size()); + ds.extents = {1, int(localY.size())}; + for (std::size_t i = 0; i < localY.size(); i++) { + ds.axis_values[0][i] = T(i); + ds.signal_values[i] = localY[i]; + } + ds.signal_ranges.push_back({T(0), T(1)}); + + expect(eq(estimators::getFrequencyEstimate(ds, 0UZ, localY.size()), value_t(0.5))) << "frequency ~ 0.5"; + expect(eq(estimators::getFrequencyEstimate(ds), value_t(0.5))) << "frequency ~ 0.5"; + } | std::tuple, gr::UncertainValue>{}; + + "getLocationMaximumGaussInterpolated"_test = [] { + using value_t = gr::meta::fundamental_base_value_type_t; + auto ds1 = generate::triangular("triag", 7UZ); // symmetric peak around index 3 + + expect(approx(gr::value(estimators::getLocationMaximumGaussInterpolated(ds1, 0UZ, 7UZ)), value_t(3), T(1e-3f))) << "Gauss interpolation ~3"; + expect(approx(gr::value(estimators::getLocationMaximumGaussInterpolated(ds1)), value_t(3), T(1e-3f))) << "Gauss interpolation ~3"; + if constexpr (gr::UncertainValueLike) { + expect(neq(gr::uncertainty(estimators::getLocationMaximumGaussInterpolated(ds1, 0UZ, 7UZ)), value_t(0))); + expect(neq(gr::uncertainty(estimators::getLocationMaximumGaussInterpolated(ds1)), value_t(0))); + } + + auto ds2 = generate::triangular("triag", 6UZ); // symmetric peak around index 2.5 + expect(approx(gr::value(estimators::getLocationMaximumGaussInterpolated(ds2, 0UZ, 6UZ)), value_t(2.5), T(1e-3f))) << "Gauss interpolation ~2.5"; + expect(approx(gr::value(estimators::getLocationMaximumGaussInterpolated(ds2)), value_t(2.5), T(1e-3f))) << "Gauss interpolation ~2.5"; + if constexpr (gr::UncertainValueLike) { + expect(neq(gr::uncertainty(estimators::getLocationMaximumGaussInterpolated(ds2, 0UZ, 6UZ)), value_t(0))); + expect(neq(gr::uncertainty(estimators::getLocationMaximumGaussInterpolated(ds2)), value_t(0))); + } + } | std::tuple, gr::UncertainValue>{}; + + "zeroCrossing test"_test = [] { + using value_t = gr::meta::fundamental_base_value_type_t; + auto ds = generate::triangular("ZC", nSamples); + + expect(approx(gr::value(estimators::getZeroCrossing(ds, value_t(0.5))), value_t(2.5), value_t(1e-3))) << "zero crossing ~2.5"; + if constexpr (gr::UncertainValueLike) { + expect(neq(gr::uncertainty(estimators::getZeroCrossing(ds, value_t(0.5))), value_t(0))); + } + } | std::tuple, gr::UncertainValue>{}; +}; + +const boost::ut::suite<"DataSet math "> _dataSetMath = [] { + using namespace boost::ut; + using namespace boost::ut::literals; + using namespace gr::dataset; + + "basic math API "_test = [] { + using value_t = gr::meta::fundamental_base_value_type_t; + auto ds1 = generate::triangular("ds1", 5); + auto ds2 = addFunction(ds1, value_t(2)); // add scalar '2' + + for (std::size_t i = 0; i < 5; i++) { + T y1 = gr::dataset::getIndexValue(ds1, dim::Y, i); + T y2 = gr::dataset::getIndexValue(ds2, dim::Y, i); + expect(approx(y2, y1 + T(2), value_t(1e-3f))); + } + + // add two dataSets + auto ds3 = addFunction(ds1, ds1); + for (std::size_t i = 0; i < 5; i++) { + T y1 = gr::dataset::getIndexValue(ds1, 1, i); + T y3 = gr::dataset::getIndexValue(ds3, 1, i); + expect(approx(y3, value_t(2) * y1, value_t(1e-3f))); + } + } | std::tuple, gr::UncertainValue>{}; + + "mathFunction variants"_test = [] { + auto ds1 = generate::triangular("triag", 5); + + // 1) dataSet + dataSet + auto dsAdd = mathFunction(ds1, ds1, MathOp::ADD); + for (std::size_t i = 0UZ; i < 5UZ; i++) { + T origVal = ds1.signal_values[i]; + T newVal = dsAdd.signal_values[i]; + expect(eq(newVal, T(2) * origVal)) << fmt::format("Add op failed at i={}, got {}", i, newVal); + } + + // 2) dataSet + double + auto dsAdd2 = mathFunction(ds1, T(2), MathOp::ADD); + for (std::size_t i = 0UZ; i < 5UZ; i++) { + T origVal = ds1.signal_values[i]; + T newVal = dsAdd2.signal_values[i]; + expect(eq(newVal, origVal + T(2))); + } + + // 3) MULTIPLY + auto dsMul = mathFunction(ds1, ds1, MathOp::MULTIPLY); + for (std::size_t i = 0UZ; i < 5UZ; i++) { + T origVal = ds1.signal_values[i]; + T newVal = dsMul.signal_values[i]; + expect(approx(newVal, origVal * origVal, T(1e-5))) << "Multiply test"; + } + + // 4) SQR => (y1 + y2)^2 => ds1 + 2 => then squared + auto dsSqr = mathFunction(ds1, T(2), MathOp::SQR); + // check a single sample + expect(eq(dsSqr.signal_values[0], T((ds1.signal_values[0] + T(2)) * (ds1.signal_values[0] + T(2))))); + + // 5) SQRT => sqrt(y1 + y2) + auto dsSqrt = mathFunction(ds1, T(2.0), MathOp::SQRT); + // if ds1.signal_values[0]=0 => sqrt(2) => ~1.4142 + expect(approx(dsSqrt.signal_values[0], T(1.4142f), T(1e-2))); + + // 6) LOG10 => 10*log10(y1 + y2) + auto dsLog = mathFunction(ds1, T(2.0), MathOp::LOG10); + // we can do a quick check for i=0 => 10 * log10(2.0) => 3.0103 dB + expect(approx(dsLog.signal_values[0], T(3.01f), T(1e-2))); + + // 7) DB => 20*log10(y1 + y2) + auto dsDb = mathFunction(ds1, T(2.0), MathOp::DB); + // i=0 => 20*log10(2.0) => 6.0206 dB + expect(approx(dsDb.signal_values[0], T(6.02f), T(1e-2))); + + // 8) INV_DB => 10^(y1/20), ignoring y2 + // if ds1.signal_values[0]=0 => => 10^(0/20)=>1 + auto dsInv = mathFunction(ds1, 2.0, MathOp::INV_DB); + expect(eq(dsInv.signal_values[0], T(1))) << "inv_db test"; + }; + + // WIP -- add more DataSet math-related functionality here +}; + +#pragma GCC diagnostic pop + +int main() { /* not needed for UT */ } diff --git a/blocks/basic/include/gnuradio-4.0/basic/StreamToDataSet.hpp b/blocks/basic/include/gnuradio-4.0/basic/StreamToDataSet.hpp index 7552501f..0b93c4eb 100644 --- a/blocks/basic/include/gnuradio-4.0/basic/StreamToDataSet.hpp +++ b/blocks/basic/include/gnuradio-4.0/basic/StreamToDataSet.hpp @@ -342,9 +342,8 @@ If multiple 'start' or 'stop' Tags arrive in a single merged tag, only one DataS dataSet.signal_names.emplace_back(signal_name); dataSet.signal_quantities.emplace_back(signal_quantity); dataSet.signal_units.emplace_back(signal_unit); - dataSet.signal_ranges.resize(1UZ); // one data set - dataSet.signal_ranges[0].resize(2UZ); // [min, max]s - dataSet.meta_information.resize(1); // one data set + dataSet.signal_ranges.resize(1UZ); // one data set + dataSet.meta_information.resize(1); // one data set dataSet.meta_information[0]["ctx"] = filter; dataSet.meta_information[0]["n_pre"] = n_pre; dataSet.meta_information[0]["n_post"] = n_post; diff --git a/blocks/basic/test/qa_DataSink.cpp b/blocks/basic/test/qa_DataSink.cpp index 3c808150..61f8d1e8 100644 --- a/blocks/basic/test/qa_DataSink.cpp +++ b/blocks/basic/test/qa_DataSink.cpp @@ -496,7 +496,7 @@ const boost::ut::suite DataSinkTests = [] { expect(eq(dataset.timing_events.size(), 1u)); expect(eq(dataset.signal_names[0], "test signal"s)); expect(eq(dataset.signal_units[0], "none"s)); - expect(eq(dataset.signal_ranges[0], std::vector{-2, +2})); + expect(eq(dataset.signal_ranges[0], gr::Range{-2, +2})); expect(eq(dataset.timing_events[0].size(), 1u)); expect(eq(dataset.timing_events[0][0].first, 3)); receivedTags.insert(receivedTags.end(), dataset.timing_events[0].begin(), dataset.timing_events[0].end()); @@ -563,7 +563,7 @@ const boost::ut::suite DataSinkTests = [] { expect(eq(dataset.timing_events.size(), 1u)); expect(eq(dataset.signal_names[0], "test signal"s)); expect(eq(dataset.signal_units[0], "no unit"s)); - expect(eq(dataset.signal_ranges[0], std::vector{-2, +2})); + expect(eq(dataset.signal_ranges[0], gr::Range{-2, +2})); expect(eq(dataset.timing_events[0].size(), 1UZ)); expect(eq(dataset.timing_events[0][0].first, 0)); receivedTags.insert(receivedTags.end(), dataset.timing_events[0].begin(), dataset.timing_events[0].end()); @@ -626,7 +626,7 @@ const boost::ut::suite DataSinkTests = [] { expect(eq(dataset.timing_events.size(), 1u)); expect(eq(dataset.signal_names[0], "test signal"s)); expect(eq(dataset.signal_units[0], "none"s)); - expect(eq(dataset.signal_ranges[0], std::vector{0, kSamples - 1})); + expect(eq(dataset.signal_ranges[0], gr::Range{0, kSamples - 1})); expect(eq(dataset.timing_events[0].size(), 1u)); expect(eq(dataset.timing_events[0][0].first, -5000)); receivedData.insert(receivedData.end(), dataset.signal_values.begin(), dataset.signal_values.end()); diff --git a/blocks/fourier/include/gnuradio-4.0/fourier/fft.hpp b/blocks/fourier/include/gnuradio-4.0/fourier/fft.hpp index 28519170..77e7f0a1 100644 --- a/blocks/fourier/include/gnuradio-4.0/fourier/fft.hpp +++ b/blocks/fourier/include/gnuradio-4.0/fourier/fft.hpp @@ -202,7 +202,6 @@ On the choice of window (mathematically aka. apodisation) functions ds.signal_ranges[i] = {*mm.first, *mm.second}; } - ds.signal_errors = {}; ds.meta_information = {{{"sample_rate", sample_rate}, {"signal_name", signal_name}, {"signal_unit", signal_unit}, {"signal_min", signal_min}, {"signal_max", signal_max}, // {"fft_size", fftSize}, {"window", window}, {"output_in_db", outputInDb}, {"output_in_deg", outputInDeg}, {"unwrap_phase", unwrapPhase}, // {"input_chunk_size", this->input_chunk_size}, {"output_chunk_size", this->output_chunk_size}, {"stride", this->stride}}}; diff --git a/blocks/fourier/test/qa_fourier.cpp b/blocks/fourier/test/qa_fourier.cpp index 0887570a..db2df078 100644 --- a/blocks/fourier/test/qa_fourier.cpp +++ b/blocks/fourier/test/qa_fourier.cpp @@ -71,8 +71,8 @@ void equalDataset(const gr::blocks::fft::FFT>& fftBlock, const for (std::size_t i = 0U; i < 5; i++) { const auto mm = std::minmax_element(std::next(ds1.signal_values.begin(), static_cast(i * N)), std::next(ds1.signal_values.begin(), static_cast((i + 1U) * N))); - expect(approx(*mm.first, ds1.signal_ranges[i][0], tolerance)); - expect(approx(*mm.second, ds1.signal_ranges[i][1], tolerance)); + expect(approx(*mm.first, ds1.signal_ranges[i].min, tolerance)); + expect(approx(*mm.second, ds1.signal_ranges[i].max, tolerance)); } } diff --git a/core/include/gnuradio-4.0/DataSet.hpp b/core/include/gnuradio-4.0/DataSet.hpp index 9c23ed3e..cc4548eb 100644 --- a/core/include/gnuradio-4.0/DataSet.hpp +++ b/core/include/gnuradio-4.0/DataSet.hpp @@ -1,21 +1,33 @@ #ifndef GNURADIO_DATASET_HPP #define GNURADIO_DATASET_HPP -#include "Tag.hpp" #include #include -#include #include #include #include #include +#include + +#include "Message.hpp" +#include "Tag.hpp" + namespace gr { struct LayoutRight {}; struct LayoutLeft {}; +template +struct Range { + T min = 0; + T max = 0; + GR_MAKE_REFLECTABLE(Range, min, max); + + auto operator<=>(const Range& other) const = default; +}; + /** * @brief a concept that describes a Packet, which is a subset of the DataSet struct. */ @@ -39,7 +51,6 @@ concept TensorLike = PacketLike && requires(T t, const std::size_t n_items) { requires std::is_same_v>; requires std::is_same_v; requires std::is_same_v>; - requires std::is_same_v>; requires std::is_same_v>; }; @@ -68,12 +79,11 @@ concept DataSetLike = TensorLike && requires(T t, const std::size_t n_items) requires std::is_same_v>; requires std::is_same_v>; requires std::is_same_v>; - requires std::is_same_v>; - requires std::is_same_v>>; + requires std::is_same_v>>; // meta data requires std::is_same_v>; - requires std::is_same_v>>>; + requires std::is_same_v>>>; }; template @@ -84,7 +94,7 @@ struct DataSet { std::int64_t timestamp = 0; // UTC timestamp [ns] // axis layout: - std::vector axis_names{}; // e.g. time, frequency, … + std::vector axis_names{}; // axis quantity, e.g. time, frequency, … std::vector axis_units{}; // axis base SI-unit std::vector> axis_values{}; // explicit axis values @@ -93,18 +103,62 @@ struct DataSet { tensor_layout_type layout{}; // row-major, column-major, “special” // signal data storage: - std::vector signal_names{}; // size = extents[0] - std::vector signal_quantities{}; // size = extents[0] - std::vector signal_units{}; // size = extents[0] - std::vector signal_values{}; // size = \PI_i extents[i] - std::vector signal_errors{}; // size = \PI_i extents[i] or '0' if not applicable - std::vector> signal_ranges{}; // [[min_0, max_0], [min_1, max_1], …] used for communicating, for example, HW limits + std::vector signal_names{}; // size = extents[0] + std::vector signal_quantities{}; // size = extents[0] + std::vector signal_units{}; // size = extents[0] + std::vector signal_values{}; // size = \PI_i extents[i] + std::vector> signal_ranges{}; // [[min_0, max_0], [min_1, max_1], …] used for communicating, for example, HW limits // meta data - std::vector meta_information{}; - std::vector>> timing_events{}; - - GR_MAKE_REFLECTABLE(DataSet, timestamp, axis_names, axis_units, axis_values, extents, layout, signal_names, signal_quantities, signal_units, signal_values, signal_errors, signal_ranges, meta_information, timing_events); + std::vector meta_information{}; + std::vector>> timing_events{}; + + GR_MAKE_REFLECTABLE(DataSet, timestamp, axis_names, axis_units, axis_values, extents, layout, signal_names, signal_quantities, signal_units, signal_values, signal_ranges, meta_information, timing_events); + + [[nodiscard]] std::size_t axisCount() const noexcept { return axis_names.size(); } + [[nodiscard]] std::string& axisName(std::size_t i = 0UZ) { return axis_names[_axCheck(i)]; } + [[nodiscard]] std::string_view axisName(std::size_t i = 0UZ) const { return axis_names[_axCheck(i)]; } + [[nodiscard]] std::string& axisUnit(std::size_t i = 0UZ) { return axis_units[_axCheck(i)]; } + [[nodiscard]] std::string_view axisUnit(std::size_t i = 0UZ) const { return axis_units[_axCheck(i)]; } + [[nodiscard]] std::span axisValues(std::size_t i = 0UZ) { return axis_values[_axCheck(i)]; } + [[nodiscard]] std::span axisValues(std::size_t i = 0UZ) const { return axis_values[_axCheck(i)]; } + + [[nodiscard]] constexpr std::size_t size() const noexcept { return signal_names.size(); } + [[nodiscard]] std::string& signalName(std::size_t i = 0UZ) { return signal_names[_idxCheck(i)]; } + [[nodiscard]] std::string_view signalName(std::size_t i = 0UZ) const { return signal_names[_idxCheck(i)]; } + [[nodiscard]] std::string& signalQuantity(std::size_t i = 0UZ) { return signal_quantities[_idxCheck(i)]; } + [[nodiscard]] std::string_view signalQuantity(std::size_t i = 0UZ) const { return signal_quantities[_idxCheck(i)]; } + [[nodiscard]] std::string& signalUnit(std::size_t i = 0UZ) { return signal_units[_idxCheck(i)]; } + [[nodiscard]] std::string_view signalUnit(std::size_t i = 0UZ) const { return signal_units[_idxCheck(i)]; } + [[nodiscard]] std::span signalValues(std::size_t i = 0UZ) { return {std::next(signal_values.data(), _idxCheckS(i) * _valsPerSigS()), _valsPerSig()}; } + [[nodiscard]] std::span signalValues(std::size_t i = 0UZ) const { return {std::next(signal_values.data(), _idxCheckS(i) * _valsPerSigS()), _valsPerSig()}; } + [[nodiscard]] Range& signalRange(std::size_t i = 0UZ) { return signal_ranges[_idxCheck(i)]; } + [[nodiscard]] const Range& signalRange(std::size_t i = 0UZ) const { return signal_ranges[_idxCheck(i)]; } + +private: + [[nodiscard]] std::size_t _axCheck(std::size_t i, std::source_location loc = std::source_location::current()) const { + if (i >= axis_names.size()) { + throw gr::exception(fmt::format("{} axis out of range: i={} >= [0, {}]", loc.function_name(), i, axis_names.size()), loc); + } + return i; + } + + [[nodiscard]] std::size_t _idxCheck(std::size_t i, std::source_location location = std::source_location::current()) const { + if (i >= size()) { + throw gr::exception(fmt::format("{} out of range: i={} >= [0, {}]", location.function_name(), i, size()), location); + } + return i; + } + + [[nodiscard]] std::ptrdiff_t _idxCheckS(std::size_t i, std::source_location location = std::source_location::current()) const { + if (i >= size()) { + throw gr::exception(fmt::format("{} out of range: i={} >= [0, {}]", location.function_name(), i, size()), location); + } + return static_cast(i); + } + + [[nodiscard]] std::size_t _valsPerSig() const noexcept { return size() == 0U ? 0U : signal_values.size() / size(); } + [[nodiscard]] std::ptrdiff_t _valsPerSigS() const noexcept { return static_cast(size() == 0U ? 0U : signal_values.size() / size()); } }; static_assert(DataSetLike>, "DataSet concept conformity"); @@ -122,12 +176,11 @@ struct Tensor { tensor_layout_type layout{}; // row-major, column-major, “special” std::vector signal_values{}; // size = \PI_i extents[i] - std::vector signal_errors{}; // size = \PI_i extents[i] or '0' if not applicable // meta data std::vector meta_information{}; - GR_MAKE_REFLECTABLE(Tensor, timestamp, extents, layout, signal_values, signal_errors, meta_information); + GR_MAKE_REFLECTABLE(Tensor, timestamp, extents, layout, signal_values, meta_information); }; static_assert(TensorLike>, "Tensor concept conformity"); diff --git a/core/include/gnuradio-4.0/Message.hpp b/core/include/gnuradio-4.0/Message.hpp index 504561c3..1692f434 100644 --- a/core/include/gnuradio-4.0/Message.hpp +++ b/core/include/gnuradio-4.0/Message.hpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include diff --git a/core/test/CMakeLists.txt b/core/test/CMakeLists.txt index 049d67d0..54c11631 100644 --- a/core/test/CMakeLists.txt +++ b/core/test/CMakeLists.txt @@ -49,6 +49,7 @@ endfunction() add_ut_test(qa_buffer) add_ut_test(qa_AtomicBitset) +add_ut_test(qa_DataSet) add_ut_test(qa_DynamicBlock) add_ut_test(qa_DynamicPort) add_ut_test(qa_HierBlock) diff --git a/core/test/qa_DataSet.cpp b/core/test/qa_DataSet.cpp new file mode 100644 index 00000000..243dfee1 --- /dev/null +++ b/core/test/qa_DataSet.cpp @@ -0,0 +1,93 @@ +#include + +#include + +const boost::ut::suite<"DataSet"> _dataSetAPI = [] { + using namespace boost::ut; + using namespace boost::ut::literals; + using namespace std::string_view_literals; + + "DataSet axis + signal access"_test = [] { + gr::DataSet ds; + + ds.axis_names = {"Time", "Frequency"}; + ds.axis_units = {"s", "Hz"}; + ds.axis_values = { + {0.f, 1.f, 2.f, 3.f}, // axis=0 + {100.f, 200.f, 300.f, 400.f} // axis=1 + }; + + ds.signal_names = {"sigA", "sigB"}; + ds.signal_quantities = {"quantityA", "quantityB"}; + ds.signal_units = {"unitA", "unitB"}; + + ds.signal_values = { + 0.f, 1.f, 2.f, 3.f, // data for sigA + 10.f, 11.f, 12.f, 13.f // data for sigB + }; + ds.signal_ranges = { + {0.f, 3.f}, // for sigA + {10.f, 13.f} // for sigB + }; + + expect(eq(ds.axisCount(), 2UZ)); + expect(eq(ds.size(), 2UZ)); + + "non-const usage - axes"_test = [&] { + ds.axisName(0) = "Time_modified"; + expect(eq(ds.axisName(0), "Time_modified"sv)); + + ds.axisUnit(1) = "kHz"; + expect(eq(ds.axisUnit(1), "kHz"sv)); + + auto timeVals = ds.axisValues(0); + expect(eq(timeVals.size(), 4UZ)); + timeVals[2] = 42.f; + expect(eq(ds.axisValues(0)[2], 42.f)); + }; + + "non-const usage - signals"_test = [&] { + ds.signalName(1) = "sigB_mod"; + expect(eq(ds.signalName(1), "sigB_mod"sv)); + + auto sA = ds.signalValues(0); + expect(eq(sA.size(), 4UZ)); + sA[3] = 99.f; + expect(eq(ds.signalValues(0)[3], 99.f)); + + ds.signalRange(1).min = 9.f; // set min + ds.signalRange(1).max = 15.f; // set max + expect(eq(ds.signalRange(1).min, 9.f)); + expect(eq(ds.signalRange(1).max, 15.f)); + }; + + "const usage - reading axes"_test = [&] { + const auto& constDs = ds; // forces use of const methods + + expect(eq(constDs.axisCount(), 2UZ)); + expect(eq(constDs.axisName(0), "Time_modified"sv)); + expect(eq(constDs.axisUnit(1), "kHz"sv)); + expect(eq(constDs.axisValues(0)[2], 42.f)); + }; + + "const usage - reading signals"_test = [&] { + const auto& constDs = ds; // forces use of const methods + expect(eq(constDs.signalName(1), "sigB_mod"sv)); + expect(eq(constDs.signalValues(0)[3], 99.f)); + expect(eq(constDs.signalRange(1).max, 15.f)); + }; + + "out-of-range checks - axes"_test = [&] { + expect(throws([&] { std::ignore = ds.axisName(99); })); + expect(throws([&] { std::ignore = ds.axisValues(2); })); + }; + + "out-of-range checks - signals"_test = [&] { + expect(throws([&] { std::ignore = ds.signalName(99); })); + expect(throws([&] { std::ignore = ds.signalValues(99); })); + expect(throws([&] { std::ignore = ds.signalRange(99); })); + }; + }; +}; + +int main() { /* tests are statically executed */ } diff --git a/meta/include/gnuradio-4.0/meta/UncertainValue.hpp b/meta/include/gnuradio-4.0/meta/UncertainValue.hpp index a01e1d59..d1df17f9 100644 --- a/meta/include/gnuradio-4.0/meta/UncertainValue.hpp +++ b/meta/include/gnuradio-4.0/meta/UncertainValue.hpp @@ -79,7 +79,7 @@ concept UncertainValueLike = gr::meta::is_instantiation_of; template requires arithmetic_or_complex_like> -[[nodiscard]] inline constexpr auto value(const T& val) noexcept { +[[nodiscard]] constexpr auto value(const T& val) noexcept { if constexpr (UncertainValueLike) { return val.value; } else { @@ -89,7 +89,7 @@ requires arithmetic_or_complex_like> template requires arithmetic_or_complex_like> -[[nodiscard]] inline constexpr auto uncertainty(const T& val) noexcept { +[[nodiscard]] constexpr auto uncertainty(const T& val) noexcept { if constexpr (UncertainValueLike) { return val.uncertainty; } else { @@ -119,7 +119,7 @@ using UncertainValueType_t = detail::UncertainValueValueType::type; template, typename ValueTypeU = UncertainValueType_t> requires(UncertainValueLike || UncertainValueLike) && std::is_same_v, meta::fundamental_base_value_type_t> -[[nodiscard]] inline constexpr auto operator+(const T& lhs, const U& rhs) noexcept { +[[nodiscard]] constexpr auto operator+(const T& lhs, const U& rhs) noexcept { if constexpr (UncertainValueLike && UncertainValueLike) { using ResultType = decltype(lhs.value + rhs.value); if constexpr (meta::complex_like || meta::complex_like) { @@ -141,13 +141,13 @@ requires(UncertainValueLike || UncertainValueLike) && std::is_same_v -inline constexpr T& operator+=(T& lhs, const U& rhs) noexcept { +constexpr T& operator+=(T& lhs, const U& rhs) noexcept { lhs = lhs + rhs; return lhs; } template> -inline constexpr T operator+(const T& val) { +constexpr T operator+(const T& val) { if constexpr (meta::complex_like) { return val; } else { @@ -157,7 +157,7 @@ inline constexpr T operator+(const T& val) { template, typename ValueTypeU = UncertainValueType_t> requires(UncertainValueLike || UncertainValueLike) && std::is_same_v, meta::fundamental_base_value_type_t> -[[nodiscard]] inline constexpr auto operator-(const T& lhs, const U& rhs) noexcept { +[[nodiscard]] constexpr auto operator-(const T& lhs, const U& rhs) noexcept { if constexpr (UncertainValueLike && UncertainValueLike) { using ResultType = decltype(lhs.value - rhs.value); if constexpr (meta::complex_like || meta::complex_like) { @@ -178,19 +178,19 @@ requires(UncertainValueLike || UncertainValueLike) && std::is_same_v -inline constexpr T& operator-=(T& lhs, const U& rhs) noexcept { +constexpr T& operator-=(T& lhs, const U& rhs) noexcept { lhs = lhs - rhs; return lhs; } template -inline constexpr T operator-(const T& val) { +constexpr T operator-(const T& val) { return {-val.value, val.uncertainty}; } template, typename ValueTypeU = UncertainValueType_t> requires(UncertainValueLike || UncertainValueLike) && std::is_same_v, meta::fundamental_base_value_type_t> -[[nodiscard]] inline constexpr auto operator*(const T& lhs, const U& rhs) noexcept { +[[nodiscard]] constexpr auto operator*(const T& lhs, const U& rhs) noexcept { if constexpr (UncertainValueLike && UncertainValueLike) { using ResultType = decltype(lhs.value * rhs.value); if constexpr (meta::complex_like || meta::complex_like) { @@ -212,14 +212,14 @@ requires(UncertainValueLike || UncertainValueLike) && std::is_same_v -inline constexpr T& operator*=(T& lhs, const U& rhs) noexcept { +constexpr T& operator*=(T& lhs, const U& rhs) noexcept { lhs = lhs * rhs; return lhs; } template, typename ValueTypeU = UncertainValueType_t> requires(UncertainValueLike || UncertainValueLike) && std::is_same_v, meta::fundamental_base_value_type_t> -[[nodiscard]] inline constexpr auto operator/(const T& lhs, const U& rhs) noexcept { +[[nodiscard]] constexpr auto operator/(const T& lhs, const U& rhs) noexcept { if constexpr (UncertainValueLike && UncertainValueLike) { using ResultType = decltype(lhs.value * rhs.value); if constexpr (meta::complex_like || meta::complex_like) { @@ -252,7 +252,7 @@ requires(UncertainValueLike || UncertainValueLike) && std::is_same_v -inline constexpr T& operator/=(T& lhs, const U& rhs) noexcept { +constexpr T& operator/=(T& lhs, const U& rhs) noexcept { lhs = lhs / rhs; return lhs; } @@ -261,9 +261,15 @@ inline constexpr T& operator/=(T& lhs, const U& rhs) noexcept { namespace gr::math { +template +requires(std::is_arithmetic_v && std::is_arithmetic_v) +[[nodiscard]] constexpr T pow(const T& base, U exponent) noexcept { + return std::pow(base, exponent); +} + template> requires std::is_same_v, U> || std::integral -[[nodiscard]] inline constexpr T pow(const T& base, U exponent) noexcept { +[[nodiscard]] constexpr T pow(const T& base, U exponent) noexcept { if (base.value == static_cast>(0)) [[unlikely]] { if (exponent == 0) [[unlikely]] { return T{1, 0}; @@ -283,8 +289,8 @@ requires std::is_same_v, U> template, typename ValueTypeU = gr::UncertainValueType_t> requires std::is_same_v, gr::meta::fundamental_base_value_type_t> -[[nodiscard]] inline constexpr T pow(const T& base, const U& exponent) noexcept { - if (base.value == 0.0) [[unlikely]] { +[[nodiscard]] constexpr T pow(const T& base, const U& exponent) noexcept { + if (base.value == ValueTypeT(0)) [[unlikely]] { if (exponent.value == static_cast(0)) [[unlikely]] { return T{1, 0}; } else { @@ -302,7 +308,7 @@ requires std::is_same_v, gr: } template -[[nodiscard]] inline constexpr T sqrt(const T& value) noexcept { +[[nodiscard]] constexpr T sqrt(const T& value) noexcept { if constexpr (gr::UncertainValueLike) { using ValueType = meta::fundamental_base_value_type_t; return gr::math::pow(value, ValueType(0.5)); @@ -312,7 +318,7 @@ template } template -[[nodiscard]] inline constexpr T sin(const T& x) noexcept { +[[nodiscard]] constexpr T sin(const T& x) noexcept { if constexpr (gr::UncertainValueLike) { return T{std::sin(x.value), std::abs(std::cos(x.value) * x.uncertainty)}; } else { @@ -321,7 +327,7 @@ template } template -[[nodiscard]] inline constexpr T cos(const T& x) noexcept { +[[nodiscard]] constexpr T cos(const T& x) noexcept { if constexpr (gr::UncertainValueLike) { return T{std::cos(x.value), std::abs(std::sin(x.value) * x.uncertainty)}; } else { @@ -330,7 +336,7 @@ template } template> -[[nodiscard]] inline constexpr T exp(const T& x) noexcept { +[[nodiscard]] constexpr T exp(const T& x) noexcept { if constexpr (gr::meta::complex_like) { return gr::math::pow(gr::UncertainValue{std::numbers::e_v, static_cast(0)}, x); } else { @@ -338,6 +344,48 @@ template +[[nodiscard]] constexpr bool isfinite(const T& value) noexcept { + if constexpr (gr::UncertainValueLike) { + return std::isfinite(gr::value(value)) && std::isfinite(gr::uncertainty(value)); + } else { + return std::isfinite(value); + } +} + +template +[[nodiscard]] constexpr T log(const T& x) noexcept { + if constexpr (UncertainValueLike) { + using base_t = gr::meta::fundamental_base_value_type_t; + auto val = std::log(gr::value(x)); + if constexpr (gr::meta::complex_like) { + constexpr auto derivative = base_t(1) / x.value; // derivative(log(z)) = 1/z + return T{val, std::abs(derivative) * gr::uncertainty(x)}; + } else { + return T{val, std::abs(gr::uncertainty(x) / gr::value(x))}; // derivative(log(x)) = 1/x + } + } else { + return std::log(x); + } +} + +template +[[nodiscard]] constexpr T log10(const T& x) noexcept { + if constexpr (UncertainValueLike) { + using base_t = gr::meta::fundamental_base_value_type_t; + auto val = std::log10(gr::value(x)); + constexpr auto ln10 = std::numbers::ln10_v; + if constexpr (gr::meta::complex_like) { + constexpr auto derivative = base_t(1) / (x.value * ln10); // derivative(log10(z)) = 1 / (z * ln(10)) + return T{val, std::abs(derivative) * gr::uncertainty(x)}; + } else { + return T{val, std::abs(gr::uncertainty(x) / (gr::value(x) * ln10))}; // derivative(log10(x)) = 1 / (x * ln(10)) + } + } else { + return std::log10(x); + } +} + } // namespace gr::math #endif // GNURADIO_UNCERTAINVALUE_HPP diff --git a/meta/include/gnuradio-4.0/meta/formatter.hpp b/meta/include/gnuradio-4.0/meta/formatter.hpp index 1accee6c..101586e0 100644 --- a/meta/include/gnuradio-4.0/meta/formatter.hpp +++ b/meta/include/gnuradio-4.0/meta/formatter.hpp @@ -7,11 +7,12 @@ #include #include #include -#include -#include #include #include +#include +#include + namespace gr { namespace time { [[nodiscard]] inline std::string getIsoTime() noexcept { @@ -98,6 +99,39 @@ struct fmt::formatter> { } }; +namespace gr { +template +std::ostream& operator<<(std::ostream& os, const T& v) { + return os << fmt::format("{}", v); +} +} // namespace gr + +// DataSet - Range formatter + +namespace gr { +template +struct Range; +} + +namespace fmt { +template +struct formatter> { + constexpr auto parse(format_parse_context& ctx) { return ctx.begin(); } + + template + auto format(const gr::Range& range, FormatContext& ctx) { + return fmt::format_to(ctx.out(), "[min: {}, max: {}]", range.min, range.max); + } +}; +} // namespace fmt + +namespace gr { +template +std::ostream& operator<<(std::ostream& os, const gr::Range& v) { + return os << fmt::format("{}", v); +} +} // namespace gr + // pmt formatter namespace gr { diff --git a/meta/test/qa_UncertainValue.cpp b/meta/test/qa_UncertainValue.cpp index 0547971b..a469a231 100644 --- a/meta/test/qa_UncertainValue.cpp +++ b/meta/test/qa_UncertainValue.cpp @@ -20,7 +20,7 @@ template } } // namespace test::detail -const boost::ut::suite uncertainValue = [] { +const boost::ut::suite<"basic UncertainValue tests"> _uncertainValue = [] { using namespace boost::ut; using namespace std::literals::complex_literals; using namespace gr; @@ -707,7 +707,7 @@ const boost::ut::suite uncertainValue = [] { }; }; -const boost::ut::suite uncertainValueTrigonometric = [] { +const boost::ut::suite<"UncertainValue trigonometric functions"> _uncertainValueTrigonometric = [] { using namespace boost::ut; using namespace std::literals::complex_literals; using namespace gr; @@ -763,7 +763,7 @@ const boost::ut::suite uncertainValueTrigonometric = [] { }; }; -const boost::ut::suite uncertainValueExpTests = [] { +const boost::ut::suite<"UncertainValue exponential functions"> uncertainValueExpTests = [] { using namespace boost::ut; using namespace std::literals::complex_literals; using namespace gr; @@ -804,4 +804,72 @@ const boost::ut::suite uncertainValueExpTests = [] { }; }; +const boost::ut::suite<"UncertainValue - additional tests"> _additionalTests = [] { + using namespace boost::ut; + using namespace std::literals::complex_literals; + using namespace gr; + using test::detail::approx; + + "math operations and isfinite"_test = [] { + using ValueType = gr::meta::fundamental_base_value_type_t; + + constexpr static auto initialize = [](typename gr::meta::fundamental_base_value_type_t value, typename gr::meta::fundamental_base_value_type_t uncertainty = {}) -> T { + if constexpr (gr::UncertainValueLike) { + return T{value, uncertainty}; + } else { + return value; + } + }; + + "log"_test = [] { + T x = initialize(ValueType(10), ValueType(0.1)); + T result = gr::math::log(x); + auto expected_value = std::log(ValueType(10)); + auto expected_uncertainty = gr::UncertainValueLike ? ValueType(0.1) / ValueType(10) : ValueType(0); + expect(approx(gr::value(result), expected_value, ValueType(1e-3f))) << "log(x): value"; + if constexpr (gr::UncertainValueLike) { + expect(approx(gr::uncertainty(result), expected_uncertainty, ValueType(1e-3f))) << "log(x): uncertainty"; + } + }; + + "log10"_test = [] { + T x = initialize(ValueType(10), ValueType(0.1)); + T result = gr::math::log10(x); + auto expected_value = std::log10(ValueType(10)); + auto expected_uncertainty = gr::UncertainValueLike ? ValueType(0.1) / (ValueType(10) * std::numbers::ln10_v) : ValueType(0); + expect(approx(gr::value(result), expected_value, ValueType(1e-3f))) << "log10(x): value"; + if constexpr (gr::UncertainValueLike) { + expect(approx(gr::uncertainty(result), expected_uncertainty, ValueType(1e-3f))) << "log10(x): uncertainty"; + } + }; + + "edge case: log and log10 with a zero value"_test = [] { + T x = initialize(ValueType(0), ValueType(0)); + auto result_log = gr::math::log(x); + auto result_log10 = gr::math::log10(x); + expect(std::isinf(gr::value(result_log)) && gr::value(result_log) < 0) << "log(0): value is -inf"; + expect(std::isinf(gr::value(result_log10)) && gr::value(result_log10) < 0) << "log10(0): value is -inf"; + }; + + "edge case: isfinite with finite and infinite inputs"_test = [] { + T x = initialize(ValueType(5), ValueType(0.1)); + expect(gr::math::isfinite(x)) << "isfinite: finite value"; + + if constexpr (gr::UncertainValueLike) { + T y{ValueType(5), std::numeric_limits::infinity()}; + expect(!gr::math::isfinite(y)) << "isfinite: infinite uncertainty"; + + T z{std::numeric_limits::infinity(), ValueType(0.1)}; + expect(!gr::math::isfinite(z)) << "isfinite: infinite value"; + + T w{std::numeric_limits::infinity(), std::numeric_limits::infinity()}; + expect(!gr::math::isfinite(w)) << "isfinite: infinite value and uncertainty"; + } else { + T y = std::numeric_limits::infinity(); + expect(!gr::math::isfinite(y)) << "isfinite: infinite value"; + } + }; + } | std::tuple, gr::UncertainValue>{}; +}; + int main() { /* tests are statically executed */ } diff --git a/meta/test/qa_formatter.cpp b/meta/test/qa_formatter.cpp index 4f137715..e06741e9 100644 --- a/meta/test/qa_formatter.cpp +++ b/meta/test/qa_formatter.cpp @@ -4,7 +4,8 @@ #include #include - +#include +#include #include namespace gr::meta::test { @@ -68,6 +69,10 @@ const boost::ut::suite uncertainValueFormatter = [] { expect(eq("(1.230 ± 0.450)"s, fmt::format("{:1.3f}", UncertainDouble{1.23, 0.45}))); expect(eq("(3.140 ± 0.010)"s, fmt::format("{:1.3f}", UncertainDouble{3.14, 0.01}))); expect(eq("(0.000 ± 0.000)"s, fmt::format("{:1.3f}", UncertainDouble{0, 0}))); + + std::stringstream ss; + ss << UncertainDouble{1.23, 0.45}; + expect(eq("(1.23 ± 0.45)"s, ss.str())); }; }; @@ -110,6 +115,20 @@ const boost::ut::suite expectedFormatter = [] { expect(eq(error, ""s)); }; +const boost::ut::suite<"Range formatter"> _rangeFormatter = [] { + using namespace boost::ut; + using namespace std::literals::string_literals; + + "fmt::formatter>"_test = [] { + gr::Range range{-2, 2}; + expect(eq("[min: -2, max: 2]"s, fmt::format("{}", range))); + }; + "fmt::formatter>"_test = [] { + gr::Range range{-2.5f, 2.5f}; + expect(eq("[min: -2.5, max: 2.5]"s, fmt::format("{}", range))); + }; +}; + } // namespace gr::meta::test int main() { /* tests are statically executed */ }