math.h

Go to the documentation of this file.

#pragma once
 
//https://github.com/stevenlovegrove/Pangolin/issues/352
#ifdef max
#undef max
#endif
 
#ifdef min
#undef min
#endif
 
#include <vector>
#include <array>
#include <cmath>
//#include <math.h>
#include <algorithm>
#include <numeric>
#include <functional>
#include <random>
#include <numbers>
#include <fstream>
 
namespace dnf_composer::tools::math
{
    // https://stackoverflow.com/questions/24518989/how-to-perform-1-dimensional-valid-convolution
    // The conv function implements the standard convolution method,
    // where the output size is the sum of input sizes minus one.
    template<typename T>
    std::vector<T> conv(std::vector<T> const& f, std::vector<T> const& g)
    {
        int const nf = f.size();
        int const ng = g.size();
        int const n = nf + ng - 1;
        std::vector<T> out(n, T());
        for (auto i(0); i < n; ++i) {
            int const jmn = (i >= ng - 1) ? i - (ng - 1) : 0;
            int const jmx = (i < nf - 1) ? i : nf - 1;
            for (auto j(jmn); j <= jmx; ++j) {
                out[i] += (f[j] * g[i - j]);
            }
        }
        return out;
    }
 
    // conv performs standard convolution, while conv_valid performs valid convolution,
    // which is generally more efficient for some applications where zero-padded areas
    // are not desired in the output.
 
    // https://stackoverflow.com/questions/24518989/how-to-perform-1-dimensional-valid-convolution
    // The conv_valid function implements the valid convolution method,
    // where the output size is the absolute difference between the input sizes plus one.
    template<typename T>
    std::vector<T> conv_valid(std::vector<T> const& f, std::vector<T> const& g)
    {
        int const nf = f.size();
        int const ng = g.size();
        std::vector<T> const& min_v = (nf < ng) ? f : g;
        std::vector<T> const& max_v = (nf < ng) ? g : f;
        int const n = std::max(nf, ng) - std::min(nf, ng) + 1;
        std::vector<T> out(n, T());
        for (auto i(0); i < n; ++i) {
            for (int j(min_v.size() - 1), k(i); j >= 0; --j) {
                out[i] += min_v[j] * max_v[k];
                ++k;
            }
        }
        return out;
    }
 
    // ChatGPT 4.0
    template<typename T>
    std::vector<T> conv_same(std::vector<T> const& f, std::vector<T> const& g) {
        int const nf = f.size();
        int const ng = g.size();
        std::vector<T> out(nf, T()); // Output size matches the size of f for 'same' mode
 
        // Calculate padding needed for 'same' mode, adjusted for shift correction
        const int pad = (ng - 1) / 2;
 
        for (int i = 0; i < nf; ++i) {
            T sum = T(); // Initialize the sum for each element of output
            for (int j = 0; j < ng; ++j) {
                int fIndex = i + j - pad; // Adjust index for 'same' mode with shift correction
                if (fIndex >= 0 && fIndex < nf) { // Check boundaries
                    sum += f[fIndex] * g[j];
                }
            }
            out[i] = sum;
        }
        return out;
    }
 
    // In-place variants — write into a caller-supplied buffer (no heap allocation).
    // `out` must be pre-sized to the correct output length before calling.
 
    template<typename T>
    void conv_valid_into(std::vector<T>& out, const std::vector<T>& f, const std::vector<T>& g)
    {
        const int nf = static_cast<int>(f.size());
        const int ng = static_cast<int>(g.size());
        const std::vector<T>& min_v = (nf < ng) ? f : g;
        const std::vector<T>& max_v = (nf < ng) ? g : f;
        const int n = std::max(nf, ng) - std::min(nf, ng) + 1;
        for (int i = 0; i < n; ++i) {
            T acc = T();
            for (int j = static_cast<int>(min_v.size()) - 1, k = i; j >= 0; --j, ++k)
                acc += min_v[j] * max_v[k];
            out[i] = acc;
        }
    }
 
    template<typename T>
    void conv_same_into(std::vector<T>& out, const std::vector<T>& f, const std::vector<T>& g)
    {
        const int nf = static_cast<int>(f.size());
        const int ng = static_cast<int>(g.size());
        const int pad = (ng - 1) / 2;
        for (int i = 0; i < nf; ++i) {
            T acc = T();
            for (int j = 0; j < ng; ++j) {
                const int fIndex = i + j - pad;
                if (fIndex >= 0 && fIndex < nf)
                    acc += f[fIndex] * g[j];
            }
            out[i] = acc;
        }
    }
 
