interpolate.hpp

Go to the documentation of this file.

#ifndef WEIGHTEDLOWESS_INTERPOLATE_HPP
#define WEIGHTEDLOWESS_INTERPOLATE_HPP
 
#include <vector>
#include <cstddef>
#include <utility>
 
#include "sanisizer/sanisizer.hpp"
 
#include "window.hpp"
#include "utils.hpp"
 
namespace WeightedLowess {
 
struct AssignedSegments {
    std::vector<std::size_t> boundaries;
};
 
template<typename Data_>
AssignedSegments assign_to_segments(
    const Data_* const x_fit,
    const PrecomputedWindows<Data_>& windows_fit,
    const std::size_t num_points_out,
    const Data_* const x_out
) {
    const auto& anchors = windows_fit.anchors;
    const auto num_anchors = anchors.size();
 
    AssignedSegments output;
    sanisizer::resize(output.boundaries, num_anchors);
 
    std::size_t counter = 0;
    {
        const auto right = x_fit[anchors[0]];
        while (counter < num_points_out && x_out[counter] < right) {
            ++counter;
        }
    }
 
    for (I<decltype(num_anchors)> i = 1; i < num_anchors; ++i) {
        output.boundaries[i - 1] = counter;
        if (counter == num_points_out) {
            continue;
        }
        const auto right = x_fit[anchors[i]];
        while (counter < num_points_out && x_out[counter] < right) {
            ++counter;
        }
    }
 
    // Anyone equal to the last anchor get assigned to the last segment.
    {
        const auto right = x_fit[anchors[num_anchors - 1]];
        while (counter < num_points_out && x_out[counter] == right) {
            ++counter;
        }
        output.boundaries[num_anchors - 1] = counter;
    }
 
    return output;
}
 
inline std::pair<std::size_t, std::size_t> get_interpolation_boundaries(const AssignedSegments& assigned_out) {
    return std::make_pair(assigned_out.boundaries.front(), assigned_out.boundaries.back());
}
 
template<typename Data_>
void interpolate(
    const Data_* const x_fit,
    const PrecomputedWindows<Data_>& windows_fit,
    const Data_* const fitted_fit,
    const Data_* const x_out,
    const AssignedSegments& assigned_out,
    Data_* const fitted_out,
    int num_threads
) {
    const auto& anchors = windows_fit.anchors;
    const auto num_anchors = anchors.size();
    const auto num_anchors_m1 = num_anchors - 1;
 
    // One would think that we should parallelize across x_out instead of anchors,
    // as this has better worksharing when x_out is not evenly distributed across anchor segments.
    // However, if we did so, we'd have to store the slope and intercept for the anchor segments first,
    // then look up the slope and intercept for each element of x_out.
    // This involves an extra memory access and is not SIMD-able.
    parallelize(num_threads, num_anchors_m1, [&](const int, const I<decltype(num_anchors_m1)> start, const I<decltype(num_anchors_m1)> length) {
        for (I<decltype(start)> s = start, end = start + length; s < end; ++s) {
            const auto run_start = assigned_out.boundaries[s];
            const auto run_end = assigned_out.boundaries[s + 1];
            if (run_start == run_end) {
                continue;
            }
 
            const auto left_anchor = windows_fit.anchors[s];
            const auto right_anchor = windows_fit.anchors[s + 1];
            const Data_ xdiff = x_fit[right_anchor] - x_fit[left_anchor];
            const Data_ ydiff = fitted_fit[right_anchor] - fitted_fit[left_anchor];
            if (xdiff > 0) {
                const Data_ slope = ydiff / xdiff;
                const Data_ intercept = fitted_fit[right_anchor] - slope * x_fit[right_anchor];
                for (auto outpt = run_start; outpt < run_end; ++outpt) {
                    fitted_out[outpt] = slope * x_out[outpt] + intercept; 
                }
            } else {
                // Protect against infinite slopes by just taking the average.
                const Data_ ave = fitted_fit[left_anchor] + ydiff / 2;
                std::fill(fitted_out + run_start, fitted_out + run_end, ave);
            }
        }
    });
}
 
template<typename Data_>
std::pair<std::size_t, std::size_t> interpolate(
    const Data_* const x_fit,
    const PrecomputedWindows<Data_>& windows_fit,
    const Data_* const fitted_fit,
    const std::size_t num_points_out,
    const Data_* const x_out,
    Data_* const fitted_out,
    int num_threads
) {
    const auto assigned_out = assign_to_segments(x_fit, windows_fit, num_points_out, x_out);
    interpolate(x_fit, windows_fit, fitted_fit, x_out, assigned_out, fitted_out, num_threads);
    return get_interpolation_boundaries(assigned_out);
}
 
}
 
#endif