    template<typename T>
    void obtainCircularVector_into(std::vector<T>& out, const std::vector<int>& indices,
                                   const std::vector<T>& contents)
    {
        for (int i = 0; i < static_cast<int>(indices.size()); ++i)
            out[i] = contents[indices[i] - 1];
    }
 
    template<typename T>
    std::vector<T> gaussNorm(const std::vector<int>& rangeX, const T& position, const T& sigma)
    {
        std::vector<T> g(rangeX.size());
        for (int i = 0; i < g.size(); i++)
            g[i] = exp(-0.5 * pow((rangeX[i] - position), 2) / pow(sigma, 2));
 
        if (!g.empty())
        {
            double sumOfG = std::reduce(g.begin(), g.end());
            for (int i = 0; i < g.size(); i++)
                g[i] = g[i] / sumOfG;
        }
 
        return g;
    }
 
    template<typename T>
    std::vector<T> gauss(const std::vector<int>& rangeX, const T& position, const T& sigma)
    {
        std::vector<T> g(rangeX.size());
 
        for (int i = 0; i < g.size(); i++)
            g[i] = exp(-0.5 * pow((rangeX[i] - position), 2) / pow(sigma, 2));
 
        return g;
    }
 
    template<typename T>
    std::vector<T> gauss(int size, const T& sigma, const T& position)
    {
        std::vector<T> g(size);
 
        for (int i = 0; i < g.size(); i++)
        {
            T x = static_cast<T>(i + 1);
            g[i] = exp(-0.5 * pow((x - position), 2) / pow(sigma, 2));
        }
 
        return g;
    }
 
    template<typename T>
    std::vector<T> nonCircularGauss(uint32_t size, const T& sigma, const T& position)
    {
        std::vector<T> g(size);
        std::vector<T> xRange(size);
        std::iota(xRange.begin(), xRange.end(), static_cast<T>(1));
 
        for (uint32_t i = 0; i < size; i++)
        {
            T distance = static_cast<T>(i + 1) - position;
            g[i] = std::exp(-0.5 * std::pow(distance, 2) / std::pow(sigma, 2));
        }
 
        return g;
    }
 
    template<typename T>
    std::vector<T> circularGauss(uint32_t size, const T& sigma, const T& position)
    {
        uint32_t l = size - 2 * 1 + 2;
        uint32_t m = 1;
        T r = position - static_cast<T>(m); // Calculate the shift with fractional part
        T rem = std::fmod(r, static_cast<T>(l)); // Calculate the remainder
        T positionShifted = rem + static_cast<T>(m); // Apply the shifted position
 
        std::vector<T> g(size);
        std::vector<T> xRange(size);
        std::iota(xRange.begin(), xRange.end(), static_cast<T>(1));
 
        std::vector<T> d(size);
        std::vector<T> lMinusd(size);
        std::transform(xRange.begin(), xRange.end(), d.begin(), [&positionShifted, &l](T element) { return std::abs(element - positionShifted); });
        std::transform(d.begin(), d.end(), lMinusd.begin(), [&l](T element) { return -1 * (element - l); });
        for (int i = 0; i < size; i++)
        {
            g[i] = std::exp(-0.5 * std::pow(std::min(d[i], lMinusd[i]), 2) / std::pow(sigma, 2));
        }
 
        return g;
    }
 
    template<typename T>
    std::vector<T> gaussDerivative(const std::vector<int>& rangeX, const T& position, const T& sigma, const T& amplitude)
    {
        std::vector<T> derivative(rangeX.size());
        T variance = sigma * sigma;
 
        for (size_t i = 0; i < derivative.size(); i++)
        {
            T x = rangeX[i];
            derivative[i] = -(x - position) / variance * amplitude * std::exp(-0.5 * std::pow((x - position), 2) / variance);
        }
 
        return derivative;
    }
 
    template<typename T>
    std::vector<T> gaussDerivativeNorm(const std::vector<int>& rangeX, const T& position, const T& sigma, const T& amplitude)
    {
        auto derivative = gaussDerivative(rangeX, position, sigma, amplitude);
 
        // Normalize by the sum of absolute values to preserve sign
        T sumOfAbs = std::accumulate(derivative.begin(), derivative.end(), static_cast<T>(0),
                                     [](T sum, T value) { return sum + std::abs(value); });
 
        if (sumOfAbs > static_cast<T>(0))
        {
            for (T& value : derivative)
                value /= sumOfAbs;
        }
 
        return derivative;
    }
 
    template<typename T>
    std::vector<T> obtainCircularVector(const std::vector<int>& indices, const std::vector<T>& contents)
    {
        std::vector<T> newContents(indices.size());
        for (int i = 0; i < indices.size(); i++)
            newContents[i] = contents[indices[i] - 1];
        return newContents;
    }
 
    template <typename T>
    std::vector<T> sigmoid(const std::vector<T>& x, T beta, T x0)
    {
        std::vector<T> s(x.size());
        for (std::size_t i = 0; i < s.size(); ++i) {
            s[i] = 1 / (1 + std::exp(-beta * (x[i] - x0)));
        }
        return s;
    }
 
    template <typename T>
    std::vector<T> absSigmoid(const std::vector<T>& x, T beta, T x0)
    {
        std::vector<T> s(x.size());
        for (std::size_t i = 0; i < s.size(); ++i) {
            const T diff = x[i] - x0;
            s[i] = static_cast<T>(0.5) * (1 + beta * diff / (1 + beta * std::abs(diff)));
        }
        return s;
    }
 
    template<typename T>
    std::vector<T> heaviside(const std::vector<T>& x, T threshold)
    {
        std::vector<T> h(x.size());
        for (int i = 0; i < static_cast<int>(h.size()); i++)
            h[i] = (x[i] > threshold) ? 1 : 0;
        return h;
    }
 
    template<typename T>
    std::vector<T> sumGauss(const std::vector<T>& gauss1, const std::vector<T>& gauss2)
    {
        std::vector<T> gaussResult(gauss1.size());
        for (int i = 0; i < gauss1.size(); i++)
            gaussResult[i] = gauss1[i] + gauss2[i];
        return gaussResult;
    }
 
    std::array<int, 2> computeKernelRange(double sigma, int cutOfFactor, int fieldSize, bool circular);
    std::vector<int> createExtendedIndex(int fieldSize, const std::array<int, 2>& kernelRange);
 
    std::vector<double> generateNormalVector(int size);
 
    template <typename T>
    std::vector<T> normalize(const std::vector<T>& vector)
    {
        static constexpr T epsilon = 1e-6;
        static constexpr T offset = 1.0;
 
        // Find the minimum and maximum values in the vector
        const T minVal = *std::ranges::min_element(vector.begin(), vector.end()) - epsilon;
 
        std::vector<T> normalizedVector = vector;
        for (T& val : normalizedVector)
            val += minVal;
 
        normalizedVector = math::sigmoid(normalizedVector, 0.2, std::abs(minVal) + offset);
        const T newMinVal = *std::ranges::min_element(normalizedVector.begin(), normalizedVector.end());
 
        for (T& val : normalizedVector)
            val -= newMinVal;
 
        return normalizedVector;
    }
 
    template <typename T>
    std::vector<T> hebbLearningRule(std::vector<T>& weights, const std::vector<T>& input, const std::vector<T>& output, double learningRate)
    {
        if (input.empty() || output.empty())
            throw std::invalid_argument("Input and output vectors cannot be empty");
 
        const size_t inputSize = input.size();
        const size_t outputSize = output.size();
        const size_t expectedSize = inputSize * outputSize;
 
        if (weights.size() != expectedSize)
            throw std::invalid_argument("Weight matrix size mismatch");
 
        for (size_t i = 0; i < inputSize; ++i)
        {
            const T scaledInput = learningRate * input[i];
            const size_t baseIndex = i * outputSize;
 
            for (size_t j = 0; j < outputSize; ++j)
                weights[baseIndex + j] += scaledInput * output[j];
        }
 
        return weights;
    }
 
    template <typename T>
    std::vector<T> ojaLearningRule(std::vector<T>& weights, const std::vector<T>& input, const std::vector<T>& output, double learningRate)
    {
        const int inputSize = input.size();
        const int outputSize = output.size();
 
        for (int i = 0; i < inputSize; i++)
            for (int j = 0; j < outputSize; j++)
            {
                int index = i * outputSize + j; // Compute the index for the flattened matrix
                weights[index] += learningRate * (input[i] * output[j] - output[j] * input[i] * weights[index]);
            }
 
        return weights;
    }
 
    template <typename T>
    std::vector<std::vector<T>> deltaLearningRuleWidrowHoff(std::vector<std::vector<T>>& weights, const std::vector<T>& input,
                                                            const std::vector<T>& actualOutput, const std::vector<T>& targetOutput, double learningRate)
    {
        const int inputSize = input.size();
        const int outputSize = targetOutput.size();
 
        // Calculate the error between the target output and the actual output
        std::vector<T> error(outputSize, 0.0);
        for (size_t j = 0; j < outputSize; ++j) {
            error[j] = targetOutput[j] - actualOutput[j];
        }
 
        // Update the weights based on the error and current activation levels of the fields
        for (size_t i = 0; i < inputSize; ++i) {
            for (size_t j = 0; j < outputSize; ++j) {
                weights[i][j] += learningRate * error[j] * input[i];
            }
        }
 
        return weights;
    }
 
    template <typename T>
    std::vector<std::vector<T>> deltaLearningRuleKroghHertz(std::vector<std::vector<T>>& weights, const std::vector<T>& input,
                                                            const std::vector<T>& targetOutput, const std::vector<T>& actualOutput,
                                                            double learningRate)
    {
        const int inputSize = input.size();
        int outputSize = targetOutput.size();
 
        //std::vector<T> actualOutput(outputSize, 0.0);
        //for (size_t j = 0; j < outputSize; ++j) {
        //  for (size_t i = 0; i < inputSize; ++i) {
        //      actualOutput[j] += input[i] * weights[i][j];
        //  }
        //}
 
        // Calculate the error between the target output and the actual output
        std::vector<T> error(outputSize, 0.0);
        for (size_t j = 0; j < outputSize; ++j) {
            error[j] = targetOutput[j] - actualOutput[j];
        }
 
        // Update the weights based on the error and current activation levels of the fields
        for (size_t i = 0; i < inputSize; ++i) {
            for (size_t j = 0; j < outputSize; ++j) {
                //weights[i][j] += learningRate * (error[j] - learningRate * weights[i][j]) * input[i]; // old
                weights[i][j] += learningRate * (error[j]) * input[i];
            }
        }
 
        return weights;
    }
 
    template <typename T>
    bool compareVectors(const std::vector<T>& vec1, const std::vector<T>& vec2, T threshold) {
        if (vec1.size() != vec2.size()) {
            return false; // Vectors are of different sizes, hence not equal
        }
 
        for (size_t i = 0; i < vec1.size(); ++i) {
            if (std::abs(vec1[i] - vec2[i]) > threshold) {
                return false; // Difference between elements at index i exceeds threshold
            }
        }
 
        return true; // Vectors are equal within the threshold
    }
 
    template <typename T>
    T calculateVectorSum(const std::vector<T>& vec) {
        T result = 0.0;
        for (T value : vec) {
            result += value;
        }
        return result;
    }
 
    template <typename T>
    T calculateVectorAvg(const std::vector<T>& vec) {
        if (vec.empty()) {
            return T(); // Return default value if vector is empty
        }
 
        T sum = T();
        for (const T& value : vec) {
            sum += value;
        }
        return sum / static_cast<T>(vec.size()); // Calculate average
    }
 
    template <typename T>
    T calculateVectorNorm(const std::vector<T>& vec) {
        T sum_of_squares = std::accumulate(vec.begin(), vec.end(), 0.0,
                                           [](T a, T b) { return a + b * b; });
        return std::sqrt(sum_of_squares);
    }
 
    inline double normalize(const double value, const double min, const double max)
    {
        if (value < min) return 0.0;
        if (value > max) return 1.0;
        return (value - min) / (max - min);
    }
 
    inline double gaussian_2d(double x, double y, double mu_x, double mu_y, double sigma_x, double sigma_y, double A) {
        const double exponent = -((std::pow(x - mu_x, 2) / (2 * std::pow(sigma_x, 2))) + (std::pow(y - mu_y, 2) / (2 * std::pow(sigma_y, 2))));
        return A * std::exp(exponent);
    }
 
    inline double circular_gaussian_2d(double x, double y, double mu_x, double mu_y, double sigma, double A) {
        const double exponent = -((std::pow(x - mu_x, 2) + std::pow(y - mu_y, 2)) / (2 * std::pow(sigma, 2)));
        return A * std::exp(exponent);
    }
 
    inline double wrap(double value, double max_value) {
        if (value < 0) return value + max_value;
        if (value >= max_value) return value - max_value;
        return value;
    }
 
    inline double gaussian_2d_periodic(double x, double y, double mu_x, double mu_y, double sigma, double A, double max_x, double max_y) {
        const double xStep = std::min(std::abs(x - mu_x), max_x - std::abs(x - mu_x));
        const double dy = std::min(std::abs(y - mu_y), max_y - std::abs(y - mu_y));
        const double exponent = -((std::pow(xStep, 2) + std::pow(dy, 2)) / (2 * std::pow(sigma, 2)));
        return A * std::exp(exponent);
    }
 
    template <typename T>
    std::vector<T> flattenMatrix(const std::vector<std::vector<T>>& matrix)
    {
        const int rows = matrix.size();
        const int cols = matrix[0].size();
        std::vector<T> flat_matrix(rows * cols);
        for (int i = 0; i < rows; ++i)
            for (int j = 0; j < cols; ++j)
                flat_matrix[i * cols + j] = static_cast<T>(matrix[i][j]);
        return flat_matrix;
    }
 
    // Linear interpolation resampling from input.size() to outputSize.
    template<typename T>
    std::vector<T> resample(const std::vector<T>& input, int outputSize)
    {
        if (input.empty() || outputSize <= 0) return {};
        const int N = static_cast<int>(input.size());
        if (N == outputSize) return input;
        if (outputSize == 1) return { input[N / 2] };
        std::vector<T> out(outputSize);
        for (int i = 0; i < outputSize; ++i)
        {
            const double pos = static_cast<double>(i) * (N - 1) / (outputSize - 1);
            const int lo = static_cast<int>(pos);
            const int hi = std::min(lo + 1, N - 1);
            const double t = pos - lo;
            out[i] = input[lo] * (1.0 - t) + input[hi] * t;
        }
        return out;
    }
 
    // In-place variant: writes into a pre-sized output buffer to avoid heap allocation.
    template<typename T>
    void resampleInto(const std::vector<T>& input, std::vector<T>& output)
    {
        const int N = static_cast<int>(input.size());
        const int outputSize = static_cast<int>(output.size());
        if (input.empty() || outputSize <= 0) return;
        if (N == outputSize) { std::copy(input.begin(), input.end(), output.begin()); return; }
        if (outputSize == 1) { output[0] = input[N / 2]; return; }
        for (int i = 0; i < outputSize; ++i)
        {
            const double pos = static_cast<double>(i) * (N - 1) / (outputSize - 1);
            const int lo = static_cast<int>(pos);
            const int hi = std::min(lo + 1, N - 1);
            const double t = pos - lo;
            output[i] = input[lo] * (1.0 - t) + input[hi] * t;
        }
    }
 
    template<typename T>
    void resampleNearestInto(const std::vector<T>& input, std::vector<T>& output)
    {
        const int N = static_cast<int>(input.size());
        const int M = static_cast<int>(output.size());
        if (input.empty() || M <= 0) return;
        if (N == M) { std::copy(input.begin(), input.end(), output.begin()); return; }
        if (M == 1) { output[0] = input[N / 2]; return; }
        for (int i = 0; i < M; ++i)
        {
            const double pos = static_cast<double>(i) * (N - 1) / (M - 1);
            output[i] = input[static_cast<int>(std::round(pos))];
        }
    }
 
    // Catmull-Rom cubic spline resampling.
    template<typename T>
    void resampleCubicInto(const std::vector<T>& input, std::vector<T>& output)
    {
        const int N = static_cast<int>(input.size());
        const int M = static_cast<int>(output.size());
        if (input.empty() || M <= 0) return;
        if (N == M) { std::copy(input.begin(), input.end(), output.begin()); return; }
        if (M == 1) { output[0] = input[N / 2]; return; }
        const auto clamp = [N](int idx) { return std::max(0, std::min(idx, N - 1)); };
        for (int i = 0; i < M; ++i)
        {
            const double pos = static_cast<double>(i) * (N - 1) / (M - 1);
            const int lo = static_cast<int>(pos);
            const double t = pos - lo;
            const double p0 = static_cast<double>(input[clamp(lo - 1)]);
            const double p1 = static_cast<double>(input[clamp(lo)]);
            const double p2 = static_cast<double>(input[clamp(lo + 1)]);
            const double p3 = static_cast<double>(input[clamp(lo + 2)]);
            output[i] = static_cast<T>(0.5 * (2.0 * p1 + (-p0 + p2) * t +
                (2.0 * p0 - 5.0 * p1 + 4.0 * p2 - p3) * t * t +
                (-p0 + 3.0 * p1 - 3.0 * p2 + p3) * t * t * t));
        }
    }
 
    // Separable 2D convolution.
    // field is stored y-major: field[y * size_x + x].
    // Applies kernel_x along the x-axis and kernel_y along the y-axis using two 1D passes.
    // extIndex_x / extIndex_y are the circular extension indices produced by createExtendedIndex
    // (pass empty vectors for non-circular mode).
    template<typename T>
    std::vector<T> conv2d_separable(
        const std::vector<T>& field,
        const std::vector<T>& kernel_x,
        const std::vector<T>& kernel_y,
        int size_x, int size_y,
        const std::vector<int>& extIndex_x,
        const std::vector<int>& extIndex_y)
    {
        const bool circular_x = !extIndex_x.empty();
        const bool circular_y = !extIndex_y.empty();
 
        // x-pass: convolve each row (fixed y) along the x-axis with kernel_x
        std::vector<T> tmp(size_x * size_y, T());
        for (int y = 0; y < size_y; ++y)
        {
            std::vector<T> row(size_x);
            for (int x = 0; x < size_x; ++x)
                row[x] = field[y * size_x + x];
 
            std::vector<T> convRow;
            if (circular_x)
            {
                const auto extRow = obtainCircularVector(extIndex_x, row);
                convRow = conv_valid(extRow, kernel_x);
            }
            else
            {
                convRow = conv_same(row, kernel_x);
            }
            for (int x = 0; x < size_x; ++x)
                tmp[y * size_x + x] = convRow[x];
        }
 
        // y-pass: convolve each column (fixed x) along the y-axis with kernel_y
        std::vector<T> result(size_x * size_y, T());
        for (int x = 0; x < size_x; ++x)
        {
            std::vector<T> col(size_y);
            for (int y = 0; y < size_y; ++y)
                col[y] = tmp[y * size_x + x];
 
            std::vector<T> convCol;
            if (circular_y)
            {
                const auto extCol = obtainCircularVector(extIndex_y, col);
                convCol = conv_valid(extCol, kernel_y);
            }
            else
            {
                convCol = conv_same(col, kernel_y);
            }
            for (int y = 0; y < size_y; ++y)
                result[y * size_x + x] = convCol[y];
        }
 
        return result;
    }
 
    // In-place 2D separable convolution — writes into caller-supplied buffers.
    // `out` and `tmp` must be pre-sized to size_x * size_y before calling.
    template<typename T>
    void conv2d_separable_into(
        std::vector<T>& out,
        std::vector<T>& tmp,
        const std::vector<T>& field,
        const std::vector<T>& kernel_x,
        const std::vector<T>& kernel_y,
        int size_x, int size_y,
        const std::vector<int>& extIndex_x,
        const std::vector<int>& extIndex_y)
    {
        const bool circular_x = !extIndex_x.empty();
        const bool circular_y = !extIndex_y.empty();
 
        std::vector<T> row(size_x);
        std::vector<T> col(size_y);
        std::vector<T> convRow(size_x);
        std::vector<T> extRow(circular_x ? extIndex_x.size() : 0);
        std::vector<T> convCol(size_y);
        std::vector<T> extCol(circular_y ? extIndex_y.size() : 0);
 
        // x-pass: convolve each row (fixed y) with kernel_x
        for (int y = 0; y < size_y; ++y)
        {
            for (int x = 0; x < size_x; ++x)
                row[x] = field[y * size_x + x];
 
            if (circular_x)
            {
                obtainCircularVector_into(extRow, extIndex_x, row);
                conv_valid_into(convRow, extRow, kernel_x);
            }
            else
            {
                conv_same_into(convRow, row, kernel_x);
            }
            for (int x = 0; x < size_x; ++x)
                tmp[y * size_x + x] = convRow[x];
        }
 
        // y-pass: convolve each column (fixed x) with kernel_y
        for (int x = 0; x < size_x; ++x)
        {
            for (int y = 0; y < size_y; ++y)
                col[y] = tmp[y * size_x + x];
 
            if (circular_y)
            {
                obtainCircularVector_into(extCol, extIndex_y, col);
                conv_valid_into(convCol, extCol, kernel_y);
            }
            else
            {
                conv_same_into(convCol, col, kernel_y);
            }
            for (int y = 0; y < size_y; ++y)
                out[y * size_x + x] = convCol[y];
        }
    }
 
    // Reduction operation used when collapsing one axis of a 2D buffer.
    enum class ReduceOp { SUM, AVERAGE, MAXIMUM, MINIMUM };
 
    // Collapse a y-major 2D buffer (field[y * size_x + x]) along one axis into a 1D
    // result, writing into a pre-sized output buffer.
    // keepX == true:  output has size_x entries; each out[x] reduces over all y.
    // keepX == false: output has size_y entries; each out[y] reduces over all x.
    template<typename T>
    void reduce2DAxis_into(std::vector<T>& out, const std::vector<T>& field,
        int size_x, int size_y, bool keepX, ReduceOp op)
    {
        // Reject non-positive dimensions or a field smaller than size_x*size_y
        // (malformed dimensions / JSON) before sizing or indexing: produce an
        // empty result rather than over-allocating or reading out of bounds.
        if (size_x <= 0 || size_y <= 0 ||
            field.size() < static_cast<std::size_t>(size_x) * static_cast<std::size_t>(size_y))
        {
            out.clear();
            return;
        }
 
        const int outSize = keepX ? size_x : size_y;
        const int reduceCount = keepX ? size_y : size_x;
        if (out.size() != static_cast<std::size_t>(outSize)) out.assign(outSize, T());
 
        for (int o = 0; o < outSize; ++o)
        {
            // First sample of the reduced line.
            auto sampleAt = [&](int k) -> T {
                const int x = keepX ? o : k;
                const int y = keepX ? k : o;
                return field[static_cast<std::size_t>(y) * size_x + x];
            };
            T acc = sampleAt(0);
            for (int k = 1; k < reduceCount; ++k)
            {
                const T v = sampleAt(k);
                switch (op)
                {
                case ReduceOp::SUM:
                case ReduceOp::AVERAGE: acc += v; break;
                case ReduceOp::MAXIMUM: acc = std::max(acc, v); break;
                case ReduceOp::MINIMUM: acc = std::min(acc, v); break;
                }
            }
            if (op == ReduceOp::AVERAGE)
                acc /= static_cast<T>(reduceCount);
            out[o] = acc;
        }
    }
 
    // Broadcast a 1D profile into a y-major 2D buffer (out[y * size_x + x]), writing
    // into a pre-sized output buffer.
    // alongX == true:  profile indexes x (size must be size_x); repeated for every y.
    // alongX == false: profile indexes y (size must be size_y); repeated for every x.
    template<typename T>
    void broadcast1DTo2D_into(std::vector<T>& out, const std::vector<T>& profile,
        int size_x, int size_y, bool alongX)
    {
        // Reject non-positive dimensions (malformed dimensions / JSON) before sizing:
        // a negative size would wrap to an enormous allocation.
        if (size_x <= 0 || size_y <= 0) { out.clear(); return; }
 
        const std::size_t total = static_cast<std::size_t>(size_x) * static_cast<std::size_t>(size_y);
        if (out.size() != total) out.assign(total, T());
        if (profile.empty()) { std::fill(out.begin(), out.end(), T()); return; }
        const int profileSize = static_cast<int>(profile.size());
        for (int y = 0; y < size_y; ++y)
            for (int x = 0; x < size_x; ++x)
            {
                // Clamp the profile index so a mismatched profile size cannot read OOB.
                const int idx = std::min(alongX ? x : y, profileSize - 1);
                out[static_cast<std::size_t>(y) * size_x + x] = profile[idx];
            }
    }
}